Multi- and Nonlocal- Continuation

In the tutorial on Bifurcation Diagrams we saw how one can create them by integrating ModelingToolkit.jl with BifurcationKit.jl. This approach is also often called continuation in the broader literature, because in essence we are "continuing" the location of individual un/stable fixed points or limit cycles in a dynamical system across a parameter axis.

Recently, an alternative continuation framework was proposed that takes a fundamentally different approach to continuation that is particularly suitable for complex systems. This framework is implemented in Attractors.jl as part of the DynamicalSystems.jl software library. This new continuation is called global continuation, while the one of BifurcationKit.jl is called local continuation.

Instead of continuing an individual fixed point or limit cycle, the global continuation finds all attractors of the dynamical system and continues all of them, in parallel, in a single continuation. It distinguishes and labels automatically the different attractors. Hence "multi-" for multiple attractors. Another key difference is that instead of estimating the local (or linear, or Jacobian) stability of the attractors, it estimates various measures of nonlocal stability (e.g, related with the size of the basins of attraction, or the size of a perturbation that would make the dynamical system state converge to an alternative attractor). Hence the "nonlocal-" component. More differences and pros & cons are discussed in the documentation of Attractors.jl.

Attractors and basins

This tutorial assumes that you have some familiarity with dynamical systems, specifically what are attractors and basins of attraction. If you don't have this yet, we recommend Chapter 1 of the textbook Nonlinear Dynamics.

Creating the `DynamicalSystem` via MTK

Let's showcase this framework by modelling a chaotic bistable dynamical system that we define via ModelingToolkit.jl, which will the be casted into a DynamicalSystem type for the DynamicalSystems.jl library. The equations of our system are

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

@variables x(t)=-4.0 y(t)=5.0 z(t)=0.0
@parameters a=5.0 b=0.1

eqs = [
    D(x) ~ y - x,
    D(y) ~ -x * z + b * abs(z),
    D(z) ~ x * y - a
]

\[ \begin{align} \frac{\mathrm{d} x\left( t \right)}{\mathrm{d}t} &= - x\left( t \right) + y\left( t \right) \\ \frac{\mathrm{d} y\left( t \right)}{\mathrm{d}t} &= b \left|z\left( t \right)\right| - x\left( t \right) z\left( t \right) \\ \frac{\mathrm{d} z\left( t \right)}{\mathrm{d}t} &= - a + x\left( t \right) y\left( t \right) \end{align} \]

Because our dynamical system is super simple, we will directly make an System and cast it in an ODEProblem as in the Systems tutorial. Since all state variables and parameters have a default value we can immediately write

@named modlorenz = System(eqs, t)
ssys = mtkcompile(modlorenz)
# The timespan given to the problem is irrelevant for DynamicalSystems.jl
prob = ODEProblem(ssys, [], (0.0, 1.0))

ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Initialization status: FULLY_DETERMINED
Non-trivial mass matrix: false
timespan: (0.0, 1.0)
u0: 3-element Vector{Float64}:
  0.0
  5.0
 -4.0

This prob can be turned to a dynamical system as simply as

using Attractors # or `DynamicalSystems`
ds = CoupledODEs(prob)

3-dimensional CoupledODEs
 deterministic: true
 discrete time: false
 in-place:      true
 dynamic rule:  GeneratedFunctionWrapper
 ODE solver:    Tsit5
 ODE kwargs:    (abstol = 1.0e-6, reltol = 1.0e-6)
 parameters:    MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}([5.0, 0.1], [0.0, -4.0, 0.0, 0.0, 0.0, 5.0], (), (), (), ())
 time:          0.0
 state:         [0.0, 5.0, -4.0]

DynamicalSystems.jl integrates fully with ModelingToolkit.jl out of the box and understands whether a given problem has been created via ModelingToolkit.jl. For example you can use the symbolic variables, or their Symbol representations, to access a system state or parameter

observe_state(ds, x)

-4.0

current_parameter(ds, :a) # or `a` directly

5.0

Finding all attractors in the state space

Attractors.jl provides an extensive interface for finding all (within a state space region and numerical accuracy) attractors of a dynamical system. This interface is structured around the type AttractorMapper and is discussed in the Attractors.jl documentation in detail. Here we will briefly mention one of the possible approaches, the recurrences-based algorithm. It finds attractors by finding locations in the state space where the trajectory returns again and again.

To use this technique, we first need to create a tessellation of the state space

grid = (
    range(-15.0, 15.0; length = 150), # x
    range(-20.0, 20.0; length = 150), # y
    range(-20.0, 20.0; length = 150) # z
)

(-15.0:0.20134228187919462:15.0, -20.0:0.2684563758389262:20.0, -20.0:0.2684563758389262:20.0)

which we then give as input to the AttractorsViaRecurrences mapper along with the dynamical system

mapper = AttractorsViaRecurrences(ds, grid;
    consecutive_recurrences = 1000,
    consecutive_lost_steps = 100
)

AttractorsViaRecurrences
 system:      CoupledODEs
 grid:        (-15.0:0.20134228187919462:15.0, -20.0:0.2684563758389262:20.0, -20.0:0.2684563758389262:20.0)
 attractors:  Dict{Int64, StateSpaceSets.StateSpaceSet{3, Float64, StaticArraysCore.SVector{3, Float64}}}()

to learn about the metaparameters of the algorithm visit the documentation of Attractors.jl.

This mapper object is incredibly powerful! It can be used to map initial conditions to attractor they converge to, while ensuring that initial conditions that converge to the same attractor are given the same label. For example, if we use the mapper as a function and give it an initial condition we get

mapper([-4.0, 5, 0])

mapper([4.0, 2, 0])

mapper([1.0, 3, 2])

The numbers returned are simply the unique identifiers of the attractors the initial conditions converged to.

DynamicalSystems.jl library is the only dynamical systems software (in any language) that provides such an infrastructure for mapping initial conditions of any arbitrary dynamical system to its unique attractors. And this is only the tip of this iceberg! The rest of the functionality of Attractors.jl is all full of brand new cutting edge progress in dynamical systems research.

The found attractors are stored in the mapper internally, to obtain them we use the function

attractors = extract_attractors(mapper)

Dict{Int64, StateSpaceSets.StateSpaceSet{3, Float64, StaticArraysCore.SVector{3, Float64}}} with 1 entry:
  1 => 3-dimensional StateSpaceSet{Float64} with 740 points

This is a dictionary that maps attractor IDs to the attractor sets themselves. StateSpaceSet is a wrapper of a vector of points and behaves exactly like a vector of points. We can plot them easily like

using CairoMakie
fig = Figure()
ax = Axis(fig[1, 1])
colors = ["#7143E0", "#191E44"]
for (id, A) in attractors
    scatter!(ax, A[:, [1, 3]]; color = colors[id])
end
fig

Basins of attraction

Estimating the basins of attraction of these attractors is a matter of a couple lines of code. First we define the state space are to estimate the basins for. Here we can re-use the grid we defined above. Then we only have to call

basins = basins_of_attraction(mapper, grid)

We won't run this in this tutorial because it is a length computation (150×150×150). We will however estimate a slice of the 3D basins of attraction. DynamicalSystems.jl allows for a rather straightforward setting of initial conditions:

ics = [Dict(:x => x, :y => 0, :z => z) for x in grid[1] for z in grid[3]]

22500-element Vector{Dict{Symbol, Real}}:
 Dict(:y => 0, :z => -20.0, :x => -15.0)
 Dict(:y => 0, :z => -19.731543624161073, :x => -15.0)
 Dict(:y => 0, :z => -19.463087248322147, :x => -15.0)
 Dict(:y => 0, :z => -19.19463087248322, :x => -15.0)
 Dict(:y => 0, :z => -18.926174496644297, :x => -15.0)
 Dict(:y => 0, :z => -18.65771812080537, :x => -15.0)
 Dict(:y => 0, :z => -18.389261744966444, :x => -15.0)
 Dict(:y => 0, :z => -18.120805369127517, :x => -15.0)
 Dict(:y => 0, :z => -17.85234899328859, :x => -15.0)
 Dict(:y => 0, :z => -17.583892617449663, :x => -15.0)
 ⋮
 Dict(:y => 0, :z => 17.85234899328859, :x => 15.0)
 Dict(:y => 0, :z => 18.120805369127517, :x => 15.0)
 Dict(:y => 0, :z => 18.389261744966444, :x => 15.0)
 Dict(:y => 0, :z => 18.65771812080537, :x => 15.0)
 Dict(:y => 0, :z => 18.926174496644297, :x => 15.0)
 Dict(:y => 0, :z => 19.19463087248322, :x => 15.0)
 Dict(:y => 0, :z => 19.463087248322147, :x => 15.0)
 Dict(:y => 0, :z => 19.731543624161073, :x => 15.0)
 Dict(:y => 0, :z => 20.0, :x => 15.0)

now we can estimate the basins of attraction on a slice on the x-z grid

fs, labels = basins_fractions(mapper, ics)
labels = reshape(labels, (length(grid[1]), length(grid[3])))

150×150 Matrix{Int64}:
 2  1  1  2  2  1  1  1  2  2  1  2  2  …  1  1  1  1  1  1  1  2  2  2  1  2
 1  1  2  2  1  1  2  2  1  2  1  1  1     1  1  1  1  1  1  1  1  1  2  2  1
 1  2  1  1  1  2  2  1  2  1  2  1  1     1  1  2  1  1  1  2  1  1  1  2  2
 2  1  1  1  2  1  1  1  2  1  1  1  1     2  2  2  2  1  2  1  1  1  1  1  2
 1  1  1  2  1  1  1  1  1  1  2  2  2     1  1  1  1  2  2  1  2  1  2  1  1
 1  1  2  2  2  2  1  1  2  2  1  1  1  …  2  1  1  1  1  1  2  2  1  1  1  1
 1  2  2  1  1  1  1  2  1  1  2  2  1     1  1  1  2  1  1  1  1  2  2  2  1
 2  1  1  1  1  2  1  1  2  1  1  1  2     2  1  1  2  1  1  2  1  1  1  2  1
 1  2  1  1  2  1  1  1  2  2  1  1  1     1  1  2  1  1  1  2  1  1  1  1  2
 2  2  1  2  1  1  1  1  1  1  2  1  1     2  2  1  1  2  1  1  1  1  2  1  1
 ⋮              ⋮              ⋮        ⋱        ⋮              ⋮           
 1  2  1  2  2  2  1  2  2  1  1  1  1     1  2  1  1  1  1  1  1  2  1  1  2
 1  1  1  1  2  1  1  1  2  1  1  1  2     2  1  2  1  2  1  1  2  1  1  1  2
 2  1  2  1  2  2  1  1  1  2  1  1  2     2  2  1  2  2  1  2  2  1  1  1  1
 2  1  2  1  1  2  1  1  1  1  2  1  2     1  2  2  1  1  2  2  1  1  1  1  1
 1  1  1  2  1  2  2  1  1  1  1  1  2  …  1  1  1  1  2  1  2  1  1  1  2  1
 1  2  1  2  1  1  2  1  1  1  1  1  1     1  1  2  2  2  1  1  1  1  2  1  1
 1  1  1  1  1  1  1  2  1  2  2  1  1     2  2  1  2  1  1  1  1  1  1  1  1
 2  1  2  1  2  2  1  1  2  2  1  1  2     1  1  2  1  1  1  1  1  2  2  1  1
 2  1  2  1  1  2  1  1  1  2  2  2  1     1  1  1  1  1  1  2  1  1  1  2  1

and visualize them

heatmap(grid[1], grid[3], labels; colormap = colors)

Global continuation

We've already outlined the principles of the global continuation, so let's just do it here! We first have to define a global continuation algorithm, which for this tutorial, it is just a wrapper of the existing mapper

ascm = AttractorSeedContinueMatch(mapper);

Attractors.AttractorSeedContinueMatch{Attractors.AttractorsViaRecurrences{DynamicalSystemsBase.CoupledODEs{true, 3, OrdinaryDiffEqCore.ODEIntegrator{OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, true, Vector{Float64}, Nothing, Float64, MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, SciMLBase.ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x40467422, 0x83ad5f21, 0x91b69766, 0xe66c70f3, 0xbcaccff9), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x0ebd4a54, 0x0f321c22, 0x031a510c, 0xf9a67687, 0x63869264), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{System}, Nothing, System, SciMLBase.OverrideInitData{NonlinearProblem{Nothing, true, MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Float64, Vector{Float64}}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, NonlinearFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xef581217, 0x116a5f2b, 0x15969938, 0xb58b247f, 0x6a5369a7), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xb07d40d3, 0x532f3061, 0xac2c47c7, 0xf5e48a4c, 0xe2a3b5a6), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{System}, Nothing, System, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, typeof(ModelingToolkit.update_initializeprob!), ComposedFunction{ComposedFunction{typeof(identity), typeof(ModelingToolkit.safe_float)}, SymbolicIndexingInterface.TimeIndependentObservedFunction{ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x0b66872e, 0xe4dc070b, 0xe5757108, 0x6482ad97, 0x9aff3e7b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x15b321ef, 0x1ac8f62a, 0xc5fc3662, 0xb920011d, 0x9a65042d), Nothing}}}}, ModelingToolkit.var"#initprobpmap_split#908"{ModelingToolkit.var"#_getter#904"{Tuple{ComposedFunction{ModelingToolkit.PConstructorApplicator{typeof(identity)}, ModelingToolkit.ObservedWrapper{false, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9c8936e5, 0x80a1ca0c, 0x496814e1, 0x6b034fc6, 0xd86b44e3), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x35e6d268, 0x0743b4ee, 0x6c68d9cd, 0xfb96b0c3, 0xac1d716d), Nothing}}}}, ComposedFunction{ModelingToolkit.PConstructorApplicator{typeof(identity)}, ModelingToolkit.ObservedWrapper{false, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x1b88a256, 0xff61e1be, 0xbfc2d27a, 0xe4ef436d, 0xf574bccd), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e418c9c, 0x000f0740, 0xd3cd8014, 0x1cbae35d, 0x79b38940), Nothing}}}}, Returns{Tuple{}}, Returns{Tuple{}}, Returns{Tuple{}}}}}, ModelingToolkit.InitializationMetadata{ModelingToolkit.ReconstructInitializeprob{ModelingToolkit.var"#_getter#904"{Tuple{ComposedFunction{ModelingToolkit.PConstructorApplicator{typeof(identity)}, ModelingToolkit.ObservedWrapper{true, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xb755fc84, 0xd201d84b, 0x3b0e1ddf, 0x6096a157, 0xad41e8b5), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xce77fac4, 0x751ccd0e, 0x8ce6e574, 0xb6c043c9, 0x157f1c52), Nothing}}}}, Returns{StaticArraysCore.SizedVector{0, Float64, Vector{Float64}}}, Returns{Tuple{}}, Returns{Tuple{}}, Returns{Tuple{}}}}, ComposedFunction{typeof(identity), SymbolicIndexingInterface.MultipleGetters{SymbolicIndexingInterface.ContinuousTimeseries, Vector{Any}}}}, ModelingToolkit.GetUpdatedU0{SymbolicIndexingInterface.MultipleGetters{SymbolicIndexingInterface.ContinuousTimeseries, Vector{SymbolicUtils.BasicSymbolic{Real}}}, SymbolicIndexingInterface.MultipleParametersGetter{SymbolicIndexingInterface.IndexerNotTimeseries, Vector{SymbolicIndexingInterface.GetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Initials, Int64}}}, Nothing}}, ModelingToolkit.SetInitialUnknowns{SymbolicIndexingInterface.MultipleSetters{Vector{SymbolicIndexingInterface.ParameterHookWrapper{SymbolicIndexingInterface.SetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Initials, Int64}}, SymbolicUtils.BasicSymbolic{Real}}}}}}, Val{true}}, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardODEProblem}, OrdinaryDiffEqTsit5.Tsit5{typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.InterpolationData{ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x40467422, 0x83ad5f21, 0x91b69766, 0xe66c70f3, 0xbcaccff9), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x0ebd4a54, 0x0f321c22, 0x031a510c, 0xf9a67687, 0x63869264), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{System}, Nothing, System, SciMLBase.OverrideInitData{NonlinearProblem{Nothing, true, MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Float64, Vector{Float64}}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, NonlinearFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xef581217, 0x116a5f2b, 0x15969938, 0xb58b247f, 0x6a5369a7), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xb07d40d3, 0x532f3061, 0xac2c47c7, 0xf5e48a4c, 0xe2a3b5a6), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{System}, Nothing, System, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, typeof(ModelingToolkit.update_initializeprob!), ComposedFunction{ComposedFunction{typeof(identity), typeof(ModelingToolkit.safe_float)}, SymbolicIndexingInterface.TimeIndependentObservedFunction{ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x0b66872e, 0xe4dc070b, 0xe5757108, 0x6482ad97, 0x9aff3e7b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x15b321ef, 0x1ac8f62a, 0xc5fc3662, 0xb920011d, 0x9a65042d), Nothing}}}}, ModelingToolkit.var"#initprobpmap_split#908"{ModelingToolkit.var"#_getter#904"{Tuple{ComposedFunction{ModelingToolkit.PConstructorApplicator{typeof(identity)}, ModelingToolkit.ObservedWrapper{false, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9c8936e5, 0x80a1ca0c, 0x496814e1, 0x6b034fc6, 0xd86b44e3), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x35e6d268, 0x0743b4ee, 0x6c68d9cd, 0xfb96b0c3, 0xac1d716d), Nothing}}}}, ComposedFunction{ModelingToolkit.PConstructorApplicator{typeof(identity)}, ModelingToolkit.ObservedWrapper{false, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x1b88a256, 0xff61e1be, 0xbfc2d27a, 0xe4ef436d, 0xf574bccd), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e418c9c, 0x000f0740, 0xd3cd8014, 0x1cbae35d, 0x79b38940), Nothing}}}}, Returns{Tuple{}}, Returns{Tuple{}}, Returns{Tuple{}}}}}, ModelingToolkit.InitializationMetadata{ModelingToolkit.ReconstructInitializeprob{ModelingToolkit.var"#_getter#904"{Tuple{ComposedFunction{ModelingToolkit.PConstructorApplicator{typeof(identity)}, ModelingToolkit.ObservedWrapper{true, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xb755fc84, 0xd201d84b, 0x3b0e1ddf, 0x6096a157, 0xad41e8b5), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xce77fac4, 0x751ccd0e, 0x8ce6e574, 0xb6c043c9, 0x157f1c52), Nothing}}}}, Returns{StaticArraysCore.SizedVector{0, Float64, Vector{Float64}}}, Returns{Tuple{}}, Returns{Tuple{}}, Returns{Tuple{}}}}, ComposedFunction{typeof(identity), SymbolicIndexingInterface.MultipleGetters{SymbolicIndexingInterface.ContinuousTimeseries, Vector{Any}}}}, ModelingToolkit.GetUpdatedU0{SymbolicIndexingInterface.MultipleGetters{SymbolicIndexingInterface.ContinuousTimeseries, Vector{SymbolicUtils.BasicSymbolic{Real}}}, SymbolicIndexingInterface.MultipleParametersGetter{SymbolicIndexingInterface.IndexerNotTimeseries, Vector{SymbolicIndexingInterface.GetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Initials, Int64}}}, Nothing}}, ModelingToolkit.SetInitialUnknowns{SymbolicIndexingInterface.MultipleSetters{Vector{SymbolicIndexingInterface.ParameterHookWrapper{SymbolicIndexingInterface.SetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Initials, Int64}}, SymbolicUtils.BasicSymbolic{Real}}}}}}, Val{true}}, Nothing}, Vector{Vector{Float64}}, Vector{Float64}, Vector{Vector{Vector{Float64}}}, Nothing, OrdinaryDiffEqTsit5.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, Nothing}, SciMLBase.DEStats, Nothing, Nothing, Nothing, Nothing}, ODEFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x40467422, 0x83ad5f21, 0x91b69766, 0xe66c70f3, 0xbcaccff9), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x0ebd4a54, 0x0f321c22, 0x031a510c, 0xf9a67687, 0x63869264), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{System}, Nothing, System, SciMLBase.OverrideInitData{NonlinearProblem{Nothing, true, MTKParameters{Vector{Float64}, StaticArraysCore.SizedVector{0, Float64, Vector{Float64}}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}, NonlinearFunction{true, SciMLBase.FullSpecialize, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xef581217, 0x116a5f2b, 0x15969938, 0xb58b247f, 0x6a5369a7), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xb07d40d3, 0x532f3061, 0xac2c47c7, 0xf5e48a4c, 0xe2a3b5a6), Nothing}}, LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, ModelingToolkit.ObservedFunctionCache{System}, Nothing, System, Nothing, Nothing}, Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}, SciMLBase.StandardNonlinearProblem}, typeof(ModelingToolkit.update_initializeprob!), ComposedFunction{ComposedFunction{typeof(identity), typeof(ModelingToolkit.safe_float)}, SymbolicIndexingInterface.TimeIndependentObservedFunction{ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x0b66872e, 0xe4dc070b, 0xe5757108, 0x6482ad97, 0x9aff3e7b), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x15b321ef, 0x1ac8f62a, 0xc5fc3662, 0xb920011d, 0x9a65042d), Nothing}}}}, ModelingToolkit.var"#initprobpmap_split#908"{ModelingToolkit.var"#_getter#904"{Tuple{ComposedFunction{ModelingToolkit.PConstructorApplicator{typeof(identity)}, ModelingToolkit.ObservedWrapper{false, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x9c8936e5, 0x80a1ca0c, 0x496814e1, 0x6b034fc6, 0xd86b44e3), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x35e6d268, 0x0743b4ee, 0x6c68d9cd, 0xfb96b0c3, 0xac1d716d), Nothing}}}}, ComposedFunction{ModelingToolkit.PConstructorApplicator{typeof(identity)}, ModelingToolkit.ObservedWrapper{false, ModelingToolkit.GeneratedFunctionWrapper{(2, 2, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x1b88a256, 0xff61e1be, 0xbfc2d27a, 0xe4ef436d, 0xf574bccd), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0x6e418c9c, 0x000f0740, 0xd3cd8014, 0x1cbae35d, 0x79b38940), Nothing}}}}, Returns{Tuple{}}, Returns{Tuple{}}, Returns{Tuple{}}}}}, ModelingToolkit.InitializationMetadata{ModelingToolkit.ReconstructInitializeprob{ModelingToolkit.var"#_getter#904"{Tuple{ComposedFunction{ModelingToolkit.PConstructorApplicator{typeof(identity)}, ModelingToolkit.ObservedWrapper{true, ModelingToolkit.GeneratedFunctionWrapper{(2, 3, true), RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xb755fc84, 0xd201d84b, 0x3b0e1ddf, 0x6096a157, 0xad41e8b5), Nothing}, RuntimeGeneratedFunctions.RuntimeGeneratedFunction{(:ˍ₋out, :__mtk_arg_1, :___mtkparameters___, :t), ModelingToolkit.var"#_RGF_ModTag", ModelingToolkit.var"#_RGF_ModTag", (0xce77fac4, 0x751ccd0e, 0x8ce6e574, 0xb6c043c9, 0x157f1c52), Nothing}}}}, Returns{StaticArraysCore.SizedVector{0, Float64, Vector{Float64}}}, Returns{Tuple{}}, Returns{Tuple{}}, Returns{Tuple{}}}}, ComposedFunction{typeof(identity), SymbolicIndexingInterface.MultipleGetters{SymbolicIndexingInterface.ContinuousTimeseries, Vector{Any}}}}, ModelingToolkit.GetUpdatedU0{SymbolicIndexingInterface.MultipleGetters{SymbolicIndexingInterface.ContinuousTimeseries, Vector{SymbolicUtils.BasicSymbolic{Real}}}, SymbolicIndexingInterface.MultipleParametersGetter{SymbolicIndexingInterface.IndexerNotTimeseries, Vector{SymbolicIndexingInterface.GetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Initials, Int64}}}, Nothing}}, ModelingToolkit.SetInitialUnknowns{SymbolicIndexingInterface.MultipleSetters{Vector{SymbolicIndexingInterface.ParameterHookWrapper{SymbolicIndexingInterface.SetParameterIndex{ModelingToolkit.ParameterIndex{SciMLStructures.Initials, Int64}}, SymbolicUtils.BasicSymbolic{Real}}}}}}, Val{true}}, Nothing}, OrdinaryDiffEqTsit5.Tsit5Cache{Vector{Float64}, Vector{Float64}, Vector{Float64}, typeof(OrdinaryDiffEqCore.trivial_limiter!), typeof(OrdinaryDiffEqCore.trivial_limiter!), Static.False}, OrdinaryDiffEqCore.DEOptions{Float64, Float64, Float64, Float64, OrdinaryDiffEqCore.PIController{Rational{Int64}}, typeof(DiffEqBase.ODE_DEFAULT_NORM), typeof(LinearAlgebra.opnorm), Bool, SciMLBase.CallbackSet{Tuple{}, Tuple{}}, typeof(DiffEqBase.ODE_DEFAULT_ISOUTOFDOMAIN), typeof(DiffEqBase.ODE_DEFAULT_PROG_MESSAGE), typeof(DiffEqBase.ODE_DEFAULT_UNSTABLE_CHECK), DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, DataStructures.BinaryHeap{Float64, DataStructures.FasterForward}, Nothing, Nothing, Float64, Tuple{}, Tuple{}, Tuple{}}, Vector{Float64}, Float64, Nothing, OrdinaryDiffEqCore.DefaultInit, Nothing}, MTKParameters{Vector{Float64}, Vector{Float64}, Tuple{}, Tuple{}, Tuple{}, Tuple{}}}, Attractors.BasinsInfo{3, Attractors.RegularGrid{3, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Float64, Float64, StaticArraysCore.SVector{3, Float64}, Attractors.SparseArray{Int64, 3}}, Attractors.RegularGrid{3, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Base.Pairs{Symbol, Int64, Tuple{Symbol, Symbol}, @NamedTuple{consecutive_recurrences::Int64, consecutive_lost_steps::Int64}}}, Attractors.MatchBySSSetDistance{StateSpaceSets.Centroid{Distances.Euclidean}, Float64}, typeof(Attractors._default_seeding)}(AttractorsViaRecurrences
 system:      CoupledODEs
 grid:        (-15.0:0.20134228187919462:15.0, -20.0:0.2684563758389262:20.0, -20.0:0.2684563758389262:20.0)
 attractors:  Dict{Int64, StateSpaceSets.StateSpaceSet{3, Float64, StaticArraysCore.SVector{3, Float64}}}(2 => 3-dimensional StateSpaceSet{Float64} with 232 points, 1 => 3-dimensional StateSpaceSet{Float64} with 740 points)
, Attractors.MatchBySSSetDistance{StateSpaceSets.Centroid{Distances.Euclidean}, Float64}(StateSpaceSets.Centroid{Distances.Euclidean}(Distances.Euclidean(0.0)), Inf, false), Attractors._default_seeding)

we need two more ingredients to perform the global continuation. One is a sampler of initial conditions in the state space. Here we'll uniformly sample initial conditions within this grid we have already defined

sampler, = statespace_sampler(grid);

(StateSpaceSets.RectangleGenerator{Float64, StaticArraysCore.SVector{3, Float64}, Random.Xoshiro}([-15.0, -20.0, -20.0], [30.0, 40.0, 40.0], [[0.0, 0.0, 0.0]], Random.Xoshiro(0xdfbcce9066bc53ec, 0x46b8c09ac669cc3b, 0x0ac5a73e7b31ec4d, 0x7b4adf67b18c20da, 0x48cf706efcc83c14)), StateSpaceSets.var"#isinside#61"{StaticArraysCore.SVector{3, Float64}, StaticArraysCore.SVector{3, Float64}}([15.0, 20.0, 20.0], [-15.0, -20.0, -20.0]))

the last ingredient is what parameter(s) to perform the continuation over. In contrast to local continuation, where we can only specify a parameter range, in global continuation one can specify an exact parameter curve to continue over. This curve can span any-dimensional parameter space, in contrast to the 1D or 2D parameter spaces supported in local continuation. Here we will use the curve

params(θ) = [:a => 5 + 0.5cos(θ), :b => 0.1 + 0.01sin(θ)]
θs = range(0, 2π; length = 101)
pcurve = params.(θs)

101-element Vector{Vector{Pair{Symbol, Float64}}}:
 [:a => 5.5, :b => 0.1]
 [:a => 5.499013364214136, :b => 0.10062790519529313]
 [:a => 5.496057350657239, :b => 0.10125333233564304]
 [:a => 5.491143625364344, :b => 0.10187381314585725]
 [:a => 5.4842915805643155, :b => 0.10248689887164855]
 [:a => 5.475528258147577, :b => 0.10309016994374948]
 [:a => 5.4648882429441255, :b => 0.10368124552684678]
 [:a => 5.45241352623301, :b => 0.10425779291565074]
 [:a => 5.438153340021932, :b => 0.10481753674101715]
 [:a => 5.422163962751007, :b => 0.10535826794978997]
 ⋮
 [:a => 5.438153340021932, :b => 0.09518246325898286]
 [:a => 5.45241352623301, :b => 0.09574220708434927]
 [:a => 5.4648882429441255, :b => 0.09631875447315323]
 [:a => 5.475528258147577, :b => 0.09690983005625053]
 [:a => 5.4842915805643155, :b => 0.09751310112835145]
 [:a => 5.491143625364344, :b => 0.09812618685414276]
 [:a => 5.496057350657239, :b => 0.09874666766435695]
 [:a => 5.499013364214136, :b => 0.09937209480470688]
 [:a => 5.5, :b => 0.1]

which makes an ellipsis over the parameter space.

We put these three ingredients together to call the global continuation

fractions_cont, attractors_cont = global_continuation(ascm, pcurve, sampler);

([Dict(2 => 0.25, 1 => 0.75), Dict(2 => 0.21568627450980393, 1 => 0.7843137254901961), Dict(2 => 0.2647058823529412, 1 => 0.7352941176470589), Dict(2 => 0.28431372549019607, 1 => 0.7156862745098039), Dict(2 => 0.3235294117647059, 1 => 0.6764705882352942), Dict(2 => 0.2549019607843137, 1 => 0.7450980392156863), Dict(2 => 0.28431372549019607, 1 => 0.7156862745098039), Dict(2 => 0.28431372549019607, 1 => 0.7156862745098039), Dict(2 => 0.22549019607843138, 1 => 0.7745098039215687), Dict(2 => 0.3333333333333333, 1 => 0.6666666666666666)  …  Dict(2 => 0.35294117647058826, 3 => 0.6470588235294118), Dict(2 => 0.27450980392156865, 3 => 0.7254901960784313), Dict(2 => 0.21568627450980393, 3 => 0.7843137254901961), Dict(2 => 0.29411764705882354, 3 => 0.7058823529411765), Dict(2 => 0.21568627450980393, 3 => 0.7843137254901961), Dict(2 => 0.29411764705882354, 3 => 0.7058823529411765), Dict(2 => 0.27450980392156865, 3 => 0.7254901960784313), Dict(2 => 0.30392156862745096, 3 => 0.696078431372549), Dict(2 => 0.27450980392156865, 3 => 0.7254901960784313), Dict(2 => 0.3235294117647059, 3 => 0.6764705882352942)], Dict{Int64, StateSpaceSets.StateSpaceSet{3, Float64, StaticArraysCore.SVector{3, Float64}}}[Dict(2 => 3-dimensional StateSpaceSet{Float64} with 276 points, 1 => 3-dimensional StateSpaceSet{Float64} with 438 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 287 points, 1 => 3-dimensional StateSpaceSet{Float64} with 442 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 279 points, 1 => 3-dimensional StateSpaceSet{Float64} with 480 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 64 points, 1 => 3-dimensional StateSpaceSet{Float64} with 524 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 286 points, 1 => 3-dimensional StateSpaceSet{Float64} with 596 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 262 points, 1 => 3-dimensional StateSpaceSet{Float64} with 603 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 259 points, 1 => 3-dimensional StateSpaceSet{Float64} with 542 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 271 points, 1 => 3-dimensional StateSpaceSet{Float64} with 548 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 254 points, 1 => 3-dimensional StateSpaceSet{Float64} with 600 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 278 points, 1 => 3-dimensional StateSpaceSet{Float64} with 640 points)  …  Dict(2 => 3-dimensional StateSpaceSet{Float64} with 261 points, 3 => 3-dimensional StateSpaceSet{Float64} with 577 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 263 points, 3 => 3-dimensional StateSpaceSet{Float64} with 564 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 244 points, 3 => 3-dimensional StateSpaceSet{Float64} with 526 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 247 points, 3 => 3-dimensional StateSpaceSet{Float64} with 521 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 275 points, 3 => 3-dimensional StateSpaceSet{Float64} with 490 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 244 points, 3 => 3-dimensional StateSpaceSet{Float64} with 427 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 275 points, 3 => 3-dimensional StateSpaceSet{Float64} with 441 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 276 points, 3 => 3-dimensional StateSpaceSet{Float64} with 438 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 263 points, 3 => 3-dimensional StateSpaceSet{Float64} with 426 points), Dict(2 => 3-dimensional StateSpaceSet{Float64} with 276 points, 3 => 3-dimensional StateSpaceSet{Float64} with 451 points)])

The output of the continuation is how the attractors and their basins fractions change over this parameter curve. We can visualize this directly using a convenience function

fig = plot_basins_attractors_curves(
    fractions_cont, attractors_cont, A -> minimum(A[:, 1]), θs
)

The top panel shows the relative basins of attractions of the attractors and the bottom panel shows their minimum x-position. The colors correspond to unique attractors. Perhaps making a video is easier to understand:

animate_attractors_continuation(
    ds, attractors_cont, fractions_cont, pcurve;
    savename = "curvecont.mp4"
);

To learn more about this global continuation and its various options, and more details about how it compares with local continuation, visit the documentation of Attractors.jl.

Multi- and Nonlocal- Continuation

Creating the DynamicalSystem via MTK

Finding all attractors in the state space

Basins of attraction

Global continuation

Creating the `DynamicalSystem` via MTK