diff --git a/.gitignore b/.gitignore
index 857fe607..531ab990 100644
--- a/.gitignore
+++ b/.gitignore
@@ -3,4 +3,5 @@ Manifest.toml
 *.scss
 *.css
 .vscode
-*style.jl
\ No newline at end of file
+*style.jl
+.DS_Store
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index 4b10e6f3..ce21a38b 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,7 +2,7 @@ name = "Attractors"
 uuid = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 authors = ["George Datseris <datseris.george@gmail.com>", "Kalel Rossi", "Alexandre Wagemakers"]
 repo = "https://github.com/JuliaDynamics/Attractors.jl.git"
-version = "1.13.6"
+version = "1.14.0"
 
 [deps]
 BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"
diff --git a/docs/make.jl b/docs/make.jl
index 5955d1f6..d2195f14 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -15,7 +15,6 @@ pages = [
     "references.md",
 ]
 
-
 import Downloads
 Downloads.download(
     "https://raw.githubusercontent.com/JuliaDynamics/doctheme/master/build_docs_with_style.jl",
@@ -31,5 +30,5 @@ bib = CitationBibliography(
 )
 
 build_docs_with_style(pages, Attractors, DynamicalSystemsBase, StateSpaceSets;
-    expandfirst = ["index.md"], bib,
+    expandfirst = ["index.md"], bib
 )
diff --git a/docs/refs.bib b/docs/refs.bib
index 2a31c152..7812108c 100644
--- a/docs/refs.bib
+++ b/docs/refs.bib
@@ -158,4 +158,65 @@ @article{Halekotte2020
   author = {Lukas Halekotte and Ulrike Feudel},
   title = {Minimal fatal shocks in multistable complex networks},
   journal = {Scientific Reports}
+}
+
+@article{Battelino1988,
+  title={Multiple coexisting attractors, basin boundaries and basic sets},
+  author={Battelino, Peter M and Grebogi, Celso and Ott, Edward and Yorke, James A and Yorke, Ellen D},
+  journal={Physica D: Nonlinear Phenomena},
+  volume={32},
+  number={2},
+  pages={296--305},
+  year={1988},
+  publisher={Elsevier}
+}
+
+@article{Skufca2006,
+  title={Edge of chaos in a parallel shear flow},
+  author={Skufca, Joseph D and Yorke, James A and Eckhardt, Bruno},
+  journal={Physical review letters},
+  volume={96},
+  number={17},
+  pages={174101},
+  year={2006},
+  publisher={APS}
+}
+
+@article{Schneider2008,
+  title={Laminar-turbulent boundary in plane Couette flow},
+  author={Schneider, Tobias M and Gibson, John F and Lagha, Maher and De Lillo, Filippo and Eckhardt, Bruno},
+  journal={Physical Review E},
+  volume={78},
+  number={3},
+  pages={037301},
+  year={2008},
+  publisher={APS}
+}
+
+@article{Wagemakers2020,
+  title={The saddle-straddle method to test for Wada basins},
+  author={Wagemakers, Alexandre and Daza, Alvar and Sanju{\'a}n, Miguel AF},
+  journal={Communications in Nonlinear Science and Numerical Simulation},
+  volume={84},
+  pages={105167},
+  year={2020},
+  publisher={Elsevier}
+}
+
+@article{Lucarini2017,
+  title={Edge states in the climate system: exploring global instabilities and critical transitions},
+  author={Lucarini, Valerio and B{\'o}dai, Tam{\'a}s},
+  journal={Nonlinearity},
+  volume={30},
+  number={7},
+  pages={R32},
+  year={2017},
+  publisher={IOP Publishing}
+}
+
+@article{Mehling2023,
+  title={Limits to predictability of the asymptotic state of the Atlantic Meridional Overturning Circulation in a conceptual climate model},
+  author={Mehling, Oliver and B{\"o}rner, Reyk and Lucarini, Valerio},
+  journal={arXiv preprint arXiv:2308.16251},
+  year={2023}
 }
\ No newline at end of file
diff --git a/docs/src/basins.md b/docs/src/basins.md
index 7b0ed558..ab36844d 100644
--- a/docs/src/basins.md
+++ b/docs/src/basins.md
@@ -29,6 +29,15 @@ basins_fractal_test
 uncertainty_exponent
 ```
 
+## Edge tracking and edge states
+The edge tracking algorithm allows to locate and construct so-called edge states (also referred to as *Melancholia states*) embedded in the basin boundary separating different basins of attraction. These could be saddle points, unstable periodic orbits or chaotic saddles. The general idea is that these sets can be found because they act as attractors when restricting to the basin boundary.
+
+```@docs
+edgetracking
+EdgeTrackingResults
+bisect_to_edge
+```
+
 ## Tipping points
 This page discusses functionality related with tipping points in dynamical systems with known rule. If instead you are interested in identifying tipping points in measured timeseries, have a look at [TransitionIndicators.jl](https://github.com/JuliaDynamics/TransitionIndicators.jl).
 
diff --git a/docs/src/continuation.md b/docs/src/continuation.md
index 7ef2a471..4c64a355 100644
--- a/docs/src/continuation.md
+++ b/docs/src/continuation.md
@@ -24,6 +24,9 @@ RecurrencesFindAndMatch
 
 ```@docs
 match_statespacesets!
+Centroid
+Hausdorff
+StrictlyMinimumDistance
 replacement_map
 set_distance
 setsofsets_distances
diff --git a/docs/src/examples.md b/docs/src/examples.md
index 55da7da0..beacb291 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -321,7 +321,7 @@ fig
 
 To achieve even better results for this kind of problematic systems than with previuosly introduced `Irregular Grids`  we provide a functionality to construct `Subdivision Based Grids` in which
 one can obtain more coarse or dense structure not only along some axis but for a specific regions where the state space flow has
-significantly different speed. [`subdivided_based_grid`](@ref) enables automatic evaluation of velocity vectors for regions of originally user specified
+significantly different speed. [`subdivision_based_grid`](@ref) enables automatic evaluation of velocity vectors for regions of originally user specified
 grid to further treat those areas as having more dense or coarse structure than others.
 
 ```@example MAIN
@@ -536,7 +536,7 @@ plot_basins_curves(aggregated_fractions, prange;
 
 ## Trivial featurizing and grouping for basins fractions
 
-This is a rather trivial example showcasing the usage of [`AttractorsViaFeaturizing`](@ref). Let us use once again the magnetic pendulum example. For it, we have a really good idea of what features will uniquely describe each attractor: the last points of a trajectory (which should be very close to the magnetic the trajectory converged to). To provide this information to the [`AttractorsviaFeaturizing`](@ref) we just create a julia function that returns this last point
+This is a rather trivial example showcasing the usage of [`AttractorsViaFeaturizing`](@ref). Let us use once again the magnetic pendulum example. For it, we have a really good idea of what features will uniquely describe each attractor: the last points of a trajectory (which should be very close to the magnetic the trajectory converged to). To provide this information to the [`AttractorsViaFeaturizing`](@ref) we just create a julia function that returns this last point
 
 ```@example MAIN
 using Attractors
@@ -871,3 +871,92 @@ animate_attractors_via_recurrences(mapper, u0s)
 <source src="https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/videos/attractors/recurrence_algorithm.mp4?raw=true" type="video/mp4">
 </video>
 ```
+
+## Edge tracking
+
+To showcase how to run the [`edgetracking`](@ref) algorithm, let us use it to find the
+saddle point of the bistable FitzHugh-Nagumo (FHN) model, a two-dimensional ODE system
+originally conceived to represent a spiking neuron.
+We define the system in the following form:
+
+```@example MAIN
+function fitzhugh_nagumo(u,p,t)
+    x, y = u
+    eps, beta = p
+    dx = (x - x^3 - y)/eps
+    dy = -beta*y + x
+    return SVector{2}([dx, dy])
+end
+
+params = [0.1, 3.0]
+ds = CoupledODEs(fitzhugh_nagumo, ones(2), params, diffeq=(;alg = Vern9(), reltol=1e-11))
+```
+
+Now, we can use Attractors.jl to compute the fixed points and basins of attraction of the
+FHN model.
+
+```@example MAIN
+xg = yg = range(-1.5, 1.5; length = 201)
+grid = (xg, yg)
+mapper = AttractorsViaRecurrences(ds, grid; sparse=false)
+basins, attractors = basins_of_attraction(mapper)
+
+for i in 1:length(attractors)
+    println(attractors[i][1])
+end
+```
+
+The `basins_of_attraction` function found three fixed points: the two stable nodes of the
+system (labeled A and B) and the saddle point at the origin. The saddle is an unstable
+equilibrium and can therefore not be found by simulation, but we can find it using the
+[`edgetracking`](@ref) algorithm. For illustration, let us initialize the algorithm from
+two initial conditions `init1` and `init2` (which must belong to different basins
+of attraction, see figure below).
+
+```julia
+attractors_AB = Dict(1 => attractors[1], 2 => attractors[2])
+init1, init2 = [-1.0, -1.0], [-1.0, 0.2]
+```
+
+Now, we run the edge tracking algorithm:
+
+```@example MAIN
+bisect_thresh, diverge_thresh, Δt, abstol = 1e-3, 2e-3, 1e-5, 1e-3
+et = edgetracking(ds, attractors_AB; u1=init1, u2=init2,
+    bisect_thresh, diverge_thresh, Δt, abstol)
+
+et.edge[end]
+```
+
+The algorithm has converged to the origin (up to the specified accuracy) where the saddle
+is located. The figure below shows how the algorithm has iteratively tracked along the basin
+boundary from the two initial conditions (red points) to the saddle (green square). Points
+of the edge track (orange) at which a re-bisection occured are marked with a white border.
+The figure also depicts two trajectories (blue) intialized on either side of the basin
+boundary at the first bisection point. We see that these trajectories follow the basin
+boundary for a while but then relax to either attractor before reaching the saddle. By
+counteracting the instability of the saddle, the edge tracking algorithm instead allows to
+track the basin boundary all the way to the saddle, or edge state.
+
+```@example MAIN
+traj1 = trajectory(ds, 2, et.track1[et.bisect_idx[1]], Δt=1e-5)
+traj2 = trajectory(ds, 2, et.track2[et.bisect_idx[1]], Δt=1e-5)
+
+fig = Figure()
+ax = Axis(fig[1,1], xlabel="x", ylabel="y")
+heatmap_basins_attractors!(ax, grid, basins, attractors, add_legend=false, labels=Dict(1=>"Attractor A", 2=>"Attractor B", 3=>"Saddle"))
+lines!(ax, traj1[1][:,1], traj1[1][:,2], color=:dodgerblue, linewidth=2, label="Trajectories")
+lines!(ax, traj2[1][:,1], traj2[1][:,2], color=:dodgerblue, linewidth=2)
+lines!(ax, et.edge[:,1], et.edge[:,2], color=:orange, linestyle=:dash)
+scatter!(ax, et.edge[et.bisect_idx,1], et.edge[et.bisect_idx,2], color=:white, markersize=15, marker=:circle, zorder=10)
+scatter!(ax, et.edge[:,1], et.edge[:,2], color=:orange, markersize=11, marker=:circle, zorder=10, label="Edge track")
+scatter!(ax, [-1.0,-1.0], [-1.0, 0.2], color=:red, markersize=15, label="Initial conditions")
+xlims!(ax, -1.2, 1.1); ylims!(ax, -1.3, 0.8)
+axislegend(ax, position=:rb)
+fig        
+```
+
+In this simple two-dimensional model, we could of course have found the saddle directly by
+computing the zeroes of the ODE system. However, the edge tracking algorithm allows finding
+edge states also in high-dimensional and chaotic systems where a simple computation of
+unstable equilibria becomes infeasible.
\ No newline at end of file
diff --git a/docs/src/index.md b/docs/src/index.md
index 63694d75..df849e41 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -10,9 +10,7 @@ using CairoMakie, Attractors
 
 ## Latest news
 
-- Our paper on the global stability analysis framework offered by Attractors.jl ([`continuation`](@ref)) and the novel continuation offered by [`RecurrencesFindAndMatch`](@ref) is published as a _Featured Article_ in Chaos (https://pubs.aip.org/aip/cha/article/33/7/073151/2904709/Framework-for-global-stability-analysis-of) and has been featured in the AIP publishing showcase (https://www.growkudos.com/publications/10.1063%25252F5.0159675/reader)
-- New function [`minimal_fatal_shock`](@ref)
-- New function [`match_continuation!`](@ref) which improves the matching during a continuation process where attractors disappear and reappear.
+- New functions [`edgetracking`](@ref) and [`bisect_to_edge`](@ref) added that implement an **edge tracking algorithm** to find saddles or *edge states* in dynamical systems, also when they are unstable chaotic sets.
 
 ## Outline of Attractors.jl
 
diff --git a/docs/src/visualization.md b/docs/src/visualization.md
index 8a0b8cdb..8cab1359 100644
--- a/docs/src/visualization.md
+++ b/docs/src/visualization.md
@@ -11,7 +11,7 @@ heatmap_basins_attractors(grid, basins, attractors; kwargs...)
 heatmap_basins_attractors!(ax, grid, basins, attractors; kwargs...)
 ```
 
-## Common plotting keywords
+## [Common plotting keywords](@id common_plot_kwargs)
 Common keywords for plotting functions in Attractors.jl are:
 
 - `ukeys`: the basin ids (unique keys, vector of integers) to use. By default all existing keys are used.
diff --git a/src/Attractors.jl b/src/Attractors.jl
index 2c074929..d651743f 100644
--- a/src/Attractors.jl
+++ b/src/Attractors.jl
@@ -18,6 +18,7 @@ include("dict_utils.jl")
 include("mapping/attractor_mapping.jl")
 include("basins/basins.jl")
 include("continuation/basins_fractions_continuation_api.jl")
+include("boundaries/edgetracking.jl")
 include("deprecated.jl")
 include("tipping/tipping.jl")
 
diff --git a/src/boundaries/edgetracking.jl b/src/boundaries/edgetracking.jl
new file mode 100644
index 00000000..07625d8f
--- /dev/null
+++ b/src/boundaries/edgetracking.jl
@@ -0,0 +1,314 @@
+export edgetracking, bisect_to_edge, EdgeTrackingResults
+
+"""
+    EdgeTrackingResults(edge, track1, track2, time, bisect_idx)
+Data type that stores output of the [`edgetracking`](@ref) algorithm.
+
+## Fields
+* `edge::StateSpaceSet`: the pseudo-trajectory representing the tracked edge segment
+  (given by the average in state space between `track1` and `track2`)
+* `track1::StateSpaceSet`: the pseudo-trajectory tracking the edge within basin 1
+* `track2::StateSpaceSet`: the pseudo-trajectory tracking the edge within basin 2
+* `time::Vector`: time points of the above `StateSpaceSet`s
+* `bisect_idx::Vector`: indices of `time` at which a re-bisection occurred
+* `success::Bool`: indicates whether the edge tracking has been successful or not
+"""
+struct EdgeTrackingResults{D, T}
+    edge::StateSpaceSet{D, T}
+    track1::StateSpaceSet{D, T}
+    track2::StateSpaceSet{D, T}
+    time::Vector{Float64}
+    bisect_idx::Vector{Int}
+    success::Bool
+end
+
+EdgeTrackingResults(nothing) = EdgeTrackingResults(StateSpaceSet([NaN]),
+    StateSpaceSet([NaN]), StateSpaceSet([NaN]), [NaN], [0], false)
+
+"""
+    edgetracking(ds::DynamicalSystem, attractors::Dict; kwargs...)
+Track along a basin boundary in a dynamical system `ds` with two or more attractors
+in order to find an *edge state*. Results are returned in the form of
+[`EdgeTrackingResults`](@ref), which contains the pseudo-trajectory `edge` representing the track on
+the basin boundary, along with additional output (see below).
+
+The system's `attractors` are specified as a `Dict` of `StateSpaceSet`s, as in
+[`AttractorsViaProximity`](@ref) or the output of [`extract_attractors`](@ref). By default, the
+algorithm is initialized from the first and second attractor in `attractors`. Alternatively,
+the initial states can be set via keyword arguments `u1`, `u2` (see below). Note that the
+two initial states must belong to different basins of attraction.
+
+## Keyword arguments
+* `bisect_thresh = 1e-7`: distance threshold for bisection
+* `diverge_thresh = 1e-6`: distance threshold for parallel integration
+* `u1`: first initial state (defaults to first point in first entry of `attractors`)
+* `u2`: second initial state (defaults to first point in second entry of `attractors`)
+* `maxiter = 100`: maximum number of iterations before the algorithm stops
+* `abstol = 0.0`: distance threshold for convergence of the updated edge state
+* `T_transient = 0.0`: transient time before the algorithm starts saving the edge track
+* `tmax = Inf`: maximum integration time of parallel trajectories until re-bisection 
+* `Δt = 0.01`: time step passed to [`step!`](@ref) when evolving the two trajectories
+* `ϵ_mapper = nothing`: `ϵ` parameter in [`AttractorsViaProximity`](@ref)
+* `show_progress = true`: if true, shows progress bar and information while running
+* `verbose = true`: if false, silences print output and warnings while running
+* `kwargs...`: additional keyword arguments to be passed to [`AttractorsViaProximity`](@ref)
+
+## Description
+The edge tracking algorithm is a numerical method to find
+an *edge state* or (possibly chaotic) saddle on the boundary between two basins of
+attraction. Introduced by [Battelino1988](@cite) and further described by
+[Skufca2006](@cite), the
+algorithm has been applied to, e.g., the laminar-turbulent boundary in plane Couette
+flow [Schneider2008](@cite), Wada basins [Wagemakers2020](@cite), as well as Melancholia
+states in conceptual [Mehling2023](@cite) and intermediate-complexity [Lucarini2017](@cite) 
+climate models. 
+Relying only on forward integration of the system, it works even in
+high-dimensional systems with complicated fractal basin boundary structures.
+
+The algorithm consists of two main steps: bisection and tracking. First, it iteratively 
+bisects along a straight line in state space between the intial states `u1` and `u2` to find
+the separating basin boundary. The bisection stops when the two updated states are less than
+`bisect_thresh` (Euclidean distance in state space) apart from each other.
+Next, a `ParallelDynamicalSystem` is initialized
+from these two updated states and integrated forward until the two trajectories diverge
+from each other by more than `diverge_thresh` (Euclidean distance). The two final states of
+the parallel integration are then used as new states `u1` and `u2` for a new bisection, and 
+so on, until a stopping criterion is fulfilled. 
+
+Two stopping criteria are implemented via the keyword arguments `maxiter` and `abstol`.
+Either the algorithm stops when the number of iterations reaches `maxiter`, or when the
+state space position of the updated edge point changes by less than `abstol` (in
+Euclidean distance) compared to the previous iteration. Convergence below `abstol` happens
+after sufficient iterations if the edge state is a saddle point. However, the edge state
+may also be an unstable limit cycle or a chaotic saddle. In these cases, the algorithm will
+never actually converge to a point but (after a transient period) continue populating the
+set constituting the edge state by tracking along it.
+
+A central idea behind this algorithm is that basin boundaries are typically the stable
+manifolds of unstable sets, namely edge states or saddles. The flow along the basin boundary 
+will thus lead to these sets, and the iterative bisection neutralizes the unstable
+direction of the flow away from the basin boundary. If the system possesses multiple edge 
+states, the algorithm will find one of them depending on where the initial bisection locates
+the boundary.
+
+## Output
+
+Returns a data type [`EdgeTrackingResults`](@ref) containing the results.
+
+Sometimes, the AttractorMapper used in the algorithm may erroneously identify both states
+`u1` and `u2` with the same basin of attraction due to being very close to the basin
+boundary. If this happens, a warning is raised and `EdgeTrackingResults.success = false`.
+"""
+function edgetracking(ds::DynamicalSystem, attractors::Dict;
+    bisect_thresh=1e-6,
+    diverge_thresh=1e-5,
+    u1=collect(values(attractors))[1][1],
+    u2=collect(values(attractors))[2][1],
+    maxiter=100,
+    abstol=0.0,
+    T_transient=0.0,
+    tmax=Inf,
+    Δt=0.01,
+    ϵ_mapper=nothing,
+    show_progress=true,
+    verbose=true,
+    kwargs...)
+    
+    pds = ParallelDynamicalSystem(ds, [u1, u2])
+    mapper = AttractorsViaProximity(ds, attractors, ϵ_mapper; kwargs...)
+    
+    edgetracking(pds, mapper;
+        bisect_thresh, diverge_thresh, maxiter, abstol, T_transient, Δt, tmax,
+        show_progress, verbose)
+end
+
+"""
+    edgetracking(pds::ParallelDynamicalSystem, mapper::AttractorMapper; kwargs...)
+Low-level function for running the edge tracking algorithm, see [`edgetracking`](@ref)
+for a description, keyword arguments and output type.
+
+`pds` is a `ParallelDynamicalSystem` with two states. The `mapper` must be an
+`AttractorMapper` of subtype `AttractorsViaProximity` or `AttractorsViaRecurrences`.
+"""
+function edgetracking(pds::ParallelDynamicalSystem, mapper::AttractorMapper;
+    bisect_thresh=1e-6,
+    diverge_thresh=1e-5,
+    maxiter=100,
+    abstol=0.0,
+    T_transient=0.0,
+    Δt=0.01,
+    tmax=Inf,
+    show_progress=true,
+    verbose=true)
+    
+    if bisect_thresh >= diverge_thresh
+        error("diverge_thresh must be larger than bisect_thresh.")
+    end
+
+    # initial bisection
+    u1, u2, success = bisect_to_edge(pds, mapper; bisect_thresh, verbose)
+    if !success
+        return EdgeTrackingResults(nothing)
+    end
+    edgestate = (u1 + u2)/2
+    track1, track2 = [u1], [u2]
+    time, bisect_idx = Float64[], Int[1]
+    progress = ProgressMeter.Progress(maxiter; desc = "Running edge tracking algorithm",
+        enabled = show_progress)
+    
+    # edge track iteration loop
+    displacement, counter, T = Inf, 1, 0.0
+    while (displacement > abstol) && (maxiter > counter)
+        t = 0
+        set_state!(pds, u1, 1)
+        set_state!(pds, u2, 2)
+        distance = diffnorm(pds)
+        # forward integration loop
+        while (distance < diverge_thresh) && (t < tmax)
+            step!(pds, Δt)
+            distance = diffnorm(pds)
+            t += Δt
+            T += Δt
+            if T >= T_transient
+                push!(track1, current_state(pds, 1))
+                push!(track2, current_state(pds, 2))
+                push!(time, T)
+            end
+        end
+        # re-bisect
+        u1, u2, success = bisect_to_edge(pds, mapper; bisect_thresh, verbose)
+        if ~success
+            track1 = StateSpaceSet(track1)
+            track2 = StateSpaceSet(track2)
+
+            return EdgeTrackingResults(
+                StateSpaceSet((vec(track1) .+ vec(track2))./2),
+                track1, track2, time, bisect_idx, false)
+        end
+        T += Δt
+        if T >= T_transient
+            push!(track1, current_state(pds, 1))
+            push!(track2, current_state(pds, 2))
+            push!(time, T)
+            push!(bisect_idx, length(time))
+        end
+        displacement = diffnorm(edgestate, (u1 + u2)/2)
+        edgestate = (u1 + u2)/2
+        counter += 1
+        
+        ProgressMeter.next!(progress;
+            showvalues = [(:Iteration, counter), (:"Edge point", edgestate)])
+        if verbose && (counter == maxiter)
+            @warn("Reached maximum number of $(maxiter) iterations.")
+        end
+    end
+
+    if verbose && (counter < maxiter)
+        println("Edge-tracking converged after $(counter) iterations.")
+    end
+    
+    track1 = StateSpaceSet(track1)
+    track2 = StateSpaceSet(track2)
+
+    return EdgeTrackingResults(
+        StateSpaceSet((vec(track1) .+ vec(track2))./2),
+        track1,
+        track2,
+        time,
+        bisect_idx,
+        true)
+end
+
+"""
+    bisect_to_edge(pds::ParallelDynamicalSystem, mapper::AttractorMapper; kwargs...) -> u1, u2
+Finds the basin boundary between two states `u1, u2 = current_states(pds)` by bisecting
+along a straight line in phase space. The states `u1` and `u2` must belong to different
+basins.
+
+Returns a triple `u1, u2, success`, where `u1, u2` are two new states located on either side
+of the basin boundary that lie less than `bisect_thresh` (Euclidean distance in state space)
+apart from each other, and `success` is a Bool indicating whether the bisection was
+successful (it may fail if the `mapper` maps both states to the same basin of attraction,
+in which case a warning is raised).
+
+## Keyword arguments
+* `bisect_thresh = 1e-7`: The maximum (Euclidean) distance between the two returned states.
+
+## Description
+`pds` is a `ParallelDynamicalSystem` with two states. The `mapper` must be an
+`AttractorMapper` of subtype `AttractorsViaProximity` or `AttractorsViaRecurrences`.
+
+!!! info
+    If the straight line between `u1` and `u2` intersects the basin boundary multiple
+    times, the method will find one of these intersection points. If more than two attractors
+    exist, one of the two returned states may belong to a different basin than the initial
+    conditions `u1` and `u2`. A warning is raised if the bisection involves a third basin.
+"""
+function bisect_to_edge(pds::ParallelDynamicalSystem, mapper::AttractorMapper;
+    bisect_thresh=1e-6,
+    verbose=true)
+
+    u1, u2 = current_states(pds)
+    idx1, idx2 = mapper(u1), mapper(u2)
+
+    if (idx1 == idx2)
+        if idx1 == -1
+            error("AttractorMapper returned label -1 (could not match the initial condition with any attractor).
+            Try changing the settings of the `AttractorMapper` or increasing bisect_thresh, diverge_thresh.")
+        else
+            if verbose
+            @warn "Both initial conditions belong to the same basin of attraction.
+                Attractor label: $(idx1)
+                u1 = $(u1)
+                u2 = $(u2)"
+            end
+            return u1, u2, false
+        end
+    end
+    
+    distance = diffnorm(u1, u2)
+    while distance > bisect_thresh
+        u_new = (u1 + u2)/2
+        idx_new = mapper(u_new)
+        # slightly shift u_new if it lands too close to the boundary
+        retry_counter = 1
+        while (idx_new == -1) && retry_counter < 3 # ToDO: make kwarg
+            if verbose
+                @warn "Shifting new point slightly because AttractorMapper returned -1"
+            end
+            u_new += bisect_thresh*(u1 - u2)
+            idx_new = mapper(u_new)
+            retry_counter += 1
+        end
+        # update u1 or u2
+        if idx_new == idx1
+            u1 = u_new
+        else
+            if idx_new != idx2
+                if idx_new == -1
+                    error("AttractorMapper returned label -1 (could not match the initial condition with any attractor.)
+                    Try changing the settings of AttractorsViaProximity or increasing bisect_thresh, diverge_thresh.")
+                else
+                    if verbose
+                        @warn "New bisection point belongs to a third basin of attraction."
+                    end
+                end
+            end
+            u2 = u_new
+        end
+        distance = diffnorm(u1, u2)
+    end
+    return u1, u2, true
+end
+
+function diffnorm(u1, u2)
+    d = zero(eltype(u1))
+    @inbounds for i in eachindex(u1)
+        d += (u1[i] - u2[i])^2
+    end
+    sqrt(d)
+end
+
+function diffnorm(pds::ParallelDynamicalSystem)
+    diffnorm(current_state(pds, 1), current_state(pds, 2))
+end
\ No newline at end of file
diff --git a/src/continuation/continuation_recurrences.jl b/src/continuation/continuation_recurrences.jl
index 10bb0326..83ea8475 100644
--- a/src/continuation/continuation_recurrences.jl
+++ b/src/continuation/continuation_recurrences.jl
@@ -56,7 +56,7 @@ on the outcome with different matching-related keywords.
 You do not need to compute attractors and basins again!
 
 Matching is a very sophisticated process that can be understood in detail by reading
-the docstrings of [`match_statespacesets!`](@ref) first and then [`match_continuation!`](@ref).
+the docstrings of [`match_statespacesets!!`](@ref) first and then [`match_continuation!`](@ref).
 Here is a short summary: attractors from previous and current parameter are matched
 based on their "distance". By default this is distance in state space, but any measure of
 "distance" may be used, such as the distance between Lyapunov spectra.
diff --git a/src/continuation/match_attractor_ids.jl b/src/continuation/match_attractor_ids.jl
index b2c68704..75ce4cda 100644
--- a/src/continuation/match_attractor_ids.jl
+++ b/src/continuation/match_attractor_ids.jl
@@ -1,11 +1,11 @@
 # Notice this file uses heavily `dict_utils.jl`!
-export match_statespacesets!, match_basins_ids!, replacement_map, match_continuation!
+export match_statespacesets!!, match_basins_ids!, replacement_map, match_continuation!
 
 ###########################################################################################
 # Matching attractors and key swapping business
 ###########################################################################################
 """
-    match_statespacesets!(a₊::AbstractDict, a₋; distance = Centroid(), threshold = Inf)
+    match_statespacesets!!(a₊::AbstractDict, a₋; distance = Centroid(), threshold = Inf)
 
 Given dictionaries `a₊, a₋` mapping IDs to `StateSpaceSet` instances,
 match the IDs in dictionary `a₊` so that its sets that are the closest to
@@ -39,7 +39,7 @@ Additionally, you can provide a `threshold` value. If the distance between two a
 is larger than this `threshold`, then it is guaranteed that the attractors will get assigned
 different key in the dictionary `a₊` (which is the next available integer).
 """
-function match_statespacesets!(a₊::AbstractDict, a₋; kwargs...)
+function match_statespacesets!!(a₊::AbstractDict, a₋; kwargs...)
     rmap = replacement_map(a₊, a₋; kwargs...)
     swap_dict_keys!(a₊, rmap)
     return rmap
@@ -52,7 +52,7 @@ end
     replacement_map(a₊, a₋; distance = Centroid(), threshold = Inf) → rmap
 
 Return a dictionary mapping keys in `a₊` to new keys in `a₋`,
-as explained in [`match_statespacesets!`](@ref).
+as explained in [`match_statespacesets!!`](@ref).
 """
 function replacement_map(a₊::Dict, a₋::Dict;
         distance = Centroid(), threshold = Inf, next_id = nothing
@@ -117,7 +117,7 @@ end
 """
     match_basins_ids!(b₊::AbstractArray, b₋; threshold = Inf)
 
-Similar to [`match_statespacesets!`](@ref) but operate on basin arrays instead
+Similar to [`match_statespacesets!!`](@ref) but operate on basin arrays instead
 (the arrays typically returned by [`basins_of_attraction`](@ref)).
 
 This method matches IDs of attractors whose basins of attraction before and after `b₋,b₊`
@@ -163,14 +163,14 @@ Loop over all entries in the given arguments (which are typically the direct out
 attractor IDs in both the attractors container and the basins fractions container.
 This means that we loop over each entry of the vectors (skipping the first),
 and in each entry we attempt to match the current dictionary keys to the keys of the
-previous dictionary using [`match_statespacesets`](@ref).
+previous dictionary using [`match_statespacesets!`](@ref).
 
-The keywords `distance, threshold` are propagated to [`match_statespacesets`](@ref).
+The keywords `distance, threshold` are propagated to [`match_statespacesets!`](@ref).
 However, there is a unique keyword for `match_continuation!`: `use_vanished::Bool`.
 If `true`, then attractors that existed before but have vanished are kept in "memory"
 when it comes to matching: the new attractors are compared to the latest instance
 of all attractors that have ever existed, and get matched to their closest ones
-as per [`match_statespacesets!`](@ref).
+as per [`match_statespacesets!!`](@ref).
 If `false`, vanished attractors are ignored. Note that in this case new attractors
 that cannot be matched to any previous attractors will get an appropriately
 incremented ID. E.g., if we started with three attractors, and attractor 3 vanished,
@@ -224,7 +224,7 @@ function _rematch_ignored!(fractions_curves, attractors_info; kwargs...)
         # and reappears, it will get a different (incremented) id as it should!
         next_id_a = max(maximum(keys(a₊)), maximum(keys(a₋))) + 1
         next_id = max(next_id+1, next_id_a)
-        rmap = match_statespacesets!(a₊, a₋; next_id, kwargs...)
+        rmap = match_statespacesets!!(a₊, a₋; next_id, kwargs...)
         swap_dict_keys!(fractions_curves[i+1], rmap)
     end
 end
@@ -240,7 +240,7 @@ function _rematch_with_past!(fractions_curves, attractors_info; kwargs...)
         for (k, A) in a₋
             latest_ghosts[k] = A
         end
-        rmap = match_statespacesets!(a₊, latest_ghosts; kwargs...)
+        rmap = match_statespacesets!!(a₊, latest_ghosts; kwargs...)
         swap_dict_keys!(fractions_curves[i+1], rmap)
     end
 end
diff --git a/src/mapping/attractor_mapping.jl b/src/mapping/attractor_mapping.jl
index 31c10b1f..24dfdb01 100644
--- a/src/mapping/attractor_mapping.jl
+++ b/src/mapping/attractor_mapping.jl
@@ -145,17 +145,18 @@ extract_attractors(::AttractorMapper) = error("not implemented")
 # It works for all mappers that define a `basins_fractions` method.
 """
     basins_of_attraction(mapper::AttractorMapper, grid::Tuple) → basins, attractors
+
 Compute the full basins of attraction as identified by the given `mapper`,
-which includes a reference to a [`GeneralizedDynamicalSystem`](@ref) and return them
+which includes a reference to a [`DynamicalSystem`](@ref) and return them
 along with (perhaps approximated) found attractors.
 
 `grid` is a tuple of ranges defining the grid of initial conditions that partition
 the state space into boxes with size the step size of each range.
 For example, `grid = (xg, yg)` where `xg = yg = range(-5, 5; length = 100)`.
 The grid has to be the same dimensionality as the state space expected by the
-integrator/system used in `mapper`. E.g., a [`projected_integrator`](@ref)
+integrator/system used in `mapper`. E.g., a [`ProjectedDynamicalSystem`](@ref)
 could be used for lower dimensional projections, etc. A special case here is
-a [`poincaremap`](@ref) with `plane` being `Tuple{Int, <: Real}`. In this special
+a [`PoincareMap`](@ref) with `plane` being `Tuple{Int, <: Real}`. In this special
 scenario the grid can be one dimension smaller than the state space, in which case
 the partitioning happens directly on the hyperplane the Poincaré map operates on.
 
diff --git a/src/mapping/grouping/attractor_mapping_featurizing.jl b/src/mapping/grouping/attractor_mapping_featurizing.jl
index 63221643..2c6ca729 100644
--- a/src/mapping/grouping/attractor_mapping_featurizing.jl
+++ b/src/mapping/grouping/attractor_mapping_featurizing.jl
@@ -27,7 +27,7 @@ end
 
 Initialize a `mapper` that maps initial conditions to attractors using a featurizing and
 grouping approach. This is a supercase of the featurizing and clustering approach that
-is utilized by bSTAB [Stender2021](@cite) and MCBB [Gelbrecht2021](@cite).
+is utilized by bSTAB [Stender2021](@cite) and MCBB [Gelbrecht2020](@cite).
 See [`AttractorMapper`](@ref) for how to use the `mapper`.
 This `mapper` also allows the syntax `mapper(u0)` but only if the `grouping_config`
 is _not_ `GroupViaClustering`.
diff --git a/src/mapping/recurrences/attractor_mapping_recurrences.jl b/src/mapping/recurrences/attractor_mapping_recurrences.jl
index b42ff9ea..b1ec8b7b 100644
--- a/src/mapping/recurrences/attractor_mapping_recurrences.jl
+++ b/src/mapping/recurrences/attractor_mapping_recurrences.jl
@@ -19,7 +19,7 @@ so that a finite state machine can operate on top of it. Possibilities are:
   `grid = (xg, yg)` with `xg = range(0, 10.0^(1/2); length = 200).^2,
   yg = range(-5, 5; length = 100)`.
 3. An instance of the special grid type
-  [`SubdividedBasedGrid`](@ref), which can be created either manually or by using
+  [`SubdivisionBasedGrid`](@ref), which can be created either manually or by using
   [`subdivision_based_grid`](@ref).
   This automatically analyzes and adapts grid discretization
   levels in accordance with state space flow speed in different regions.
@@ -178,7 +178,7 @@ extract_attractors(m::AttractorsViaRecurrences) = m.bsn_nfo.attractors
     basins_of_attraction(mapper::AttractorsViaRecurrences; show_progress = true)
 
 This is a special method of `basins_of_attraction` that using recurrences does
-_exactly_ what is described in the paper by Datseris & Wagemakers [Datseris2022](@ref).
+_exactly_ what is described in the paper by Datseris & Wagemakers [Datseris2022](@cite).
 By enforcing that the internal grid of `mapper` is the same as the grid of initial
 conditions to map to attractors, the method can further utilize found exit and attraction
 basins, making the computation faster as the grid is processed more and more.
diff --git a/src/plotting.jl b/src/plotting.jl
index 9f384808..7d9e1aa3 100644
--- a/src/plotting.jl
+++ b/src/plotting.jl
@@ -9,7 +9,7 @@ Plot a heatmap of found (2-dimensional) `basins` of attraction and corresponding
 
 ## Keyword arguments
 
-- All the [common plotting keywords](@ref).
+- All the [common plotting keywords](@ref common_plot_kwargs).
 """
 function heatmap_basins_attractors end
 function heatmap_basins_attractors! end
@@ -40,7 +40,7 @@ while `ds, prange, pidx` are the input to [`continuation`](@ref).
 - `markersize = 10`
 - `Δt, T`: propagated to `trajectory` for evolving an initial condition sampled
   from an attractor
-- Also all [common plotting keywords](@ref).
+- Also all [common plotting keywords](@ref common_plot_kwargs).
 """
 function animate_attractors_continuation end
 export animate_attractors_continuation
@@ -60,7 +60,7 @@ i.e., visualize the output of [`continuation`](@ref).
   if the style is `:band`
 - `axislegend_kwargs = (position = :lt,)`: propagated to `axislegend` if a legend is added
 - `series_kwargs = NamedTuple()`: propagated to the band or scatterline plot
-- Also all [common plotting keywords](@ref).
+- Also all [common plotting keywords](@ref common_plot_kwargs).
 """
 function plot_basins_curves end
 function plot_basins_curves! end
@@ -74,7 +74,7 @@ the parameter instead of their fraction.
 The function `attractor_to_real` takes as input a `StateSpaceSet` (attractor)
 and returns a real number so that it can be plotted versus the parameter axis.
 
-Same keywords as [`plot_basins_curves`](@ref).
+Same keywords as [`plot_basins_curves`](@ref common_plot_kwargs).
 """
 function plot_attractors_curves end
 function plot_attractors_curves! end
diff --git a/src/tipping/tipping_probabilities.jl b/src/tipping/tipping_probabilities.jl
index c0e915d7..e98b3eca 100644
--- a/src/tipping/tipping_probabilities.jl
+++ b/src/tipping/tipping_probabilities.jl
@@ -4,14 +4,14 @@ export tipping_probabilities
     tipping_probabilities(basins_before, basins_after) → P
 
 Return the tipping probabilities of the computed basins before and after a change
-in the system parameters (or time forcing), according to the definition of [Kaszás2019](@cite).
+in the system parameters (or time forcing), according to the definition of [Kaszas2019](@cite).
 
 The input `basins` are integer-valued arrays, where the integers enumerate the attractor,
 e.g. the output of [`basins_of_attraction`](@ref).
 
 ## Description
 Let ``\\mathcal{B}_i(p)`` denote the basin of attraction of attractor ``A_i`` at
-parameter(s) ``p``. Kaszás et al [Kaszás2019](@cite) define the tipping probability
+parameter(s) ``p``. Kaszás et al [Kaszas2019](@cite) define the tipping probability
 from ``A_i`` to ``A_j``, given a parameter change in the system of ``p_- \\to p_+``, as
 
 ```math
diff --git a/test/boundaries/edgetracking_tests.jl b/test/boundaries/edgetracking_tests.jl
new file mode 100644
index 00000000..72a5bee8
--- /dev/null
+++ b/test/boundaries/edgetracking_tests.jl
@@ -0,0 +1,77 @@
+using Attractors
+using Test
+using OrdinaryDiffEq
+using LinearAlgebra
+
+@testset "Saddle point of cubic map" begin
+    cubicmap(u, p, n) = SVector{1}(p[1]*u[1] - u[1]^3)
+    ds = DeterministicIteratedMap(cubicmap, [1.0], [2.0])
+    attrs = Dict(1 => StateSpaceSet([1.0]), 2 => StateSpaceSet([-1.0]))
+    saddle = edgetracking(ds, attrs; Δt=1, abstol=1e-8).edge[end]
+    @test saddle[1] < 1e-5
+end
+
+@testset "Saddle point of FitzHugh-Nagumo system" begin
+    fhn(u, p, t) = SVector{2}([10*(u[1] - u[1]^3 - u[2]), -3*u[2] + u[1]])
+    ds = CoupledODEs(fhn, ones(2), diffeq=(;alg = Vern9()))
+    fp = [sqrt(2/3), sqrt(2/27)]
+    attrs = Dict([1 => StateSpaceSet([fp]), 2 => StateSpaceSet([-fp])])
+    bisect_thresh, diverge_thresh, maxiter, abstol = 1e-8, 1e-7, 100, 1e-9
+    edge = edgetracking(ds, attrs; u1=[-1.0, 0.2], u2=[1.0, 0.2],
+        bisect_thresh, diverge_thresh, maxiter, abstol).edge
+    println(edge[end], edge[end-1])
+    @test sqrt(sum(abs, (edge[end]-zeros(2)).^2)) < 1e-6
+end
+
+@testset "Thomas' rule" begin
+    # Chaotic dynamical system
+    function thomas_rule(u, p, t)
+        x,y,z = u
+        b = p[1]
+        xdot = sin(y) - b*x
+        ydot = sin(z) - b*y
+        zdot = sin(x) - b*z
+        return SVector{3}(xdot, ydot, zdot)
+    end
+    ds = CoupledODEs(thomas_rule, [1.0, 0, 0], [0.16]; diffeq=(reltol=1e-12,))
+
+    # Find attractors on a 3D grid
+    xg = yg = yz = range(-6.0, 6.0; length = 101)
+    grid = (xg, yg, yz)
+    mapper = AttractorsViaRecurrences(ds, grid; consecutive_recurrences = 1000)
+    sampler, = statespace_sampler(grid)
+    basins_fractions(mapper, sampler)
+    attractors = extract_attractors(mapper)
+
+    # Run edgetracking between pairs of points lying on different attractors
+    n_sample = 25
+    pairs12, pairs13, pairs23 = [], [], []
+    for i in 1:Int(sqrt(n_sample))
+        _pairs12, _pairs13, _pairs23 = [], [], []
+        for j in 1:Int(sqrt(n_sample))
+            et12 = edgetracking(ds, attractors;
+                    u1=attractors[1][i], u2=attractors[2][j], bisect_thresh=1e-4,
+                    diverge_thresh=1e-3, maxiter=10000, abstol=1e-5, verbose=false)
+            et13 = edgetracking(ds, attractors;
+                    u1=attractors[1][i], u2=attractors[3][j], bisect_thresh=1e-4,
+                    diverge_thresh=1e-3, maxiter=10000, abstol=1e-5, verbose=false)
+            et23 = edgetracking(ds, attractors;
+                    u1=attractors[2][i], u2=attractors[3][j], bisect_thresh=1e-4,
+                    diverge_thresh=1e-3, maxiter=10000, abstol=1e-5, verbose=false)
+
+            et12.success ? push!(_pairs12, et12.edge[end]) : nothing
+            et13.success ? push!(_pairs13, et13.edge[end]) : nothing
+            et23.success ? push!(_pairs23, et23.edge[end]) : nothing
+        end
+        push!(pairs12, _pairs12)
+        push!(pairs13, _pairs13)
+        push!(pairs23, _pairs23)
+    end
+    edgestates = reduce(vcat, [pairs12 pairs13 pairs23])
+
+    # Verify that all found edge states have the same Euclidean norm of `norm_value`
+    # (due to the symmetry of the system `ds`)
+    norm_value = 4.06585
+    norm_deviations = [norm(edgestates[i]) - norm_value for i in 1:length(edgestates)]
+    @test maximum(norm_deviations) < 1e-3
+end
diff --git a/test/continuation/matching_attractors.jl b/test/continuation/matching_attractors.jl
index d6b1b4e5..18fdcb21 100644
--- a/test/continuation/matching_attractors.jl
+++ b/test/continuation/matching_attractors.jl
@@ -8,7 +8,7 @@ DO_EXTENSIVE_TESTS = get(ENV, "ATTRACTORS_EXTENSIVE_TESTS", "false") == "true"
     @testset "infinite threshold" begin
         a_afte = Dict(2 => [SVector(0.0, 0.0)], 1 => [SVector(2.0, 2.0)])
         a_afte = Dict(keys(a_afte) .=> StateSpaceSet.(values(a_afte)))
-        rmap = match_statespacesets!(a_afte, a_befo)
+        rmap = match_statespacesets!!(a_afte, a_befo)
         @test rmap == Dict(1 => 2, 2 => 1)
         @test a_afte[1] == a_befo[1] == StateSpaceSet([SVector(0.0, 0.0)])
         @test haskey(a_afte, 2)
@@ -17,7 +17,7 @@ DO_EXTENSIVE_TESTS = get(ENV, "ATTRACTORS_EXTENSIVE_TESTS", "false") == "true"
     @testset "separating threshold" begin
         a_afte = Dict(2 => [SVector(0.0, 0.0)], 1 => [SVector(2.0, 2.0)])
         a_afte = Dict(keys(a_afte) .=> StateSpaceSet.(values(a_afte)))
-        rmap = match_statespacesets!(a_afte, a_befo; threshold = 0.1)
+        rmap = match_statespacesets!!(a_afte, a_befo; threshold = 0.1)
         @test rmap == Dict(1 => 3, 2 => 1)
         @test a_afte[1] == a_befo[1] == StateSpaceSet([SVector(0.0, 0.0)])
         @test !haskey(a_afte, 2)
@@ -149,7 +149,7 @@ if DO_EXTENSIVE_TESTS
         d, α, ω = 0.3, 0.2, 0.5
         ds = magnetic_pendulum(; d, α, ω, diffeq)
         xg = yg = range(-3, 3, length = 100)
-        ds = projected_integrator(ds, 1:2, [0.0, 0.0])
+        ds = ProjectedDynamicalSystem(ds, 1:2, [0.0, 0.0])
         mapper = AttractorsViaRecurrences(ds, (xg, yg); sparse = false, Δt = 1.0)
         b₋, a₋ = basins_of_attraction(mapper; show_progress = false)
         # still 3 attractors at γ3 = 0.2, but only 2 at 0.1
@@ -158,7 +158,7 @@ if DO_EXTENSIVE_TESTS
             mapper = AttractorsViaRecurrences(ds, (xg, yg); sparse = false, Δt = 1.0)
             b₊, a₊ = basins_of_attraction(mapper; show_progress = false)
             @testset "distances match" begin
-                rmap = match_statespacesets!(a₊, a₋)
+                rmap = match_statespacesets!!(a₊, a₋)
                 for k in keys(a₊)
                     dist = minimum(norm(x .- y) for x ∈ a₊[k] for y ∈ a₋[k])
                     @test dist < 0.2
diff --git a/test/runtests.jl b/test/runtests.jl
index 5a7bec47..10312aff 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -34,4 +34,8 @@ testfile(file, testname=defaultname(file)) = @testset "$testname" begin; include
     @testset "mfs" begin
         testfile("mfs/mfstest.jl")
     end
+
+    @testset "boundaries" begin
+        testfile("boundaries/edgetracking_tests.jl")
+    end
 end