From 4538692ba8bf97ce6a00b748cb7a3f5c4247614d Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Fri, 7 Mar 2025 20:28:20 +0000
Subject: [PATCH 01/12] Informative error on nonexisting symbolic parameter

---
 CHANGELOG.md                          |  4 ++++
 Project.toml                          |  2 +-
 src/core/dynamicalsystem_interface.jl | 10 ++++++++--
 test/mtk_integration.jl               |  8 ++++++++
 4 files changed, 21 insertions(+), 3 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index abead5e7..45153690 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,7 @@
+# v3.14.0
+
+- `set_parameter!` and `current_parameter` will now throw an informative error message explicitly naming the parameter if the user tries to get/set a symbolic parameter that does not exist in the MTK-generated dynamical system.
+
 # v3.13.0
 
 - `jacobian` will now re-use symbolically generated Jacobian functions if the dynamical system is made through a ModelingToolkit.jl-generated `DEProblem`.
diff --git a/Project.toml b/Project.toml
index a15e627f..47d6de69 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "DynamicalSystemsBase"
 uuid = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 repo = "https://github.com/JuliaDynamics/DynamicalSystemsBase.jl.git"
-version = "3.13.2"
+version = "3.14.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
diff --git a/src/core/dynamicalsystem_interface.jl b/src/core/dynamicalsystem_interface.jl
index 95b26102..0290d446 100644
--- a/src/core/dynamicalsystem_interface.jl
+++ b/src/core/dynamicalsystem_interface.jl
@@ -247,7 +247,10 @@ function current_parameter(ds::DynamicalSystem, index, p = current_parameters(ds
     prob = referrenced_sciml_prob(ds)
     if !has_referrenced_model(prob)
         return _get_parameter(p, index)
-    else # symbolic dispatch
+    else # symbolic parameter
+        if !SymbolicIndexingInterface.is_parameter(prob, index)
+            error("Symbolic parameter with name $(index) does not exist in the system.")
+        end
         i = SymbolicIndexingInterface.getp(prob, index)
         return i(p)
     end
@@ -431,7 +434,10 @@ function _set_parameter!(ds::DynamicalSystem, index, value, p = current_paramete
         else
             setproperty!(p, index, value)
         end
-    else
+    else # symbolic parameter
+        if !SymbolicIndexingInterface.is_parameter(prob, index)
+            error("Symbolic parameter with name $(index) does not exist in the system.")
+        end
         set! = SymbolicIndexingInterface.setp(prob, index)
         set!(p, value)
     end
diff --git a/test/mtk_integration.jl b/test/mtk_integration.jl
index 629316cc..00e888b2 100644
--- a/test/mtk_integration.jl
+++ b/test/mtk_integration.jl
@@ -200,6 +200,14 @@ end
 
 end
 
+@testset "informative errors" begin
+    prob = ODEProblem(roessler_model)
+    ds = CoupledODEs(prob)
+    @parameters XOXO = 0.5
+    @test_throws "XOXO" current_parameter(ds, :XOXO)
+    @test_throws "XOXO" current_parameter(ds, XOXO)
+end
+
 # %% Trajectory with mixed and time dependent indexing
 η2 = 1.0
 η3 = 0.3

From 71ce2b1b3fff18dce5460a5ea4057d1bc7b9e9ee Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Fri, 7 Mar 2025 20:43:18 +0000
Subject: [PATCH 02/12] mention init cond as parameters

---
 src/core/dynamicalsystem_interface.jl | 18 ++++++++----------
 1 file changed, 8 insertions(+), 10 deletions(-)

diff --git a/src/core/dynamicalsystem_interface.jl b/src/core/dynamicalsystem_interface.jl
index 0290d446..c66de2f2 100644
--- a/src/core/dynamicalsystem_interface.jl
+++ b/src/core/dynamicalsystem_interface.jl
@@ -56,14 +56,12 @@ The referenced MTK model corresponding to the dynamical system can be obtained w
 
 See also the DynamicalSystems.jl tutorial online for an example.
 
-!!! warn "ModelingToolkit.jl v9"
-    In ModelingToolkit.jl v9 the default `split` behavior of the parameter container
-    is `true`. This means that the parameter container is no longer a `Vector{Float64}`
-    by default, which means that you cannot use integers to access parameters.
-    It is recommended to keep `split = true` (default) and only access
-    parameters via their symbolic parameter binding.
-    Use `structural_simplify(sys; split = false)` to allow accessing parameters
-    with integers again.
+!!! warn "Initial conditions and parameters"
+    ModelingToolkit.jl treats initial conditions of all variables as additional parameters.
+    This is because its terminology is aimed primarily at initial value problems rather
+    than dynamical systems. It is strongly recommended when making a dynamical system from
+    an MTK model to only access parameters via symbols, and to not use the `split = false`
+    keyword when creating the problem type from the MTK model.
 
 ## API
 
@@ -79,8 +77,8 @@ unless when developing new algorithm implementations that use dynamical systems.
 ### API - obtain information
 
 - `ds(t)` with `ds` an instance of `DynamicalSystem`: return the state of `ds` at time `t`.
-  For continuous time systems this interpolates and extrapolates,
-  while for discrete time systems it only works if `t` is the current time.
+  For continuous time systems this interpolates current and previous time step and is very accurate;
+  for discrete time systems it only works if `t` is the current time.
 - [`current_state`](@ref)
 - [`initial_state`](@ref)
 - [`observe_state`](@ref)

From b6a3e20416d36011b8253574f1014cbe91fc5eb1 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Fri, 7 Mar 2025 20:45:16 +0000
Subject: [PATCH 03/12] correct typo

---
 src/core/dynamicalsystem_interface.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/core/dynamicalsystem_interface.jl b/src/core/dynamicalsystem_interface.jl
index c66de2f2..9b34df9f 100644
--- a/src/core/dynamicalsystem_interface.jl
+++ b/src/core/dynamicalsystem_interface.jl
@@ -192,7 +192,7 @@ Return the state `u` of `ds` _observed_ at "index" `i`. Possibilities are:
   symbolic variables.
 
 For [`ProjectedDynamicalSystem`](@ref), this function assumes that the
-state of the system is the full state space state, not the projected one
+state of the system is the original state space state, not the projected one
 (this makes the most sense for allowing MTK-based indexing).
 
 Use [`state_name`](@ref) for an accompanying name.

From ca161e5ede320b8f25d239875feb7843e3485fd0 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Mon, 10 Mar 2025 12:50:55 +0000
Subject: [PATCH 04/12] add a new function for named current parameters

---
 src/core/dynamicalsystem_interface.jl | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/src/core/dynamicalsystem_interface.jl b/src/core/dynamicalsystem_interface.jl
index 9b34df9f..22ee8500 100644
--- a/src/core/dynamicalsystem_interface.jl
+++ b/src/core/dynamicalsystem_interface.jl
@@ -329,6 +329,17 @@ successful_step(ds::DiscreteTimeDynamicalSystem) = all(x -> (isfinite(x) && !isn
 # Generic implementation, most types re-define it as compile-time info
 StateSpaceSets.dimension(ds::DynamicalSystem) = length(current_state(ds))
 
+# TODO: Document it
+function named_current_parameters(ds::DynamicalSystem)
+    mtk = referrenced_sciml_model(ebm)
+    # check if I can use `SymbolicIndexingInterface.parameter_symbols(mtk)` instead.
+    # Does it include initials or not?
+    params_names = Symbol.(ModelingToolkit.parameters(mtk))
+    params_values = current_parameter.(ds, params_names)
+    params = Dict(params_names .=> params_values)
+    return params
+end
+
 ###########################################################################################
 # API - altering status of the system
 ###########################################################################################

From 7066d655263a0f87e5bf0f253913dd4ac9063b3d Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Mon, 10 Mar 2025 20:14:22 +0000
Subject: [PATCH 05/12] organize MTK tests better

---
 test/mtk_integration.jl | 166 +++++++++++++++++++++-------------------
 1 file changed, 87 insertions(+), 79 deletions(-)

diff --git a/test/mtk_integration.jl b/test/mtk_integration.jl
index 00e888b2..fbbbe9f6 100644
--- a/test/mtk_integration.jl
+++ b/test/mtk_integration.jl
@@ -7,7 +7,7 @@ D = Differential(t)
 
 function fol_factory(separate = false; name)
     @parameters τ
-    @variables t x(t) f(t) RHS(t)
+    @variables x(t) f(t) RHS(t)
 
     eqs = separate ? [RHS ~ (f - x) / τ,
         D(x) ~ RHS] :
@@ -35,67 +35,73 @@ p = [fol_1.τ => 2.0,
 prob = ODEProblem(sys, u0, (0.0, 10.0), p)
 ds = CoupledODEs(prob)
 
-# parameters
-@test current_parameter(ds, 1) == 2.0
-@test current_parameter(ds, fol_1.τ) == 2.0
-@test current_parameter(ds, 2) == 4.0
-@test current_parameter(ds, fol_2.τ) == 4.0
+@testset "parameter get/set" begin
+    @test current_parameter(ds, 1) == 2.0
+    @test current_parameter(ds, fol_1.τ) == 2.0
+    @test current_parameter(ds, 2) == 4.0
+    @test current_parameter(ds, fol_2.τ) == 4.0
 
-set_parameter!(ds, 1, 3.0)
-@test current_parameter(ds, 1) == 3.0
-@test current_parameter(ds, fol_1.τ) == 3.0
+    set_parameter!(ds, 1, 3.0)
+    @test current_parameter(ds, 1) == 3.0
+    @test current_parameter(ds, fol_1.τ) == 3.0
 
-set_parameter!(ds, fol_1.τ, 2.0)
-@test current_parameter(ds, 1) == 2.0
-@test current_parameter(ds, fol_1.τ) == 2.0
+    set_parameter!(ds, fol_1.τ, 2.0)
+    @test current_parameter(ds, 1) == 2.0
+    @test current_parameter(ds, fol_1.τ) == 2.0
+end
 
 # pure parameter container
-pp = deepcopy(current_parameters(ds))
-set_parameter!(ds, fol_1.τ, 4.0, pp)
-@test current_parameter(ds, fol_1.τ, pp) == 4.0
-
-# states and observed variables
-@test observe_state(ds, 1) == -0.5
-@test observe_state(ds, fol_1.x) == -0.5
-@test observe_state(ds, fol_2.RHS) == -0.375
-
-set_state!(ds, 1.5, 1)
-@test observe_state(ds, 1) == 1.5
-set_state!(ds, -0.5, fol_1.x)
-@test observe_state(ds, 1) == -0.5
-
-# test that derivative dynamical systems also work as execpted
-u1 = current_state(ds)
-pds = ParallelDynamicalSystem(ds, [u1, copy(u1)])
-
-set_parameter!(pds, fol_1.τ, 4.0)
-@test current_parameter(pds, 1) == 4.0
-@test current_parameter(pds, fol_1.τ) == 4.0
-@test observe_state(pds, fol_1.x) == -0.5
-@test observe_state(pds, fol_2.RHS) == -0.375
-
-sds = StroboscopicMap(ds, 1.0)
-set_parameter!(sds, fol_1.τ, 2.0)
-@test current_parameter(sds, 1) == 2.0
-@test current_parameter(sds, fol_1.τ) == 2.0
-@test observe_state(sds, fol_1.x) == -0.5
-@test observe_state(sds, fol_2.RHS) == -0.375
-
-prods = ProjectedDynamicalSystem(ds, [1], [0.0])
-set_parameter!(prods, fol_1.τ, 3.0)
-@test current_parameter(prods, 1) == 3.0
-@test current_parameter(prods, fol_1.τ) == 3.0
-@test observe_state(prods, fol_1.x) == -0.5
-@test observe_state(prods, fol_2.RHS) == -0.375
-
-# notice this evolves the dynamical system!!!
-pmap = PoincareMap(ds, (1, 0.0))
-set_parameter!(pmap, fol_1.τ, 4.0)
-@test current_parameter(pmap, 1) == 4.0
-@test current_parameter(pmap, fol_1.τ) == 4.0
-@test observe_state(pmap, fol_1.x) ≈ 0 atol = 1e-3 rtol = 0
-
-# test with split
+@testset "pure parameter container" begin
+    pp = deepcopy(current_parameters(ds))
+    set_parameter!(ds, fol_1.τ, 4.0, pp)
+    @test current_parameter(ds, fol_1.τ, pp) == 4.0
+end
+
+@testset "states and observed variables" begin
+    @test observe_state(ds, 1) == -0.5
+    @test observe_state(ds, fol_1.x) == -0.5
+    @test observe_state(ds, fol_2.RHS) == -0.375
+
+    set_state!(ds, 1.5, 1)
+    @test observe_state(ds, 1) == 1.5
+    set_state!(ds, -0.5, fol_1.x)
+    @test observe_state(ds, 1) == -0.5
+end
+
+@testset "derivative dynamical systems" begin
+    u1 = current_state(ds)
+    pds = ParallelDynamicalSystem(ds, [u1, copy(u1)])
+
+    set_parameter!(pds, fol_1.τ, 4.0)
+    @test current_parameter(pds, 1) == 4.0
+    @test current_parameter(pds, fol_1.τ) == 4.0
+    @test observe_state(pds, fol_1.x) == -0.5
+    @test observe_state(pds, fol_2.RHS) == -0.375
+
+    sds = StroboscopicMap(ds, 1.0)
+    set_parameter!(sds, fol_1.τ, 2.0)
+    @test current_parameter(sds, 1) == 2.0
+    @test current_parameter(sds, fol_1.τ) == 2.0
+    @test observe_state(sds, fol_1.x) == -0.5
+    @test observe_state(sds, fol_2.RHS) == -0.375
+
+    prods = ProjectedDynamicalSystem(ds, [1], [0.0])
+    set_parameter!(prods, fol_1.τ, 3.0)
+    @test current_parameter(prods, 1) == 3.0
+    @test current_parameter(prods, fol_1.τ) == 3.0
+    @test observe_state(prods, fol_1.x) == -0.5
+    @test observe_state(prods, fol_2.RHS) == -0.375
+
+    # notice this evolves the dynamical system!!!
+    pmap = PoincareMap(ds, (1, 0.0))
+    set_parameter!(pmap, fol_1.τ, 4.0)
+    @test current_parameter(pmap, 1) == 4.0
+    @test current_parameter(pmap, fol_1.τ) == 4.0
+    @test observe_state(pmap, fol_1.x) ≈ 0 atol = 1e-3 rtol = 0
+end
+
+
+@testset "split = true" begin
 sys = structural_simplify(connected; split = true)
 
 u0 = [fol_1.x => -0.5,
@@ -104,32 +110,34 @@ u0 = [fol_1.x => -0.5,
 p = [fol_1.τ => 2.0,
     fol_2.τ => 4.0]
 
-prob = ODEProblem(sys, u0, (0.0, 10.0), p)
-ds = CoupledODEs(prob)
-
-@test current_parameter(ds, fol_1.τ) == 2.0
-set_parameter!(ds, fol_1.τ, 3.0)
-@test current_parameter(ds, fol_1.τ) == 3.0
+    prob = ODEProblem(sys, u0, (0.0, 10.0), p)
+    ds = CoupledODEs(prob)
 
-# test without sys
-function lorenz!(du, u, p, t)
-    du[1] = p[1] * (u[2] - u[1])
-    du[2] = u[1] * (28.0 - u[3]) - u[2]
-    du[3] = u[1] * u[2] - (8 / 3) * u[3]
+    @test current_parameter(ds, fol_1.τ) == 2.0
+    set_parameter!(ds, fol_1.τ, 3.0)
+    @test current_parameter(ds, fol_1.τ) == 3.0
 end
-u0 = [1.0; 0.0; 0.0]
-tspan = (0.0, 100.0)
-p0 = [10.0]
-prob = ODEProblem(lorenz!, u0, tspan, p0)
-ds = CoupledODEs(prob)
 
-@test current_parameter(ds, 1) == 10.0
-set_parameter!(ds, 1, 2.0)
-@test current_parameter(ds, 1) == 2.0
+@testset "no MTK sys" begin
+    function lorenz!(du, u, p, t)
+        du[1] = p[1] * (u[2] - u[1])
+        du[2] = u[1] * (28.0 - u[3]) - u[2]
+        du[3] = u[1] * u[2] - (8 / 3) * u[3]
+    end
+    u0 = [1.0; 0.0; 0.0]
+    tspan = (0.0, 100.0)
+    p0 = [10.0]
+    prob = ODEProblem(lorenz!, u0, tspan, p0)
+    ds = CoupledODEs(prob)
 
-@test observe_state(ds, 1) == 1.0
+    @test current_parameter(ds, 1) == 10.0
+    set_parameter!(ds, 1, 2.0)
+    @test current_parameter(ds, 1) == 2.0
 
-@test_throws ErrorException observe_state(ds, fol_1.f)
+    @test observe_state(ds, 1) == 1.0
+
+    @test_throws ErrorException observe_state(ds, fol_1.f)
+end
 
 # Test that remake works also without anything initial
 

From 0425c0e4ab59f16ea3a9f3003bb3a599f9262f87 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Fri, 14 Mar 2025 13:55:36 +0000
Subject: [PATCH 06/12] add reinit_dae = false

See https://github.com/SciML/ModelingToolkit.jl/issues/3451#
---
 src/core_systems/continuous_time_ode.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/core_systems/continuous_time_ode.jl b/src/core_systems/continuous_time_ode.jl
index e5bb8894..1c21e994 100644
--- a/src/core_systems/continuous_time_ode.jl
+++ b/src/core_systems/continuous_time_ode.jl
@@ -145,7 +145,7 @@ function SciMLBase.reinit!(ds::ContinuousTimeDynamicalSystem, u::AbstractArray =
         p = current_parameters(ds), t0 = initial_time(ds)
     )
     set_parameters!(ds, p)
-    reinit!(ds.integ, u; reset_dt = true, t0)
+    reinit!(ds.integ, u; reinit_dae = false, reset_dt = true, t0)
     return ds
 end
 

From be41a0faad8873e76113a7cf1bfbfae5fef606e8 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Sun, 16 Mar 2025 08:15:15 +0000
Subject: [PATCH 07/12] remove `named_parametersd`

Needs MTK to be defined
---
 src/core/dynamicalsystem_interface.jl | 11 -----------
 1 file changed, 11 deletions(-)

diff --git a/src/core/dynamicalsystem_interface.jl b/src/core/dynamicalsystem_interface.jl
index 22ee8500..9b34df9f 100644
--- a/src/core/dynamicalsystem_interface.jl
+++ b/src/core/dynamicalsystem_interface.jl
@@ -329,17 +329,6 @@ successful_step(ds::DiscreteTimeDynamicalSystem) = all(x -> (isfinite(x) && !isn
 # Generic implementation, most types re-define it as compile-time info
 StateSpaceSets.dimension(ds::DynamicalSystem) = length(current_state(ds))
 
-# TODO: Document it
-function named_current_parameters(ds::DynamicalSystem)
-    mtk = referrenced_sciml_model(ebm)
-    # check if I can use `SymbolicIndexingInterface.parameter_symbols(mtk)` instead.
-    # Does it include initials or not?
-    params_names = Symbol.(ModelingToolkit.parameters(mtk))
-    params_values = current_parameter.(ds, params_names)
-    params = Dict(params_names .=> params_values)
-    return params
-end
-
 ###########################################################################################
 # API - altering status of the system
 ###########################################################################################

From ecda1568630822184201571e9872a93636bc64a8 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Tue, 9 Sep 2025 15:03:02 +0100
Subject: [PATCH 08/12] Named state space sets

---
 CHANGELOG.md                          |  5 +++++
 Project.toml                          |  4 ++--
 docs/src/index.md                     |  1 +
 src/core/dynamicalsystem_interface.jl | 22 +++++++++++++++++++--
 src/core/trajectory.jl                | 28 +++++++++++++++++----------
 5 files changed, 46 insertions(+), 14 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 45153690..80616725 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,3 +1,8 @@
+# v3.15
+
+- New function `named_variables(ds)` for getting the variable names of an MTK-generated dynamical system.
+- `trajectory` function now automatically makes a named state space set if the dynamical system is MTK-generated.
+
 # v3.14.0
 
 - `set_parameter!` and `current_parameter` will now throw an informative error message explicitly naming the parameter if the user tries to get/set a symbolic parameter that does not exist in the MTK-generated dynamical system.
diff --git a/Project.toml b/Project.toml
index 47d6de69..c28b50dc 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "DynamicalSystemsBase"
 uuid = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 repo = "https://github.com/JuliaDynamics/DynamicalSystemsBase.jl.git"
-version = "3.14.0"
+version = "3.15.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -21,7 +21,7 @@ OrdinaryDiffEqTsit5 = "1.1"
 Reexport = "1"
 Roots = "1, 2"
 SciMLBase = "1.19.5, 2"
-StateSpaceSets = "2"
+StateSpaceSets = "2.5"
 Statistics = "1"
 StochasticDiffEq = "6.66.0"
 SymbolicIndexingInterface = "0.3.4"
diff --git a/docs/src/index.md b/docs/src/index.md
index 342a9f67..f382ff27 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -30,6 +30,7 @@ initial_time
 isinplace(::DynamicalSystem)
 successful_step
 referrenced_sciml_model
+named_variables
 ```
 
 ```@docs
diff --git a/src/core/dynamicalsystem_interface.jl b/src/core/dynamicalsystem_interface.jl
index 9b34df9f..abece90a 100644
--- a/src/core/dynamicalsystem_interface.jl
+++ b/src/core/dynamicalsystem_interface.jl
@@ -44,7 +44,7 @@ from an experimental measurement of a dynamical system with an unknown dynamic r
 
 See also the DynamicalSystems.jl tutorial online for examples making dynamical systems.
 
-## Integration with ModelingToolkit.jl
+## Integration with ModelingToolkit.jl (MTK)
 
 Dynamical systems that have been constructed from `DEProblem`s that themselves
 have been constructed from ModelingToolkit.jl keep a reference to the symbolic
@@ -93,6 +93,7 @@ unless when developing new algorithm implementations that use dynamical systems.
 - [`isinplace`](@ref)
 - [`successful_step`](@ref)
 - [`referrenced_sciml_model`](@ref)
+- [`named_variables`](@ref)
 
 ### API - alter status
 
@@ -126,7 +127,7 @@ errormsg(ds) = error("Not yet implemented for dynamical system of type $(nameof(
 
 export current_state, initial_state, current_parameters, current_parameter, initial_parameters, isinplace,
     current_time, initial_time, successful_step, isdeterministic, isdiscretetime, dynamic_rule,
-    reinit!, set_state!, set_parameter!, set_parameters!, step!, observe_state, referrenced_sciml_model
+    reinit!, set_state!, set_parameter!, set_parameters!, step!, observe_state, referrenced_sciml_model, named_variables
 
 ###########################################################################################
 # Symbolic support
@@ -148,10 +149,27 @@ referrenced_sciml_model(::Nothing) = nothing
 
 # return true if there is an actual referrenced system
 has_referrenced_model(prob::SciMLBase.DEProblem) = has_referrenced_model(referrenced_sciml_model(prob))
+has_referrenced_model(prob::DynamicalSystem) = has_referrenced_model(referrenced_sciml_model(prob))
 has_referrenced_model(::Nothing) = false
 has_referrenced_model(::SymbolicIndexingInterface.SymbolCache{Nothing, Nothing, Nothing}) = false
 has_referrenced_model(sys) = true
 
+"""
+    named_variables(ds::DynamicalSystem)
+
+If `ds` is constructed via MTK, return a vector of the variable names (as symbols).
+Otherwise return `nothing`. If `X` is a `StateSpaceSet`, you can always do
+```julia
+X = StateSpaceSet(X; names = named_variables(ds))
+```
+in downstream functions to name a set coming from `ds` (if possible).
+"""
+function named_variables(ds::DynamicalSystem)
+    mtk = referrenced_sciml_model(ds)
+    isnothing(mtk) && return nothing
+    return SymbolicIndexingInterface.getname.(SymbolicIndexingInterface.variable_symbols(mtk))
+end
+
 ###########################################################################################
 # API - obtaining information from the system
 ###########################################################################################
diff --git a/src/core/trajectory.jl b/src/core/trajectory.jl
index acf59c9c..f2e8b464 100644
--- a/src/core/trajectory.jl
+++ b/src/core/trajectory.jl
@@ -5,18 +5,15 @@ export trajectory
 
 Evolve `ds` for a total time of `T` and return its trajectory `X`, sampled at equal time
 intervals, and corresponding time vector. `X` is a `StateSpaceSet`.
-Optionally provide a starting state `u0` which is `current_state(ds)` by default.
-
 The returned time vector is `t = (t0+Ttr):Δt:(t0+Ttr+T)`.
 
-If time evolution diverged, or in general failed, before `T`,
+Optionally provide a starting state `u0` which is `current_state(ds)` by default.
+
+If time evolution diverged or in general failed before `T`,
 the remaining of the trajectory is set to the last valid point.
 
-`trajectory` is a very simple function provided for convenience.
-For continuous time systems, it doesn't play well with callbacks,
-use `DifferentialEquations.solve` if you want a trajectory/timeseries
-that works with callbacks, or in general you want more flexibility in the generated
-trajectory (but remember to convert the output of `solve` to a `StateSpaceSet`).
+The dimensions of `X` are automatically named if `ds` referrences an MTK model,
+see [`DynamicalSystem`](@ref).
 
 ## Keyword arguments
 
@@ -32,6 +29,14 @@ trajectory (but remember to convert the output of `solve` to a `StateSpaceSet`).
   Defaults to `1:dimension(ds)` (all dynamic variables).
   Note: if you mix integer and symbolic indexing be sure to initialize the array
   as `Any` so that integers `1, 2, ...` are not converted to symbolic expressions.
+
+## Description
+
+`trajectory` is a very simple function provided for convenience.
+For continuous time systems, it doesn't play well with callbacks,
+use `DifferentialEquations.solve` if you want a trajectory
+that works with callbacks, or in general you want more flexibility in the generated
+trajectory (but remember to convert the output of `solve` to a `StateSpaceSet`).
 """
 function trajectory(ds::DynamicalSystem, T, u0 = current_state(ds);
         save_idxs = nothing, t0 = initial_time(ds), kwargs...
@@ -39,10 +44,13 @@ function trajectory(ds::DynamicalSystem, T, u0 = current_state(ds);
     accessor = svector_access(save_idxs)
     reinit!(ds, u0; t0)
     if isdiscretetime(ds)
-        trajectory_discrete(ds, T; accessor, kwargs...)
+        X, tvec = trajectory_discrete(ds, T; accessor, kwargs...)
     else
-        trajectory_continuous(ds, T; accessor, kwargs...)
+        X, tvec = trajectory_continuous(ds, T; accessor, kwargs...)
     end
+    # name automatically if possible
+    X = StateSpaceSet(X; names = named_variables(ds))
+    return X, tvec
 end
 
 function trajectory_discrete(ds, T;

From 76cb085ca3a29007ef0bf1229abb2c9104d54fec Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Tue, 9 Sep 2025 15:16:53 +0100
Subject: [PATCH 09/12] add tests for trajectory and named variables

---
 src/core/trajectory.jl  | 13 +++++++------
 test/mtk_integration.jl | 16 +++++++++++-----
 2 files changed, 18 insertions(+), 11 deletions(-)

diff --git a/src/core/trajectory.jl b/src/core/trajectory.jl
index f2e8b464..143c72db 100644
--- a/src/core/trajectory.jl
+++ b/src/core/trajectory.jl
@@ -12,8 +12,8 @@ Optionally provide a starting state `u0` which is `current_state(ds)` by default
 If time evolution diverged or in general failed before `T`,
 the remaining of the trajectory is set to the last valid point.
 
-The dimensions of `X` are automatically named if `ds` referrences an MTK model,
-see [`DynamicalSystem`](@ref).
+The dimensions of `X` are automatically named if `ds` referrences an MTK model
+and if `save_idxs` remains at its default value.
 
 ## Keyword arguments
 
@@ -24,9 +24,8 @@ see [`DynamicalSystem`](@ref).
 * `t0 = initial_time(ds)`: Starting time.
 * `container = SVector`: Type of vector that will represent the state space points
   that will be included in the `StateSpaceSet` output. See `StateSpaceSet` for valid options.
-* `save_idxs::AbstractVector`: Which variables to output in `X`. It can be
-  any type of index that can be given to [`observe_state`](@ref).
-  Defaults to `1:dimension(ds)` (all dynamic variables).
+* `save_idxs`: Which variables to output in `X`. By default it is `nothing` (all variables).
+  It can be a vector of any type of index that can be given to [`observe_state`](@ref).
   Note: if you mix integer and symbolic indexing be sure to initialize the array
   as `Any` so that integers `1, 2, ...` are not converted to symbolic expressions.
 
@@ -49,7 +48,9 @@ function trajectory(ds::DynamicalSystem, T, u0 = current_state(ds);
         X, tvec = trajectory_continuous(ds, T; accessor, kwargs...)
     end
     # name automatically if possible
-    X = StateSpaceSet(X; names = named_variables(ds))
+    if isnothing(save_idxs)
+        X = StateSpaceSet(X; names = named_variables(ds))
+    end
     return X, tvec
 end
 
diff --git a/test/mtk_integration.jl b/test/mtk_integration.jl
index fbbbe9f6..2617b345 100644
--- a/test/mtk_integration.jl
+++ b/test/mtk_integration.jl
@@ -68,6 +68,12 @@ end
     @test observe_state(ds, 1) == -0.5
 end
 
+@testset "trajectory naming" begin
+    @test named_variables(ds) == [:fol_1₊x, :fol_2₊x]
+    X, tvec = trajectory(ds, 10)
+    @test X.names == [:fol_1₊x, :fol_2₊x]
+end
+
 @testset "derivative dynamical systems" begin
     u1 = current_state(ds)
     pds = ParallelDynamicalSystem(ds, [u1, copy(u1)])
@@ -102,13 +108,13 @@ end
 
 
 @testset "split = true" begin
-sys = structural_simplify(connected; split = true)
+    sys = structural_simplify(connected; split = true)
 
-u0 = [fol_1.x => -0.5,
-    fol_2.x => 1.0]
+    u0 = [fol_1.x => -0.5,
+        fol_2.x => 1.0]
 
-p = [fol_1.τ => 2.0,
-    fol_2.τ => 4.0]
+    p = [fol_1.τ => 2.0,
+        fol_2.τ => 4.0]
 
     prob = ODEProblem(sys, u0, (0.0, 10.0), p)
     ds = CoupledODEs(prob)

From 5090ce1e813f56d4543df219b946b43b8e76a0e8 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Wed, 10 Sep 2025 10:09:25 +0100
Subject: [PATCH 10/12] update tests to MTK v10

---
 test/Project.toml       |  2 +-
 test/mtk_integration.jl | 77 ++++++++++++++++-------------------------
 2 files changed, 31 insertions(+), 48 deletions(-)

diff --git a/test/Project.toml b/test/Project.toml
index 1a3a01b4..08533fad 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -10,4 +10,4 @@ StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
-ModelingToolkit = "9"
+ModelingToolkit = "10"
diff --git a/test/mtk_integration.jl b/test/mtk_integration.jl
index 2617b345..e219172a 100644
--- a/test/mtk_integration.jl
+++ b/test/mtk_integration.jl
@@ -22,9 +22,9 @@ end
 connections = [fol_1.f ~ 1.5,
     fol_2.f ~ fol_1.x]
 
-connected = compose(ODESystem(connections, t; name = :connected), fol_1, fol_2)
+connected = compose(System(connections, t; name = :connected), fol_1, fol_2)
 
-sys = structural_simplify(connected; split = false)
+sys = mtkcompile(connected; split = false)
 
 u0 = [fol_1.x => -0.5,
     fol_2.x => 1.0]
@@ -32,46 +32,34 @@ u0 = [fol_1.x => -0.5,
 p = [fol_1.τ => 2.0,
     fol_2.τ => 4.0]
 
-prob = ODEProblem(sys, u0, (0.0, 10.0), p)
+initials = vcat(u0, p)
+
+prob = ODEProblem(sys, initials, (0.0, 10.0))
 ds = CoupledODEs(prob)
 
 @testset "parameter get/set" begin
-    @test current_parameter(ds, 1) == 2.0
     @test current_parameter(ds, fol_1.τ) == 2.0
-    @test current_parameter(ds, 2) == 4.0
     @test current_parameter(ds, fol_2.τ) == 4.0
-
-    set_parameter!(ds, 1, 3.0)
-    @test current_parameter(ds, 1) == 3.0
-    @test current_parameter(ds, fol_1.τ) == 3.0
-
-    set_parameter!(ds, fol_1.τ, 2.0)
-    @test current_parameter(ds, 1) == 2.0
-    @test current_parameter(ds, fol_1.τ) == 2.0
+    set_parameter!(ds, fol_1.τ, 4.0)
+    @test current_parameter(ds, fol_1.τ) == 4.0
 end
 
 # pure parameter container
 @testset "pure parameter container" begin
     pp = deepcopy(current_parameters(ds))
-    set_parameter!(ds, fol_1.τ, 4.0, pp)
-    @test current_parameter(ds, fol_1.τ, pp) == 4.0
+    set_parameter!(ds, fol_1.τ, 2.0, pp)
+    @test current_parameter(ds, fol_1.τ, pp) == 2.0
 end
 
 @testset "states and observed variables" begin
-    @test observe_state(ds, 1) == -0.5
     @test observe_state(ds, fol_1.x) == -0.5
     @test observe_state(ds, fol_2.RHS) == -0.375
 
-    set_state!(ds, 1.5, 1)
-    @test observe_state(ds, 1) == 1.5
+    x0 = observe_state(ds, 1)
+    set_state!(ds, x0, 1)
+    @test observe_state(ds, 1) == x0
     set_state!(ds, -0.5, fol_1.x)
-    @test observe_state(ds, 1) == -0.5
-end
-
-@testset "trajectory naming" begin
-    @test named_variables(ds) == [:fol_1₊x, :fol_2₊x]
-    X, tvec = trajectory(ds, 10)
-    @test X.names == [:fol_1₊x, :fol_2₊x]
+    @test observe_state(ds, fol_1.x) == -0.5
 end
 
 @testset "derivative dynamical systems" begin
@@ -79,46 +67,41 @@ end
     pds = ParallelDynamicalSystem(ds, [u1, copy(u1)])
 
     set_parameter!(pds, fol_1.τ, 4.0)
-    @test current_parameter(pds, 1) == 4.0
     @test current_parameter(pds, fol_1.τ) == 4.0
     @test observe_state(pds, fol_1.x) == -0.5
     @test observe_state(pds, fol_2.RHS) == -0.375
 
     sds = StroboscopicMap(ds, 1.0)
     set_parameter!(sds, fol_1.τ, 2.0)
-    @test current_parameter(sds, 1) == 2.0
     @test current_parameter(sds, fol_1.τ) == 2.0
     @test observe_state(sds, fol_1.x) == -0.5
     @test observe_state(sds, fol_2.RHS) == -0.375
 
     prods = ProjectedDynamicalSystem(ds, [1], [0.0])
     set_parameter!(prods, fol_1.τ, 3.0)
-    @test current_parameter(prods, 1) == 3.0
     @test current_parameter(prods, fol_1.τ) == 3.0
     @test observe_state(prods, fol_1.x) == -0.5
     @test observe_state(prods, fol_2.RHS) == -0.375
 
     # notice this evolves the dynamical system!!!
-    pmap = PoincareMap(ds, (1, 0.0))
+    pmap = PoincareMap(ds, (2, 0.0))
     set_parameter!(pmap, fol_1.τ, 4.0)
-    @test current_parameter(pmap, 1) == 4.0
     @test current_parameter(pmap, fol_1.τ) == 4.0
     @test observe_state(pmap, fol_1.x) ≈ 0 atol = 1e-3 rtol = 0
 end
 
+@testset "trajectory naming" begin
+    vars = named_variables(ds)
+    @test length(vars) == 2
+    @test sort!(string.(vars)) == ["fol_1₊x", "fol_2₊x"]
+    X, tvec = trajectory(ds, 10)
+    @test X.names == vars
+end
 
 @testset "split = true" begin
-    sys = structural_simplify(connected; split = true)
-
-    u0 = [fol_1.x => -0.5,
-        fol_2.x => 1.0]
-
-    p = [fol_1.τ => 2.0,
-        fol_2.τ => 4.0]
-
-    prob = ODEProblem(sys, u0, (0.0, 10.0), p)
+    sys = mtkcompile(connected; split = true)
+    prob = ODEProblem(sys, initials, (0.0, 10.0))
     ds = CoupledODEs(prob)
-
     @test current_parameter(ds, fol_1.τ) == 2.0
     set_parameter!(ds, fol_1.τ, 3.0)
     @test current_parameter(ds, fol_1.τ) == 3.0
@@ -142,7 +125,7 @@ end
 
     @test observe_state(ds, 1) == 1.0
 
-    @test_throws ErrorException observe_state(ds, fol_1.f)
+    @test_throws ArgumentError observe_state(ds, fol_1.f)
 end
 
 # Test that remake works also without anything initial
@@ -169,11 +152,11 @@ D = Differential(t)
     end
 end
 
-@mtkbuild roessler_model = Roessler()
+@mtkcompile roessler_model = Roessler()
 
 @testset "type: $(iip)" for iip in (true, false)
     if iip
-        prob = ODEProblem(roessler_model)
+        prob = ODEProblem(roessler_model, nothing, (0.0, Inf))
     else
         prob = ODEProblem{false}(roessler_model, nothing, (0.0, Inf); u0_constructor = x->SVector(x...))
     end
@@ -215,7 +198,7 @@ end
 end
 
 @testset "informative errors" begin
-    prob = ODEProblem(roessler_model)
+    prob = prob = ODEProblem(roessler_model, nothing, (0.0, Inf))
     ds = CoupledODEs(prob)
     @parameters XOXO = 0.5
     @test_throws "XOXO" current_parameter(ds, :XOXO)
@@ -237,10 +220,10 @@ Differential(t)(DS) ~ η2 - η3*DS - abs(DT - DS)*DS,
 η1 ~ η1_0 + r_η*t, # this symbolic variable has its own equation!
 ]
 
-sys = ODESystem(eqs, t; name = :stommel)
-sys = structural_simplify(sys; split = false)
+sys = System(eqs, t; name = :stommel)
+sys = mtkcompile(sys; split = false)
 
-prob = ODEProblem(sys)
+prob = ODEProblem(sys, nothing, (0.0, Inf))
 ds = CoupledODEs(prob)
 
 X, tvec = trajectory(ds, 10.0; Δt = 0.1, save_idxs = Any[1, 2, η1])

From 85119f732582d5f10dd159984b12b962f9e6c906 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Wed, 10 Sep 2025 10:20:41 +0100
Subject: [PATCH 11/12] some typos

---
 test/mtk_integration.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/test/mtk_integration.jl b/test/mtk_integration.jl
index e219172a..a726dbf7 100644
--- a/test/mtk_integration.jl
+++ b/test/mtk_integration.jl
@@ -13,7 +13,7 @@ function fol_factory(separate = false; name)
         D(x) ~ RHS] :
           D(x) ~ (f - x) / τ
 
-    ODESystem(eqs, t; name)
+    System(eqs, t; name)
 end
 
 @named fol_1 = fol_factory()
@@ -112,6 +112,7 @@ end
         du[1] = p[1] * (u[2] - u[1])
         du[2] = u[1] * (28.0 - u[3]) - u[2]
         du[3] = u[1] * u[2] - (8 / 3) * u[3]
+        return nothing
     end
     u0 = [1.0; 0.0; 0.0]
     tspan = (0.0, 100.0)

From 2b0af76cbf7cf31ac3f5e6a40758b6dbf96943e3 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Wed, 10 Sep 2025 13:07:17 +0100
Subject: [PATCH 12/12] fix jacobian

---
 test/jacobian.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/jacobian.jl b/test/jacobian.jl
index 56503c7c..bf33f9e3 100644
--- a/test/jacobian.jl
+++ b/test/jacobian.jl
@@ -30,7 +30,7 @@ end
 
 @testset "MTK Jacobian" begin
     using ModelingToolkit
-    using ModelingToolkit: Num, GeneratedFunctionWrapper
+    using ModelingToolkit: GeneratedFunctionWrapper
     using DynamicalSystemsBase: SciMLBase
     @independent_variables t
     @variables u(t)[1:2]
@@ -38,15 +38,15 @@ end
 
     eqs = [D(u[1]) ~ 3.0 * u[1],
            D(u[2]) ~ -3.0 * u[2]]
-    @named sys = ODESystem(eqs, t)
-    sys = structural_simplify(sys)
+    @named sys = System(eqs, t)
+    sys = mtkcompile(sys)
 
     prob = ODEProblem(sys, [1.0, 1.0], (0.0, 1.0); jac=true)
     ds = CoupledODEs(prob)
 
     jac = jacobian(ds)
     @test jac isa GeneratedFunctionWrapper
-    @test jac([1.0, 1.0], [], 0.0) == [3 0;0 -3]
+    @test jac([1.0, 1.0], [], 0.0) == [-3 0;0 3]
 
     @testset "CoupledSDEs" begin
         # just to check if MTK @brownian does not give any problems
@@ -61,6 +61,6 @@ end
 
         jac = jacobian(sde)
         @test jac isa GeneratedFunctionWrapper
-        @test jac([1.0, 1.0], [], 0.0) == [3 0; 0 -3]
+        @test jac([1.0, 1.0], [], 0.0) == [-3 0; 0 3]
     end
 end