diff --git a/CHANGELOG.md b/CHANGELOG.md
index 06300e1c..0e8550b7 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,9 @@
+# v1.19
+
+- Global continuation can now be performed across any arbitrary curve
+  in parameter space. This is something completely novel
+  and we'll likely work on a paper on this!
+
 # v1.18
 
 This is a big release, with (hopefully) nothing breaking, but lots of deprecations!
diff --git a/Project.toml b/Project.toml
index 3884b66d..755a7981 100644
--- a/Project.toml
+++ b/Project.toml
@@ -2,8 +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.18.5"
-
+version = "1.19.4"
 
 [deps]
 BlackBoxOptim = "a134a8b2-14d6-55f6-9291-3336d3ab0209"
@@ -33,7 +32,7 @@ Clustering = "0.15"
 ComplexityMeasures = "2.3, 3"
 Distances = "0.7, 0.8, 0.9, 0.10"
 Distributions = "0.21, 0.22, 0.23, 0.24, 0.25"
-DynamicalSystemsBase = "3.0.3"
+DynamicalSystemsBase = "3.9.1"
 Makie = "≥ 0.19"
 Neighborhood = "0.2.2"
 Optim = "1.7"
diff --git a/docs/src/bfkit_comparison.jl b/docs/src/bfkit_comparison.jl
index 43f36f45..315b2622 100644
--- a/docs/src/bfkit_comparison.jl
+++ b/docs/src/bfkit_comparison.jl
@@ -210,7 +210,7 @@ plot_attractors_curves(
 
 # ## Discussion and comparison
 
-# On the one hand, Attractors.jl found not only a single limit cycle,
+# Attractors.jl found not only a single limit cycle,
 # but also all system attractors, including chaotic ones.
 # It didn't require any initial guess regarding the limit cycle or its period,
 # but only a state space box that may contain attractors.
@@ -224,9 +224,17 @@ plot_attractors_curves(
 # And finally, Attractors.jl estimates a more general nonlocal measure of stability,
 # in the sense that if a set is nonlocally stable, it is guaranteed to be locally stable,
 # however the other way around isn't guaranteed.
-
-# On the other hand, traditional local continuation can track the unstable
-# branches, and automatically detect and label local bifurcations
+# Moreover, due to the orthogonality of finding and matching attractors,
+# as well as finding _all_ attractors, the global continuation of Attractors.jl
+# can continue along arbitrary user-defined curves in parameter space; not just
+# along a single parameter axis. This is possible because it is completely fine
+# for some attractors to stop existing during the global continuation,
+# while local continuation stops when attractors (and their unstable version)
+# stop existing.
+
+# Traditional local continuation can track the unstable
+# branches, and automatically detect and label local bifurcations,
+# both of which are not possible in Attractors.jl
 # (note we didn't bother to plot any of the detected bifurcations - if any are found).
 # In our experience having the local bifurcations is always useful.
 # Now, whether the unstable branches are of a limit cycle are useful or not,
@@ -248,6 +256,7 @@ plot_attractors_curves(
 # Using a slightly incorrect initial period guess of 20.0 instead of 19.0 also fails.
 # We imagine that this sensitivity would apply also to some other of the
 # several meta-parameters that enter a traditional continuation routine,
+# for example the thresholds and accuracy parameters related to Newton convergence,
 # but we didn't check further. This is exactly what we were alluding to
 # in the comparison we did in [Datseris2023](@cite), that traditional
 # local continuation "requires expertise and constant interventions".
diff --git a/docs/src/examples.md b/docs/src/examples.md
index 9545a7f2..c45c1e75 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -2,7 +2,7 @@
 
 Note that the examples utilize some convenience plotting functions offered by Attractors.jl which come into scope when using `Makie` (or any of its backends such as `CairoMakie`), see the [visualization utilities](@ref) for more.
 
-## Newton's fractal (basins of 2D map)
+## Newton's fractal (basins of a 2D map)
 
 ```@example MAIN
 using Attractors
diff --git a/docs/src/index.md b/docs/src/index.md
index c84492be..d62e7f0a 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -10,12 +10,8 @@ using CairoMakie, Attractors
 
 ## Latest news
 
-- Brand new [tutorial](@ref tutorial) for the whole repo!
-- `continuation` has been renamed to `global_continuation`
-- Brand new continuation algorithm `AttractorSeedContinueMatch`
-  that generalizes existing infrastructure and makes Attractors.jl
-  more composable and more extendable.
-- New "matchers" and extendable interface for them.
+- Global continuation can now be performed across any arbitrary curve
+  in parameter space.
 - See the CHANGELOG.md (at the GitHub repo) for more!
 
 ## Getting started
diff --git a/docs/src/tutorial.jl b/docs/src/tutorial.jl
index b911177a..3de773db 100644
--- a/docs/src/tutorial.jl
+++ b/docs/src/tutorial.jl
@@ -1,5 +1,13 @@
 # # [Attractors.jl Tutorial](@id tutorial)
 
+# ```@raw html
+# <video width="auto" controls loop>
+# <source src="../attracont.mp4" type="video/mp4">
+# </video>
+# ```
+
+#nb @doc Attractors
+
 # [`Attractors`](@ref) is a component of the **DynamicalSystems.jl** library.
 # This tutorial will walk you through its main functionality.
 # That is, given a `DynamicalSystem` instance, find all its attractors and their basins
@@ -13,6 +21,18 @@
 # but we won't cover anything else in this introductory tutorial.
 # See the [examples](@ref examples) page instead.
 
+
+# ### Package versions used
+
+import Pkg
+
+#nb # Activate an environment in the folder containing the notebook
+#nb Pkg.activate(dirname(@__DIR__))
+#nb Pkg.add(["DynamicalSystems", "CairoMakie", "GLMakie", "OrdinaryDiffEq", "BenchmarkTools"])
+
+Pkg.status(["Attractors", "CairoMakie", "OrdinaryDiffEq"])
+
+
 # ## Tutorial - copy-pasteable version
 
 # _Gotta go fast!_
@@ -294,7 +314,7 @@ pidx = 1 # index of the parameter
 # Then, we may call the [`global_continuation`](@ref) function.
 # We have to provide a continuation algorithm, which itself references an [`AttractorMapper`](@ref).
 # In this example we will re-use the `mapper` to create the "flagship product" of Attractors.jl
-# which is the generic [`AttractorSeedContinueMatch`](@ref).
+# which is the geenral [`AttractorSeedContinueMatch`](@ref).
 # This algorithm uses the `mapper` to find all attractors at each parameter value
 # and from the found attractors it continues them along a parameter axis
 # using a seeding process (see its documentation string).
@@ -486,6 +506,45 @@ for ax in (ax1, ax2,); hidexdecorations!(ax; grid = false); end
 resize!(fig, 500, 500)
 fig
 
+# ## Continuation along arbitrary parameter curves
+
+# One of the many advantages of the global continuation is that we can choose
+# what parameters to continue over. We can provide any arbitrary curve
+# in parameter space. This is possible because (1) finding and matching attractors
+# are two completely orthogonal steps, and (2) it is completely fine for
+# attractors to dissapear (and perhaps re-appear) during a global continuation.
+
+#For example, we can probe an elipsoid defined as
+
+params(θ) = [1 => 5 + 0.5cos(θ), 2 => 0.1 + 0.01sin(θ)]
+pcurve = params.(range(0, 2π; length = 101))
+
+# here each component maps the parameter index to its value.
+# We can just give this `pcurve` to the global continuation,
+# using the same mapper and continuation algorithm,
+# but adjusting the matching process so that vanished attractors
+# are kept in "memory"
+
+matcher = MatchBySSSetDistance(use_vanished = true)
+
+ascm = AttractorSeedContinueMatch(mapper, matcher)
+
+fractions_cont, attractors_cont = global_continuation(
+	ascm, pcurve, sampler; samples_per_parameter = 1_000
+)
+
+# and animate the result
+animate_attractors_continuation(
+    ds, attractors_cont, fractions_cont, pcurve;
+    savename = "curvecont.mp4"
+);
+
+# ```@raw html
+# <video width="auto" controls loop>
+# <source src="../curvecont.mp4" type="video/mp4">
+# </video>
+# ```
+
 # ## Conclusion and comparison with traditional local continuation
 
 # We've reached the end of the tutorial! Some aspects we haven't highlighted is
diff --git a/ext/plotting.jl b/ext/plotting.jl
index 3d66f76b..b2c68812 100644
--- a/ext/plotting.jl
+++ b/ext/plotting.jl
@@ -395,7 +395,13 @@ end
 # Videos
 ##########################################################################################
 function Attractors.animate_attractors_continuation(
-        ds::DynamicalSystem, attractors_cont, fractions_cont, prange, pidx;
+    ds::DynamicalSystem, attractors_cont, fractions_cont, prange, pidx; kw...)
+    pcurve = [[pidx => p] for p in prange]
+    return animate_attractors_continuation(ds, attractors_cont, fractions_cont, pcurve; kw...)
+end
+
+function Attractors.animate_attractors_continuation(
+        ds::DynamicalSystem, attractors_cont, fractions_cont, pcurve;
         savename = "attracont.mp4", access = SVector(1, 2),
         limits = auto_attractor_lims(attractors_cont, access),
         framerate = 4, markersize = 10,
@@ -406,10 +412,11 @@ function Attractors.animate_attractors_continuation(
         T = 100,
         figure = NamedTuple(), axis = NamedTuple(), fracaxis = NamedTuple(),
         legend = NamedTuple(),
+        add_legend = length(ukeys) ≤ 6
     )
     length(access) ≠ 2 && error("Need two indices to select two dimensions of `ds`.")
     K = length(ukeys)
-    fig = Figure(figure...)
+    fig = Figure(; figure...)
     ax = Axis(fig[1,1]; limits, axis...)
     fracax = Axis(fig[1,2]; width = 50, limits = (0,1,0,1), ylabel = "fractions",
         yaxisposition = :right, fracaxis...
@@ -423,19 +430,21 @@ function Attractors.animate_attractors_continuation(
     for k in ukeys
         plotf!(ax, att_obs[k]; color = (colors[k], 0.75), label = "$k", markersize, marker = markers[k])
     end
-    axislegend(ax; legend...)
+    if add_legend
+        axislegend(ax; legend...)
+    end
 
     # setup fractions axis
     heights = Observable(fill(0.1, K))
     barcolors = [colors[k] for k in ukeys]
     barplot!(fracax, fill(0.5, K), heights; width = 1, gap = 0, stack=1:K, color = barcolors)
 
-    record(fig, savename, eachindex(prange); framerate) do i
-        p = prange[i]
-        ax.title = "p = $p"
+    record(fig, savename, eachindex(pcurve); framerate) do i
+        p = pcurve[i]
+        ax.title = "p: $p" # TODO: Add compat printing here.
         attractors = attractors_cont[i]
         fractions = fractions_cont[i]
-        set_parameter!(ds, pidx, p)
+        set_parameters!(ds, p)
         heights[] = [get(fractions, k, 0) for k in ukeys]
 
         for (k, att) in attractors
diff --git a/src/basins/basins_utilities.jl b/src/basins/basins_utilities.jl
index 72889843..33a7b093 100644
--- a/src/basins/basins_utilities.jl
+++ b/src/basins/basins_utilities.jl
@@ -1,3 +1,46 @@
+# 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 [`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 [`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
+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.
+
+`basins_of_attraction` function is a convenience 5-lines-of-code wrapper which uses the
+`labels` returned by [`basins_fractions`](@ref) and simply assigns them to a full array
+corresponding to the state space partitioning indicated by `grid`.
+
+See also [`convergence_and_basins_of_attraction`](@ref).
+"""
+function basins_of_attraction(mapper::AttractorMapper, grid::Tuple; kwargs...)
+    basins = zeros(Int32, map(length, grid))
+    I = CartesianIndices(basins)
+    A = StateSpaceSet([generate_ic_on_grid(grid, i) for i in vec(I)])
+    fs, labels = basins_fractions(mapper, A; kwargs...)
+    attractors = extract_attractors(mapper)
+    vec(basins) .= vec(labels)
+    return basins, attractors
+end
+
+# Type-stable generation of an initial condition given a grid array index
+@generated function generate_ic_on_grid(grid::NTuple{B, T}, ind) where {B, T}
+    gens = [:(grid[$k][ind[$k]]) for k=1:B]
+    quote
+        Base.@_inline_meta
+        @inbounds return SVector{$B, Float64}($(gens...))
+    end
+end
+
 """
     basins_fractions(basins::AbstractArray [,ids]) → fs::Dict
 
diff --git a/src/basins/fractality_of_basins.jl b/src/basins/fractality_of_basins.jl
index be465ed1..58e0b3af 100644
--- a/src/basins/fractality_of_basins.jl
+++ b/src/basins/fractality_of_basins.jl
@@ -3,20 +3,32 @@ export uncertainty_exponent, basins_fractal_dimension, basins_fractal_test, basi
 """
     basin_entropy(basins::Array, ε = 20) -> Sb, Sbb
 
-Compute the basin entropy [Daza2016](@cite) `Sb` and basin boundary entropy `Sbb`
-of the given `basins` of attraction by considering `ε` boxes along each dimension.
+Return the basin entropy [Daza2016](@cite) `Sb` and basin boundary entropy `Sbb`
+of the given `basins` of attraction by considering `ε`-sized boxes along each dimension.
 
 ## Description
 
-First, the input `basins`
-is divided regularly into n-dimensional boxes of side `ε` (along all dimensions).
-Then `Sb` is simply the average of the Gibbs entropy computed over these boxes. The
-function returns the basin entropy `Sb` as well as the boundary basin entropy `Sbb`.
-The later is the average of the entropy only for boxes that contains at least two
+First, the n-dimensional input `basins`
+is divided regularly into n-dimensional boxes of side `ε`.
+If `ε` is an integer, the same size is used for all dimensions, otherwise `ε` can be
+a tuple with the same size as the dimensions of `basins`.
+Assuming that there are ``N`` `ε`-boxes that cover the `basins`, the basin entropy is estimated
+as [Daza2016](@cite)
+
+```math
+S_b = \\tfrac{1}{N}\\sum_{i=1}^{N}\\sum_{j=1}^{m_i}-p_{ij}\\log(p_{ij})
+```
+where ``m`` is the number of unique IDs (integers of `basins`) in box ``i``
+and ``p_{ij}`` is the relative frequency (probability) to obtain ID ``j``
+in the ``i`` box (simply the count of IDs ``j`` divided by the total in the box).
+
+`Sbb` is the boundary basin entropy.
+This follows the same definition as ``S_b``, but now averaged over only
+only boxes that contains at least two
 different basins, that is, for the boxes on the boundaries.
 
 The basin entropy is a measure of the uncertainty on the initial conditions of the basins.
-It is maximum at the value `log(n_att)` being `n_att` the number of attractors. In
+It is maximum at the value `log(n_att)` being `n_att` the number of unique IDs in `basins`. In
 this case the boundary is intermingled: for a given initial condition we can find
 another initial condition that lead to another basin arbitrarily close. It provides also
 a simple criterion for fractality: if the boundary basin entropy `Sbb` is above `log(2)`
@@ -25,35 +37,40 @@ have a fractal boundary, for a more precise test see [`basins_fractal_test`](@re
 An important feature of the basin entropy is that it allows
 comparisons between different basins using the same box size `ε`.
 """
-function basin_entropy(basins, ε = 20)
-    dims = size(basins)
-    vals = unique(basins)
-    Sb = 0; Nb = 0; N = 0
-    bx_tuple = ntuple(i -> range(1, dims[i] - rem(dims[i],ε), step = ε), length(dims))
-    box_indices = CartesianIndices(bx_tuple)
-    for box in box_indices
-        # compute the range of indices for the current box
-        I = CartesianIndices(ntuple(i -> range(box[i], box[i]+ε-1, step = 1), length(dims)))
-        box_values = [basins[k] for k in I]
-        N = N + 1
-        Nb = Nb + (length(unique(box_values)) > 1)
-        Sb = Sb + _box_entropy(box_values)
+function basin_entropy(basins::AbstractArray{Int, D}, ε::Int = 20) where {D}
+    es = ntuple(i -> ε, Val(D))
+    return basin_entropy(basins, es)
+end
+
+function basin_entropy(basins::AbstractArray{Int, D}, es::Dims{D}) where {D}
+    Sb = 0.0; Nb = 0
+    εranges = map((d, ε) -> 1:ε:d, size(basins), es)
+    box_iterator = Iterators.product(εranges...)
+    for box_start in box_iterator
+        box_ranges = map((d, ε) -> d:(d+ε-1), box_start, es)
+        box_values = view(basins, box_ranges...)
+        uvals = unique(box_values)
+        if length(uvals) > 1
+            Nb += 1
+            # we only need to estimate entropy for boxes with more than 1 val,
+            # because in other cases the entropy is zero
+            Sb = Sb + _box_entropy(box_values, uvals)
+        end
     end
-    return Sb/N, Sb/Nb
+    return Sb/length(box_iterator), Sb/Nb
 end
 
-function _box_entropy(box_values)
-    h = 0.
-    for (k,v) in enumerate(unique(box_values))
-        p = count( x -> (x == v), box_values)/length(box_values)
-        h += p*log(1/p)
+function _box_entropy(box_values, unique_vals = unique(box_values))
+    h = 0.0
+    for v in unique_vals
+        p = count(x -> (x == v), box_values)/length(box_values)
+        h += -p*log(p)
     end
     return h
 end
 
 
 
-
 """
     basins_fractal_test(basins; ε = 20, Ntotal = 1000) -> test_res, Sbb
 
@@ -87,15 +104,12 @@ the estimated value of the boundary basin entropy with the sampling method.
 """
 function basins_fractal_test(basins; ε = 20, Ntotal = 1000)
     dims = size(basins)
-    vals = unique(basins)
-    S=Int(length(vals))
-
     # Sanity check.
     if minimum(dims)/ε < 50
         @warn "Maybe the size of the grid is not fine enough."
     end
     if Ntotal < 100
-        error("Ntotal must be larger than 1000 to gather enough statistics.")
+        error("Ntotal must be larger than 100 to gather enough statistics.")
     end
 
     v_pts = zeros(Float64, length(dims), prod(dims))
@@ -105,7 +119,7 @@ function basins_fractal_test(basins; ε = 20, Ntotal = 1000)
     end
     tree = searchstructure(KDTree, v_pts, Euclidean())
     # Now get the values in the boxes.
-    Nb = 1; N = 1; Sb = 0;
+    Nb = 1
     N_stat = zeros(Ntotal)
     while Nb < Ntotal
         p = [rand()*(sz-ε)+ε for sz in dims]
diff --git a/src/continuation/basins_fractions_continuation_api.jl b/src/continuation/basins_fractions_continuation_api.jl
index 1b4fbe4e..5194af8d 100644
--- a/src/continuation/basins_fractions_continuation_api.jl
+++ b/src/continuation/basins_fractions_continuation_api.jl
@@ -14,6 +14,7 @@ abstract type GlobalContinuationAlgorithm end
 
 """
     global_continuation(gca::GlobalContinuationAlgorithm, prange, pidx, ics; kwargs...)
+    global_continuation(gca::GlobalContinuationAlgorithm, pcurve, ics; kwargs...)
 
 Find and continue attractors (or representations of attractors)
 and the fractions of their basins of attraction across a parameter range.
@@ -27,7 +28,15 @@ are given when creating `gca`.
 
 The basin fractions and the attractors (or some representation of them) are continued
 across the parameter range `prange`, for the parameter of the system with index `pidx`
-(any index valid in `DynamicalSystems.set_parameter!` can be used).
+(any index valid in `DynamicalSystems.set_parameter!` can be used). In contrast to
+traditional continuation (see online Tutorial for a comparison), global continuation
+can be performed over arbitrary user-defined curves in parameter space.
+The second call signature with `pcurve` allows for this possibility. In this case
+`pcurve` is a vector of iterables, where each itereable maps parameter indices
+to parameter values. These iterables can be dictionaries, named tuples, `Vector{Pair}`,
+etc., and the sequence of the iterables defines a curve in parameter space.
+In fact, the version with `prange, pidx` simply defines
+`pcurve = [[pidx => p] for p in prange]` and calls the second method.
 
 `ics` are the initial conditions to use when globally sampling the state space.
 Like in [`basins_fractions`](@ref) it can be either a set vector of initial conditions,
@@ -35,27 +44,29 @@ or a 0-argument function that generates random initial conditions.
 
 Possible subtypes of `GlobalContinuationAlgorithm` are:
 
-- [`RecurrencesFindAndMatch`](@ref)
+- [`AttractorSeedContinueMatch`](@ref)
 - [`FeaturizeGroupAcrossParameter`](@ref)
 
 ## Return
 
 1. `fractions_cont::Vector{Dict{Int, Float64}}`. The fractions of basins of attraction.
    `fractions_cont[i]` is a dictionary mapping attractor IDs to their basin fraction
-   at the `i`-th parameter.
-2. `attractors_cont::Vector{Dict{Int, <:Any}}`. Information about the attractors.
-   `attractors_cont[i]` is a dictionary mapping attractor ID to information about the
-   attractor at the `i`-th parameter.
-   The type of information stored depends on the chosen global continuation type,
-   but typically it is the attractors themselves as `StateSpaceSet`s.
+   at the `i`-th parameter combination.
+2. `attractors_cont::Vector{Dict{Int, <:Any}}`. The continued attractors.
+   `attractors_cont[i]` is a dictionary mapping attractor ID to the
+   attractor set at the `i`-th parameter combination.
 
 ## Keyword arguments
 
 - `show_progress = true`: display a progress bar of the computation.
 - `samples_per_parameter = 100`: amount of initial conditions sampled at each parameter
-  from `ics` if `ics` is a function instead of set initial conditions.
+  combination from `ics` if `ics` is a function instead of set initial conditions.
 """
-function global_continuation end
+function global_continuation(alg::GlobalContinuationAlgorithm, prange::AbstractVector, pidx, sampler; kw...)
+    # everything is propagated to the curve setting
+    pcurve = [[pidx => p] for p in prange]
+    return global_continuation(alg, pcurve, sampler; kw...)
+end
 
 include("continuation_ascm_generic.jl")
 include("continuation_recurrences.jl")
diff --git a/src/continuation/continuation_ascm_generic.jl b/src/continuation/continuation_ascm_generic.jl
index b4058084..b95ddf31 100644
--- a/src/continuation/continuation_ascm_generic.jl
+++ b/src/continuation/continuation_ascm_generic.jl
@@ -8,8 +8,8 @@ using Random: MersenneTwister
 A global continuation method for [`global_continuation`](@ref).
 `mapper` is any subtype of [`AttractorMapper`](@ref) which implements
 [`extract_attractors`](@ref), i.e., it finds the actual attractors.
-`matcher` is a configuration of how to match attractor IDs,
-and at the moment can only be an instance of [`MatchBySSSetDistance`](@ref).
+`matcher` is a configuration of how to match attractor IDs, see [`IDMatcher`](@ref)
+for more options.
 
 ## Description
 
@@ -116,17 +116,17 @@ function _default_seeding(attractor::AbstractStateSpaceSet)
     return (attractor[1],) # must be iterable
 end
 
-function global_continuation(acam::AttractorSeedContinueMatch, prange, pidx, ics;
+function global_continuation(acam::AttractorSeedContinueMatch, pcurve, ics;
         samples_per_parameter = 100, show_progress = true,
     )
     N = samples_per_parameter
-    progress = ProgressMeter.Progress(length(prange);
+    progress = ProgressMeter.Progress(length(pcurve);
         desc = "Continuing attractors and basins:", enabled=show_progress
     )
     mapper = acam.mapper
     reset_mapper!(mapper)
     # first parameter is run in isolation, as it has no prior to seed from
-    set_parameter!(referenced_dynamical_system(mapper), pidx, prange[1])
+    set_parameters!(referenced_dynamical_system(mapper), pcurve[1])
     if ics isa Function
         fs = basins_fractions(mapper, ics; show_progress = false, N = samples_per_parameter)
     else # we ignore labels in this continuation algorithm
@@ -137,10 +137,10 @@ function global_continuation(acam::AttractorSeedContinueMatch, prange, pidx, ics
     # The attractors are also stored (and are the primary output)
     prev_attractors = deepcopy(extract_attractors(mapper))
     attractors_cont = [deepcopy(prev_attractors)] # we need the copy
-    ProgressMeter.next!(progress; showvalues = [("previous parameter", prange[1]),])
+    ProgressMeter.next!(progress)
     # Continue loop over all remaining parameters
-    for (j, p) in enumerate(prange[2:end])
-        set_parameter!(referenced_dynamical_system(mapper), pidx, p)
+    for p in @view(pcurve[2:end])
+        set_parameters!(referenced_dynamical_system(mapper), p)
         reset_mapper!(mapper)
         # Seed initial conditions from previous attractors
         # Notice that one of the things that happens here is some attractors have
@@ -161,10 +161,10 @@ function global_continuation(acam::AttractorSeedContinueMatch, prange, pidx, ics
         push!(attractors_cont, current_attractors)
         # This is safe due to the deepcopies
         overwrite_dict!(prev_attractors, current_attractors)
-        ProgressMeter.next!(progress; showvalues = [("previous parameter", p),])
+        ProgressMeter.next!(progress)
     end
     rmaps = match_sequentially!(
-        attractors_cont, acam.matcher; prange, ds = referenced_dynamical_system(mapper), pidx
+        attractors_cont, acam.matcher; pcurve, ds = referenced_dynamical_system(mapper)
     )
     match_sequentially!(fractions_cont, rmaps)
     return fractions_cont, attractors_cont
diff --git a/src/continuation/continuation_grouping.jl b/src/continuation/continuation_grouping.jl
index 037c7a04..6e717de0 100644
--- a/src/continuation/continuation_grouping.jl
+++ b/src/continuation/continuation_grouping.jl
@@ -54,20 +54,21 @@ function mean_across_features(fs)
 end
 
 function global_continuation(
-        continuation::FeaturizeGroupAcrossParameter, prange, pidx, ics;
+        continuation::FeaturizeGroupAcrossParameter, pcurve, ics;
         show_progress = true, samples_per_parameter = 100
     )
     (; mapper, info_extraction, par_weight) = continuation
-    spp, n = samples_per_parameter, length(prange)
-    features = _get_features_prange(mapper, ics, n, spp, prange, pidx, show_progress)
+    spp, n = samples_per_parameter, length(pcurve)
+    features = _get_features_pcurve(mapper, ics, n, spp, pcurve, show_progress)
 
     # This is a special clause for implementing the MCBB algorithm (weighting
     # also by parameter value, i.e., making the parameter value a feature)
     # It calls a special `group_features` function that also incorporates the
     # parameter value (see below). Otherwise, we call normal `group_features`.
+    # TODO: We have deprecated this special clause. In the next version we need to cleanup
+    # the source code and remove the `par_weight` and its special treatment in `group_features`.
     if mapper.group_config isa GroupViaClustering && par_weight ≠ 0
         labels = group_features(features, mapper.group_config; par_weight, plength = n, spp)
-
     else
         labels = group_features(features, mapper.group_config)
     end
@@ -76,7 +77,7 @@ function global_continuation(
     return fractions_cont, attractors_cont
 end
 
-function _get_features_prange(mapper::AttractorsViaFeaturizing, ics, n, spp, prange, pidx, show_progress)
+function _get_features_pcurve(mapper::AttractorsViaFeaturizing, ics, n, spp, pcurve, show_progress)
     progress = ProgressMeter.Progress(n;
         desc="Generating features", enabled=show_progress, offset = 2,
     )
@@ -84,8 +85,8 @@ function _get_features_prange(mapper::AttractorsViaFeaturizing, ics, n, spp, pra
     feature = extract_features(mapper, ics; N = 1)
     features = Vector{typeof(feature[1])}(undef, n*spp)
     # Collect features
-    for (i, p) in enumerate(prange)
-        set_parameter!(mapper.ds, pidx, p)
+    for (i, p) in enumerate(pcurve)
+        set_parameters!(mapper.ds, p)
         current_features = extract_features(mapper, ics; show_progress, N = spp)
         features[((i - 1)*spp + 1):i*spp] .= current_features
         ProgressMeter.next!(progress)
diff --git a/src/dict_utils.jl b/src/dict_utils.jl
index c5d78425..fb765fd1 100644
--- a/src/dict_utils.jl
+++ b/src/dict_utils.jl
@@ -65,6 +65,7 @@ end
 
 """
     retract_keys_to_consecutive(v::Vector{<:Dict}) → rmap
+
 Given a vector of dictionaries with various positive integer keys, retract all keys so that
 consecutive integers are used. So if the dictionaries have overall keys 2, 3, 42,
 then they will transformed to 1, 2, 3.
diff --git a/src/mapping/attractor_mapping.jl b/src/mapping/attractor_mapping.jl
index bf04ae5e..77e6d93c 100644
--- a/src/mapping/attractor_mapping.jl
+++ b/src/mapping/attractor_mapping.jl
@@ -174,52 +174,6 @@ be used instead of [`basins_fractions`](@ref).
 """
 function convergence_time end
 
-#########################################################################################
-# Generic basins of attraction method structure definition
-#########################################################################################
-# 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 [`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 [`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
-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.
-
-`basins_of_attraction` function is a convenience 5-lines-of-code wrapper which uses the
-`labels` returned by [`basins_fractions`](@ref) and simply assigns them to a full array
-corresponding to the state space partitioning indicated by `grid`.
-
-See also [`convergence_and_basins_of_attraction`](@ref).
-"""
-function basins_of_attraction(mapper::AttractorMapper, grid::Tuple; kwargs...)
-    basins = zeros(Int32, map(length, grid))
-    I = CartesianIndices(basins)
-    A = StateSpaceSet([generate_ic_on_grid(grid, i) for i in vec(I)])
-    fs, labels = basins_fractions(mapper, A; kwargs...)
-    attractors = extract_attractors(mapper)
-    vec(basins) .= vec(labels)
-    return basins, attractors
-end
-
-# Type-stable generation of an initial condition given a grid array index
-@generated function generate_ic_on_grid(grid::NTuple{B, T}, ind) where {B, T}
-    gens = [:(grid[$k][ind[$k]]) for k=1:B]
-    quote
-        Base.@_inline_meta
-        @inbounds return SVector{$B, Float64}($(gens...))
-    end
-end
-
 #########################################################################################
 # Includes
 #########################################################################################
diff --git a/src/mapping/grouping/attractor_mapping_featurizing.jl b/src/mapping/grouping/attractor_mapping_featurizing.jl
index c8d7b87f..561eb215 100644
--- a/src/mapping/grouping/attractor_mapping_featurizing.jl
+++ b/src/mapping/grouping/attractor_mapping_featurizing.jl
@@ -166,6 +166,10 @@ function extract_features_single(mapper, ics; show_progress = true, N = 1000)
     return feature_vector
 end
 
+function (mapper::AttractorsViaFeaturizing)(u0)
+   f = extract_features_single(mapper, [u0]) 
+   return feature_to_group(f[1], mapper.group_config) 
+end
 
 # TODO: We need an alternative to deep copying integrators that efficiently
 # initializes integrators for any given kind of system. But that can be done
diff --git a/src/matching/basin_enclosure.jl b/src/matching/basin_enclosure.jl
index 742287e1..52ef418b 100644
--- a/src/matching/basin_enclosure.jl
+++ b/src/matching/basin_enclosure.jl
@@ -44,7 +44,8 @@ which is estimated with the `distance` keyword, which can be anything
 [`MatchBySSSetDistance`](@ref) accepts. The closest `₊` attractor gets the
 ID of the `₋` closest attractor that converge to it.
 
-Basin enclosure is a concept similar to "basin instability" in [Ritchie2023](@cite).
+Basin enclosure is a concept similar to "basin (in)stability" in [Ritchie2023](@cite):
+attractors that quantify as "basin stable" are matched.
 """
 @kwdef struct MatchByBasinEnclosure{E, D, S, T} <: IDMatcher
     ε::E = nothing
@@ -56,14 +57,14 @@ end
 
 function matching_map(
         current_attractors, prev_attractors, matcher::MatchByBasinEnclosure;
-        ds, pidx, p, pprev, next_id = next_free_id(current_attractors, prev_attractors)
+        ds, p, pprev = nothing, next_id = next_free_id(current_attractors, prev_attractors)
     )
     if matcher.ε === nothing
         e = ε_from_centroids(current_attractors)
     else
         e = matcher.ε
     end
-    set_parameter!(ds, pidx, p)
+    set_parameters!(ds, p)
     proximity = AttractorsViaProximity(ds, current_attractors, e;
         horizon_limit = Inf, Ttr = 0, consecutive_lost_steps = matcher.consecutive_lost_steps
     )
@@ -71,6 +72,7 @@ function matching_map(
     # to where they flowed to in current attractors
     # (notice that `proximity(u)` returns IDs of current attractors)
     flow = Dict(k => proximity(matcher.seeding(A)) for (k, A) in prev_attractors)
+
     # of course, the matching map is the inverse of `flow`
     rmap = Dict{Int, Int}()
     # but we need to take care of diverging and co-flowing attractors.
@@ -84,18 +86,20 @@ function matching_map(
         end
     end
     # next up are the co-flowing attractors
-    grouped_flows = _grouped_flows(flow)
+    allnewids = keys(current_attractors) #needed because flow only includes new ids that prev attractors flowed onto
+    grouped_flows = _grouped_flows(flow, allnewids) #map all current ids to prev ids that flowed onto them (if the att is new, no prev ids flowed, and the corresponding entry is empty)
+
     # notice the keys of `grouped_flows` are new IDs, same as with `rmap`.
     for (new_ID, old_flowed_to_same) in grouped_flows
-        if length(old_flowed_to_same) == 0
-            continue # none of the old IDs converged to the current `new_ID`
+        if length(old_flowed_to_same) == 0 # none of the old IDs converged to the current `new_ID`
+            rmap[new_ID] = next_id #necessary to make rmap complete and avoid skipping keys
+            next_id += 1
         elseif length(old_flowed_to_same) == 1
             rmap[new_ID] = only(old_flowed_to_same)
         else # need to resolve coflowing using distances
             a₊ = Dict(new_ID => current_attractors[new_ID])
             a₋ = Dict(old_ID => prev_attractors[old_ID] for old_ID in old_flowed_to_same)
             ssmatcher = MatchBySSSetDistance(; distance = matcher.distance)
-            @show length(a₊), length(a₋)
             matched_rmap = matching_map(a₊, a₋, ssmatcher)
             # this matcher has only one entry, so we use it to match
             # (we don't care what happens to the rest of the old_IDs, as the `rmap`
@@ -128,12 +132,8 @@ function ε_from_centroids(attractors::AbstractDict)
 end
 
 # group flows so that all old IDs that go to same new ID are in one vector
-function _grouped_flows(flows) # separated into
-    grouped = Dict{Int, Vector{Int}}()
-    oldids = collect(keys(flows))
-    for k in values(flows)
-        grouped[k] = findall(isequal(k), oldids)
-    end
+# i.e., new id => [prev ids that flowed to new id]
+function _grouped_flows(flows, allnewids) # separated into
+    grouped = Dict(newid=>[k for (k,v) in flows if v==newid] for newid in allnewids)
     return grouped
-    # return Dict(k => findall(isequal(k), flows) for k in values(flows))
 end
\ No newline at end of file
diff --git a/src/matching/basin_overlap.jl b/src/matching/basin_overlap.jl
index 6e8381e5..aa25f1e1 100644
--- a/src/matching/basin_overlap.jl
+++ b/src/matching/basin_overlap.jl
@@ -50,21 +50,19 @@ attraction and whose elements are the IDs.
 
 See [`MatchByBasinOverlap`](@ref) for how matching works.
 """
-function matching_map(b₊::AbstractArray, b₋::AbstractArray, matcher::MatchByBasinOverlap; i = nothing)
+function matching_map(b₊::AbstractArray, b₋::AbstractArray, matcher::MatchByBasinOverlap; kw...)
     a₊, a₋ = _basin_to_dict.((b₊, b₋))
-    matching_map(a₊, a₋, matcher; i)
+    matching_map(a₊, a₋, matcher; kw...)
 end
 
-function matching_map!(b₊::AbstractArray, b₋::AbstractArray, matcher::MatchByBasinOverlap; i = nothing)
-    rmap = matching_map(b₊, b₋, matcher; i)
+function matching_map!(b₊::AbstractArray, b₋::AbstractArray, matcher::MatchByBasinOverlap; kw...)
+    rmap = matching_map(b₊, b₋, matcher; kw...)
     replace!(b₊, rmap...)
     return rmap
 end
 
 # actual implementation
-function matching_map(a₊::AbstractDict, a₋, matcher::MatchByBasinOverlap;
-        i = nothing, kw...
-    )
+function matching_map(a₊::AbstractDict, a₋, matcher::MatchByBasinOverlap; kw...)
     # input checks
     if !(valtype(a₊) <: Vector{<:CartesianIndex})
         throw(ArgumentError("Incorrect input given. For matcher `MatchByBasinOverlap`,
@@ -99,5 +97,3 @@ function _basin_to_dict(b::AbstractArray{Int})
     d = Dict(k => findall(isequal(k), b) for k in ukeys)
     return d
 end
-
-# TODO: test that it works also with vector of basins of attraction
diff --git a/src/matching/matching_interface.jl b/src/matching/matching_interface.jl
index 758f9d79..e4b3f698 100644
--- a/src/matching/matching_interface.jl
+++ b/src/matching/matching_interface.jl
@@ -27,7 +27,7 @@ abstract type IDMatcher end
 """
     matching_map(
         a₊::Dict, a₋::Dict, matcher;
-        ds::DynamicalSystem, p, pprev, pidx, next_id
+        ds::DynamicalSystem, p, pprev, next_id
     ) → rmap
 
 Given dictionaries `a₊, a₋` mapping IDs to values,
@@ -47,14 +47,15 @@ Typically the +,- mean after and before some change of parameter of a dynamical
 ## Keyword arguments
 
 - `ds`: the dynamical system that generated `a₊, a₋`.
-- `p, pprev, pidx`: the parameter values corresponding to `a₊, a₋` and the index the
-  parameter has in `ds`.
+- `p, pprev`: the parameters corresponding to `a₊, a₋`. Both need to be iterables mapping
+  parameter index to parameter value (such as `Dict, Vector{Pair}`, etc., so whatever
+  can be given as input to `DynamicalSystems.set_parameters!`).
 - `next_id = next_free_id(a₊, a₋)`: the ID to give to values of  `a₊` that cannot be
   matched to `a₋` and hence must obtain a new unique ID.
 
-Some matchers like [`MatchBySSSetDistance`](@ref) do not utilize `ds, p, pprev, pidx` in any way
+Some matchers like [`MatchBySSSetDistance`](@ref) do not utilize `ds, p, pprev` in any way
 while other matchers like [`MatchByBasinEnclosure`](@ref) do, and those require
-expliticly giving values to `ds, p, pprev, pidx` as their default values
+expliticly giving values to `ds, p, pprev` as their default values
 is just `nothing`.
 """
 function matching_map(a₊, a₋, matcher::IDMatcher; kw...)
@@ -92,9 +93,10 @@ i.e., the pairs of `old => new` IDs.
 
 ## Keyword arguments
 
-- `prange = eachindex(attractors)`: the range of parameters from which to extract
-  the `p, pprev` values given to [`matching_map`](@ref).
-- `pidx, ds`: both propagated to [`matching_map`](@ref) and are `nothing` by default.
+- `pcurve = nothing`: the curve of parameters along which the continuation occured,
+  from which to extract the `p, pprev` values given to [`matching_map`](@ref).
+  See [`global_continuation`](@ref) if you are unsure what this means.
+- `ds = nothing`: propagated to [`matching_map`](@ref).
 - `retract_keys::Bool = true`: If `true` at the end the function will "retract" keys (i.e., make the
   integers smaller integers) so that all unique IDs
   are the 1-incremented positive integers. E.g., if the IDs where 1, 6, 8, they will become
@@ -159,27 +161,27 @@ end
 
 # Concrete implementation of `match_sequentially!`:
 function _rematch_ignored!(attractors_cont, matcher;
-        pidx = nothing, ds = nothing, prange = eachindex(attractors_cont),
+        ds = nothing, pcurve = eachindex(attractors_cont),
     )
     next_id = 1
     rmaps = Dict{keytype(attractors_cont[1]), keytype(attractors_cont[1])}[]
     for i in 1:length(attractors_cont)-1
         a₊, a₋ = attractors_cont[i+1], attractors_cont[i]
-        p, pprev = prange[i+1], prange[i]
+        p, pprev = pcurve[i+1], pcurve[i]
         # If there are no attractors, skip the matching
         (isempty(a₊) || isempty(a₋)) && continue
         # Here we always compute a next id. In this way, if an attractor disappears
         # and reappears, it will get a different (incremented) ID as it should!
         next_id_a = max(maximum(keys(a₊)), maximum(keys(a₋)))
         next_id = max(next_id, next_id_a) + 1
-        rmap = matching_map!(a₊, a₋, matcher; next_id, pidx, ds, p, pprev)
+        rmap = matching_map!(a₊, a₋, matcher; next_id, ds, p, pprev)
         push!(rmaps, rmap)
     end
     return rmaps
 end
 
 function _rematch_with_past!(attractors_cont, matcher;
-        pidx = nothing, ds = nothing, prange = eachindex(attractors_cont),
+        ds = nothing, pcurve = eachindex(attractors_cont),
     )
     # this dictionary stores all instances of previous attractors and is updated
     # at every step. It is then given to the matching function as if it was
@@ -188,12 +190,12 @@ function _rematch_with_past!(attractors_cont, matcher;
     rmaps = Dict{keytype(attractors_cont[1]), keytype(attractors_cont[1])}[]
     for i in 1:length(attractors_cont)-1
         a₊, a₋ = attractors_cont[i+1], attractors_cont[i]
-        p, pprev = prange[i+1], prange[i]
+        p, pprev = pcurve[i+1], pcurve[i]
         # update ghosts
         for (k, A) in a₋
             latest_ghosts[k] = A
         end
-        rmap = matching_map!(a₊, latest_ghosts, matcher; pprev, p, ds, pidx)
+        rmap = matching_map!(a₊, latest_ghosts, matcher; pprev, p, ds)
         push!(rmaps, rmap)
     end
     return rmaps
diff --git a/src/matching/sssdistance.jl b/src/matching/sssdistance.jl
index 580456f1..79fc7bf2 100644
--- a/src/matching/sssdistance.jl
+++ b/src/matching/sssdistance.jl
@@ -108,9 +108,5 @@ function _matching_map_distances(keys₊, keys₋, distances::Dict, threshold;
         end
     end
 
-    # the final step is to filter out equivalent mappings where key and value are the same
-    filter!(p -> p.first ≠ p.second, rmap)
     return rmap
 end
-
-
diff --git a/src/plotting.jl b/src/plotting.jl
index a3ccf4e5..a3cd0da4 100644
--- a/src/plotting.jl
+++ b/src/plotting.jl
@@ -61,7 +61,7 @@ export shaded_basins_heatmap, shaded_basins_heatmap!
 ##########################################################################################
 """
     animate_attractors_continuation(
-        ds::DynamicalSystem, attractors_cont, fractions_cont, prange, pidx;
+        ds::DynamicalSystem, attractors_cont, fractions_cont, pcurve;
         kwargs...
     )
 
@@ -72,7 +72,7 @@ and output of the [`global_continuation`](@ref) function into a video output.
 The input dynamical system `ds` is used to evolve initial conditions sampled from the
 found attractors, so that the attractors are better visualized.
 `attractors_cont, fractions_cont` are the output of [`global_continuation`](@ref)
-while `ds, prange, pidx` are the input to [`global_continuation`](@ref).
+while `ds, pcurve` are the input to [`global_continuation`](@ref).
 
 ## Keyword arguments
 
@@ -84,14 +84,15 @@ while `ds, prange, pidx` are the input to [`global_continuation`](@ref).
 - `figure, axis, fracaxis, legend`: named tuples propagated as keyword arguments to the
   creation of the `Figure`, the `Axis`, the "bar-like" axis containing the fractions,
   and the `axislegend` that adds the legend (if `add_legend = true`).
+- `add_legend = true`: whether to display the axis legend.
 """
 function animate_attractors_continuation end
 export animate_attractors_continuation
 
 """
-    plot_basins_curves(fractions_cont [, prange]; kwargs...)
+    plot_basins_curves(fractions_cont [, prange]; kw...)
 
-Plot the fractions of basins of attraction versus a parameter range,
+Plot the fractions of basins of attraction versus a parameter range/curve,
 i.e., visualize the output of [`global_continuation`](@ref).
 See also [`plot_basins_attractors_curves`](@ref) and
 [`plot_continuation_curves`](@ref).
@@ -112,10 +113,10 @@ function plot_basins_curves! end
 export plot_basins_curves, plot_basins_curves!
 
 """
-    plot_attractors_curves(attractors_cont, attractor_to_real, prange = 1:length(); kwargs...)
+    plot_attractors_curves(attractors_cont, attractor_to_real [, prange]; kw...)
 
-Same as in [`plot_basins_curves`](@ref) but visualizes the attractor dependence on
-the parameter instead of their fraction.
+Same as in [`plot_basins_curves`](@ref) but visualize the attractor dependence on
+the parameter(s) instead of their basin 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.
 See also [`plot_basins_attractors_curves`](@ref).
@@ -130,7 +131,7 @@ export plot_attractors_curves, plot_attractors_curves!
 """
     plot_continuation_curves(continuation_info [, prange]; kwargs...)
 
-Same as in [`plot_basins_curves`](@ref) but visualizes any arbitrary quantity characterizing
+Same as in [`plot_basins_curves`](@ref) but visualize any arbitrary quantity characterizing
 the continuation. Hence, the `continuation_info` is of exactly the same format as
 `fractions_cont`: a vector of dictionaries, each dictionary mapping attractor IDs to real numbers.
 `continuation_info` is meant to accompany `attractor_info` in [`plot_attractors_curves`](@ref).
diff --git a/test/continuation/matching_attractors.jl b/test/continuation/matching_attractors.jl
index a65fb00a..e15fe521 100644
--- a/test/continuation/matching_attractors.jl
+++ b/test/continuation/matching_attractors.jl
@@ -39,6 +39,22 @@ end
     end
 end
 
+@testset "synthetic multistable matching" begin
+    attractors_cont_simple = Dict{Int64, SVector{1, Float64}}[Dict(1 => [0.0]), Dict(1 => [0.0]), Dict(2 => [2.0], 1 => [0.0]), Dict(2 => [0.0], 1 => [2.0]), Dict(2 => [2.0], 3 => [4.0], 1 => [-5.0]), Dict(2 => [4.0], 3 => [-5.0], 1 => [2.0]), Dict(4 => [6.0], 2 => [0.0], 3 => [2.0], 1 => [4.0]), Dict(4 => [4.0], 2 => [0.0], 3 => [2.0], 1 => [6.0]), Dict(5 => [8.0], 4 => [6.0], 2 => [0.0], 3 => [2.0], 1 => [4.0])]
+    attractors_cont = [Dict(k=>StateSpaceSet(Vector(v)) for (k,v) in atts) for atts in attractors_cont_simple]
+
+    mapped_atts = deepcopy(attractors_cont)
+    rmaps = match_sequentially!(mapped_atts, default)
+
+    fractions_cont = [Dict(1 => 1.0), Dict(1 => 1.0), Dict(2 => 0.8, 1 => 0.2), Dict(2 => 0.2, 1 => 0.8), Dict(2 => 0.2, 3 => 0.6, 1 => 0.2), Dict(2 => 0.6, 3 => 0.2, 1 => 0.2), Dict(4 => 0.4, 2 => 0.2, 3 => 0.2, 1 => 0.2), Dict(4 => 0.2, 2 => 0.2, 3 => 0.2, 1 => 0.4), Dict(5 => 0.2, 4 => 0.2, 2 => 0.2, 3 => 0.2, 1 => 0.2)]
+    mapped_fracs = deepcopy(fractions_cont)
+    match_sequentially!(mapped_fracs, rmaps)
+
+    @test all(keys.(attractors_cont) .== keys.(mapped_atts)    )
+    @test all(keys.(attractors_cont) .== keys.(mapped_fracs)    )
+    @test all(Set.(values.(fractions_cont)) .== Set.(values.(mapped_fracs)))
+end
+
 @testset "global_continuation matching" begin
     # Make fake attractors with points that become more "separated" as "parameter"
     # is increased
@@ -164,4 +180,59 @@ end # Matcher by distance tests
     @test length(rmap) == 2
     @test rmap[2] == 3
     @test rmap[3] == 2
-end
\ No newline at end of file
+end
+
+@testset "BasinEncloure" begin
+
+    @testset "synthetic multistable continuation" begin
+        function dummy_multistable_equilibrium!(dx, x, p, n)
+            r = p[1]
+            num_max_atts = 5
+            x_max_right = 10
+            x_pos_atts = [x_max_right*(i-1)/num_max_atts for i=1:num_max_atts]
+            if 3 <=r < 4 && (x[1] <= 2)
+                x_pos_atts[1] = -5
+            end
+            num_atts = r < (num_max_atts + 1) ? floor(Int, r) : ceil(Int, 2*num_max_atts - r)
+            x_atts = [x_pos_atts[i] for i=1:num_atts]
+            att_of_x = findlast(xatt->xatt<=x[1], x_atts)
+            if x[1] < 0 att_of_x = 1 end
+            dx .= x_atts[att_of_x]
+
+            return nothing
+        end
+
+        ds = DeterministicIteratedMap(dummy_multistable_equilibrium!, [0.], [1.0])
+        featurizer(A,t) = A[end]
+        grouping_config = GroupViaPairwiseComparison(; threshold=0.2)
+        mapper = AttractorsViaFeaturizing(ds, featurizer, grouping_config)
+
+        xg = range(0, 10, length = 100)
+        grid = (xg,)
+        sampler, = statespace_sampler(grid, 1234)
+        samples_per_parameter = 1000
+        ics = Dataset([deepcopy(sampler()) for _ in 1:samples_per_parameter])
+
+        rrange = range(1, 9.5; step=0.5)
+        ridx = 1
+
+        mapper = AttractorsViaFeaturizing(ds, featurizer, grouping_config; T=10, Ttr=1)
+        matcher  = MatchByBasinEnclosure(;ε=0.1)
+        assc = AttractorSeedContinueMatch(mapper, matcher)
+        fs_curves, atts_all = global_continuation(assc, rrange, ridx, ics; show_progress = true)
+
+        atts_all_endpoint_solution = Dict{Int64, SVector{1, Float64}}[Dict(1 => [0.0]), Dict(1 => [0.0]), Dict(2 => [2.0], 1 => [0.0]), Dict(2 => [2.0], 1 => [0.0]), Dict(2 => [2.0], 3 => [4.0], 1 => [-5.0]), Dict(2 => [2.0], 3 => [4.0], 1 => [-5.0]), Dict(4 => [6.0], 2 => [2.0], 3 => [4.0], 1 => [0.0]), Dict(4 => [6.0], 2 => [2.0], 3 => [4.0], 1 => [0.0]), Dict(5 => [8.0], 4 => [6.0], 2 => [2.0], 3 => [4.0], 1 => [0.0]), Dict(5 => [8.0], 4 => [6.0], 2 => [2.0], 3 => [4.0], 1 => [0.0]), Dict(4 => [6.0], 2 => [2.0], 3 => [4.0], 1 => [0.0]), Dict(4 => [6.0], 2 => [2.0], 3 => [4.0], 1 => [0.0]), Dict(2 => [2.0], 3 => [4.0], 1 => [0.0]), Dict(2 => [2.0], 3 => [4.0], 1 => [0.0]), Dict(2 => [2.0], 1 => [0.0]), Dict(2 => [2.0], 1 => [0.0]), Dict(1 => [0.0]), Dict(1 => [0.0])]
+        atts_all_endpoint = [Dict(k=>v[end] for (k,v) in atts) for atts in atts_all]
+
+        @test atts_all_endpoint == atts_all_endpoint_solution
+
+        fs_curves_solution = [Dict(1 => 1.0), Dict(1 => 1.0), Dict(2 => 0.8, 1 => 0.2), Dict(2 => 0.8, 1 => 0.2), Dict(2 => 0.2, 3 => 0.6, 1 => 0.2), Dict(2 => 0.2, 3 => 0.6, 1 => 0.2), Dict(4 => 0.4, 2 => 0.2, 3 => 0.2, 1 => 0.2), Dict(4 => 0.4, 2 => 0.2, 3 => 0.2, 1 => 0.2), Dict(5 => 0.2, 4 => 0.2, 2 => 0.2, 3 => 0.2, 1 => 0.2), Dict(5 => 0.2, 4 => 0.2, 2 => 0.2, 3 => 0.2, 1 => 0.2), Dict(4 => 0.4, 2 => 0.2, 3 => 0.2, 1 => 0.2), Dict(4 => 0.4, 2 => 0.2, 3 => 0.2, 1 => 0.2), Dict(2 => 0.2, 3 => 0.6, 1 => 0.2), Dict(2 => 0.2, 3 => 0.6, 1 => 0.2), Dict(2 => 0.8, 1 => 0.2), Dict(2 => 0.8, 1 => 0.2), Dict(1 => 1.0), Dict(1 => 1.0)]
+        @test all(keys.(fs_curves_solution) .== keys.(fs_curves)    )
+        for (fs_curve, fs_curve_solution) in zip(fs_curves, fs_curves_solution)
+            for (k, fs) in fs_curve
+                @test isapprox(fs, fs_curve_solution[k], atol=1e-1)
+            end
+        end
+    end
+
+end
diff --git a/test/continuation/seed_continue_generic.jl b/test/continuation/seed_continue_generic.jl
index 3db26dec..664e27bb 100644
--- a/test/continuation/seed_continue_generic.jl
+++ b/test/continuation/seed_continue_generic.jl
@@ -9,7 +9,7 @@ using Random
     # for r > 0.5. It has analytically resolved fractions for any box.
     function dumb_map(dz, z, p, n)
         x, y = z
-        r = p[1]
+        r, q = p
 
         if r < 0.5
             dz[1] = dz[2] = 0.0
@@ -25,11 +25,11 @@ using Random
         return
     end
 
-    r = 1.0
-    ds = DeterministicIteratedMap(dumb_map, [0., 0.], [r])
+    r = 1.0; q = 0.5
+    ds = DeterministicIteratedMap(dumb_map, [0., 0.], [r, q])
     yg = xg = range(-10., 10, length = 100)
     grid = (xg,yg)
-    mapper = AttractorsViaRecurrences(ds, grid; sparse = true, show_progress = false)
+    mapper1 = AttractorsViaRecurrences(ds, grid; sparse = true, show_progress = false)
 
     sampler, = statespace_sampler(grid, 1234)
 
@@ -43,15 +43,29 @@ using Random
     # in all honesty, we don't have to test 2 grouping configs,
     # as the algorithm is agnostic to the grouping. But oh well!
     group1 = GroupViaClustering(optimal_radius_method = 0.1)
+    mapper2 = AttractorsViaFeaturizing(ds, featurizer, group1; Ttr = 2, T = 2)
+
     group2 = GroupViaPairwiseComparison(threshold = 0.1)
+    mapper3 = AttractorsViaFeaturizing(ds, featurizer, group2; Ttr = 2, T = 2)
+
+    mappers = [mapper1, mapper2, mapper3, mapper1]
 
-    @testset "grouping: $(nameof(typeof(group)))" for group in (group1, group2)
-        mapper = AttractorsViaFeaturizing(ds, featurizer, group; Ttr = 2, T = 2)
+    @testset "case: $(i)" for (i, mapper) in enumerate(mappers)
         algo = AttractorSeedContinueMatch(mapper)
-        fractions_cont, a = global_continuation(
-            algo, rrange, ridx, sampler;
-            show_progress = false, samples_per_parameter = 1000
-        )
+
+        if i < 4
+            fractions_cont, a = global_continuation(
+                algo, rrange, ridx, sampler;
+                show_progress = false, samples_per_parameter = 1000
+            )
+        else # test parameter curve version
+            pcurve = [[1 => r, 2 => 1.1] for r in rrange]
+            fractions_cont, a = global_continuation(
+                algo, pcurve, sampler;
+                show_progress = false, samples_per_parameter = 1000
+            )
+        end
+
 
         for (i, r) in enumerate(rrange)
 
diff --git a/test/mapping/attractor_mapping.jl b/test/mapping/attractor_mapping.jl
index 80c92b92..fab75ea9 100644
--- a/test/mapping/attractor_mapping.jl
+++ b/test/mapping/attractor_mapping.jl
@@ -115,6 +115,8 @@ function test_basins(ds, u0s, grid, expected_fs_raw, featurizer;
 
         config = GroupViaNearestFeature(templates; max_distance)
         mapper = AttractorsViaFeaturizing(ds, featurizer, config; Ttr=500)
+        # test the functionality mapper(u0) -> label
+        @test isinteger(mapper(get_state(ds))) == true
         test_basins_fractions(mapper; err = ferr, single_u_mapping = false)
     end