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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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/19] 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 e2fde7c9d166769f7267e84ec8a5b234d93b5b3c Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Tue, 1 Apr 2025 11:54:40 +0100
Subject: [PATCH 08/19] allow kw propagation in reinit

---
 src/core_systems/continuous_time_ode.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/core_systems/continuous_time_ode.jl b/src/core_systems/continuous_time_ode.jl
index 1c21e994..48aa58a0 100644
--- a/src/core_systems/continuous_time_ode.jl
+++ b/src/core_systems/continuous_time_ode.jl
@@ -142,10 +142,10 @@ set_state!(ds::CoupledODEs, u::AbstractArray) = (set_state!(ds.integ, u); ds)
 SciMLBase.step!(ds::CoupledODEs, args...) = (step!(ds.integ, args...); ds)
 
 function SciMLBase.reinit!(ds::ContinuousTimeDynamicalSystem, u::AbstractArray = initial_state(ds);
-        p = current_parameters(ds), t0 = initial_time(ds)
+        p = current_parameters(ds), t0 = initial_time(ds), kw...
     )
     set_parameters!(ds, p)
-    reinit!(ds.integ, u; reinit_dae = false, reset_dt = true, t0)
+    reinit!(ds.integ, u; reinit_dae = false, reset_dt = true, t0, kw...)
     return ds
 end
 

From dd98da65f59e5cfcb44986127969e6b6947b4393 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Tue, 1 Apr 2025 12:09:41 +0100
Subject: [PATCH 09/19] fix reinit for SDE

---
 src/core_systems/continuous_time_ode.jl |  2 +-
 src/core_systems/continuous_time_sde.jl | 12 ++++++++++--
 2 files changed, 11 insertions(+), 3 deletions(-)

diff --git a/src/core_systems/continuous_time_ode.jl b/src/core_systems/continuous_time_ode.jl
index 48aa58a0..b2b6b480 100644
--- a/src/core_systems/continuous_time_ode.jl
+++ b/src/core_systems/continuous_time_ode.jl
@@ -141,7 +141,7 @@ set_state!(ds::CoupledODEs, u::AbstractArray) = (set_state!(ds.integ, u); ds)
 # so that `ds` is printed
 SciMLBase.step!(ds::CoupledODEs, args...) = (step!(ds.integ, args...); ds)
 
-function SciMLBase.reinit!(ds::ContinuousTimeDynamicalSystem, u::AbstractArray = initial_state(ds);
+function SciMLBase.reinit!(ds::ContinuousTimeDynamicalSystem, u = initial_state(ds);
         p = current_parameters(ds), t0 = initial_time(ds), kw...
     )
     set_parameters!(ds, p)
diff --git a/src/core_systems/continuous_time_sde.jl b/src/core_systems/continuous_time_sde.jl
index bcd9a789..ce6155c8 100644
--- a/src/core_systems/continuous_time_sde.jl
+++ b/src/core_systems/continuous_time_sde.jl
@@ -13,7 +13,7 @@ where
 ``\\mathbf{u}(t)`` is the state vector at time ``t``, ``\\mathbf{f}`` describes the
 deterministic dynamics, and the noise term
 ``\\mathbf{g}(\\mathbf{u}, p, t) \\text{d}\\mathcal{N}_t`` describes
-the stochastic forcing in terms of a noise function (or *diffusion function*) 
+the stochastic forcing in terms of a noise function (or *diffusion function*)
 ``\\mathbf{g}`` and a noise process ``\\mathcal{N}_t``. The parameters of the functions
 ``\\mathcal{f}`` and ``\\mathcal{g}`` are contained in the vector ``p``.
 
@@ -55,7 +55,7 @@ Commonly, this is a matrix or sparse matrix. If this is not given, it
 defaults to `nothing`, which means the `g` should be interpreted as being diagonal.
 
 The noise process can be specified via `noise_process`. It defaults to a standard Wiener
-process (Gaussian white noise). 
+process (Gaussian white noise).
 For details on defining noise processes, see the docs of [DiffEqNoiseProcess.jl
 ](https://docs.sciml.ai/DiffEqNoiseProcess/stable/). A complete list of the pre-defined
 processes can be found [here](https://docs.sciml.ai/DiffEqNoiseProcess/stable/noise_processes/).
@@ -111,4 +111,12 @@ struct CoupledSDEs{IIP,D,I,P} <: ContinuousTimeDynamicalSystem
     p0::P
     diffeq # isn't parameterized because it is only used for display
     noise_type::NamedTuple
+end
+
+function SciMLBase.reinit!(ds::CoupledSDEs, u = initial_state(ds);
+    p = current_parameters(ds), t0 = initial_time(ds), kw...
+  )
+  set_parameters!(ds, p)
+  reinit!(ds.integ, u; t0, kw...)
+  return ds
 end
\ No newline at end of file

From 1c44a77c3e52c996634fdf834efde76463738d5c Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Thu, 3 Apr 2025 15:44:39 +0100
Subject: [PATCH 10/19] simplify ODE test dependency

---
 src/core_systems/continuous_time_ode.jl |  2 +-
 test/Project.toml                       |  3 ++-
 test/continuous.jl                      | 18 ++++--------------
 3 files changed, 7 insertions(+), 16 deletions(-)

diff --git a/src/core_systems/continuous_time_ode.jl b/src/core_systems/continuous_time_ode.jl
index b2b6b480..48aa58a0 100644
--- a/src/core_systems/continuous_time_ode.jl
+++ b/src/core_systems/continuous_time_ode.jl
@@ -141,7 +141,7 @@ set_state!(ds::CoupledODEs, u::AbstractArray) = (set_state!(ds.integ, u); ds)
 # so that `ds` is printed
 SciMLBase.step!(ds::CoupledODEs, args...) = (step!(ds.integ, args...); ds)
 
-function SciMLBase.reinit!(ds::ContinuousTimeDynamicalSystem, u = initial_state(ds);
+function SciMLBase.reinit!(ds::ContinuousTimeDynamicalSystem, u::AbstractArray = initial_state(ds);
         p = current_parameters(ds), t0 = initial_time(ds), kw...
     )
     set_parameters!(ds, p)
diff --git a/test/Project.toml b/test/Project.toml
index 90f1118d..e6332ab3 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -3,7 +3,8 @@ DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/continuous.jl b/test/continuous.jl
index d66dfd38..c19b27a2 100644
--- a/test/continuous.jl
+++ b/test/continuous.jl
@@ -1,6 +1,7 @@
 using DynamicalSystemsBase, Test
 
-using OrdinaryDiffEq: Vern9, ODEProblem, Rodas5, Tsit5
+using OrdinaryDiffEqTsit5: ODEProblem, Tsit5
+using OrdinaryDiffEqRosenbrock: Rodas5
 
 include("test_system_function.jl")
 
@@ -25,14 +26,9 @@ p0 = [10, 28, 8/3]
 
 lorenz_oop = CoupledODEs(lorenz_rule, u0, p0)
 lorenz_iip = CoupledODEs(ODEProblem(lorenz_rule_iip, copy(u0), (0.0, Inf), p0))
-lorenz_vern = CoupledODEs(lorenz_rule, u0, p0;
-    diffeq = (alg = Vern9(), abstol = 1e-9, reltol = 1e-9)
-)
 
-for (ds, iip) in zip((lorenz_oop, lorenz_iip, lorenz_vern), (false, true, false))
-
-    name = (ds === lorenz_vern) ? "lorvern" : "lorenz"
-    @testset "$(name) IIP=$(iip)" begin
+for (ds, iip) in zip((lorenz_oop, lorenz_iip), (false, true))
+    @testset "IIP=$(iip)" begin
         @test dynamic_rule(ds) == (iip ? lorenz_rule_iip : lorenz_rule)
         test_dynamical_system(ds, u0, p0; idt = false, iip=iip)
     end
@@ -42,12 +38,6 @@ end
     lorenz_oop = CoupledODEs(lorenz_rule, u0, p0)
     @test lorenz_oop.integ.alg isa Tsit5
 
-    lorenz_vern = CoupledODEs(lorenz_rule, u0, p0;
-        diffeq = (alg = Vern9(), verbose = false, abstol = 1e-9, reltol = 1e-9)
-    )
-    @test lorenz_vern.integ.alg isa Vern9
-    @test lorenz_vern.integ.opts.verbose == false
-
     # also test ODEproblem creation
     prob = lorenz_vern.integ.sol.prob
 

From 8eadc1c0023963efd298fc8368b669c448e05204 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Thu, 3 Apr 2025 15:45:12 +0100
Subject: [PATCH 11/19] continuous tests pass

---
 test/continuous.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/continuous.jl b/test/continuous.jl
index c19b27a2..6a99c4a8 100644
--- a/test/continuous.jl
+++ b/test/continuous.jl
@@ -39,7 +39,7 @@ end
     @test lorenz_oop.integ.alg isa Tsit5
 
     # also test ODEproblem creation
-    prob = lorenz_vern.integ.sol.prob
+    prob = lorenz_oop.integ.sol.prob
 
     ds = CoupledODEs(prob, (alg=Rodas5(autodiff=false), abstol=0.0, reltol=1e-6, verbose=false))
 

From d6d908b1fcbb1c3e370e6cb298dc2c2b21fb47c3 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Thu, 3 Apr 2025 15:49:44 +0100
Subject: [PATCH 12/19] all tests pass?

---
 test/Project.toml    | 1 +
 test/stochastic.jl   | 2 +-
 test/stroboscopic.jl | 2 +-
 3 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/test/Project.toml b/test/Project.toml
index e6332ab3..8a40e594 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -5,6 +5,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/stochastic.jl b/test/stochastic.jl
index 8d2bebd0..e45039ec 100644
--- a/test/stochastic.jl
+++ b/test/stochastic.jl
@@ -1,5 +1,5 @@
 using DynamicalSystemsBase, Test
-using OrdinaryDiffEq: Tsit5
+using OrdinaryDiffEqTsit5: Tsit5
 using StochasticDiffEq: SDEProblem, SRA, SOSRA, LambaEM, CorrelatedWienerProcess, EM
 
 StochasticSystemsBase = Base.get_extension(DynamicalSystemsBase, :StochasticSystemsBase)
diff --git a/test/stroboscopic.jl b/test/stroboscopic.jl
index 06856cf9..e1245678 100644
--- a/test/stroboscopic.jl
+++ b/test/stroboscopic.jl
@@ -1,6 +1,6 @@
 using DynamicalSystemsBase, Test
 
-using OrdinaryDiffEq: Vern9
+using OrdinaryDiffEqVerner: Vern9
 
 include("test_system_function.jl")
 

From 19a07f621842fd506fa9dbbb4849f449ee8c7ca9 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Thu, 3 Apr 2025 15:52:15 +0100
Subject: [PATCH 13/19] don't forget array for SDEs

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

diff --git a/src/core_systems/continuous_time_sde.jl b/src/core_systems/continuous_time_sde.jl
index ce6155c8..32b83775 100644
--- a/src/core_systems/continuous_time_sde.jl
+++ b/src/core_systems/continuous_time_sde.jl
@@ -113,7 +113,7 @@ struct CoupledSDEs{IIP,D,I,P} <: ContinuousTimeDynamicalSystem
     noise_type::NamedTuple
 end
 
-function SciMLBase.reinit!(ds::CoupledSDEs, u = initial_state(ds);
+function SciMLBase.reinit!(ds::CoupledSDEs, u::AbstractArray = initial_state(ds);
     p = current_parameters(ds), t0 = initial_time(ds), kw...
   )
   set_parameters!(ds, p)

From bd2b04afeb0b005e1a0d77bb05ef2e56be876329 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Thu, 3 Apr 2025 15:56:18 +0100
Subject: [PATCH 14/19] abvoid method overwrites in tests

---
 test/stochastic_characterisation.jl | 14 +++++++-------
 1 file changed, 7 insertions(+), 7 deletions(-)

diff --git a/test/stochastic_characterisation.jl b/test/stochastic_characterisation.jl
index 7dcb1cf4..870882ba 100644
--- a/test/stochastic_characterisation.jl
+++ b/test/stochastic_characterisation.jl
@@ -25,7 +25,7 @@ f!(du, u, p, t) = du .= 1.01u
     @test CoupledSDEs(sciml_prob(corr)).noise_type == types
 
     # invertible correlated Wiener noise
-    g!(du, u, p, t) = (du .= [1 0.3; 0.3 1]; return nothing)
+    g! = (du, u, p, t) -> (du .= [1 0.3; 0.3 1]; return nothing)
     corr_alt = CoupledSDEs(f!, zeros(2); g = g!, noise_prototype = zeros(2, 2))
     types = corr_alt.noise_type
     @test values(types) == (true, true, true, true)
@@ -38,14 +38,14 @@ f!(du, u, p, t) = du .= 1.01u
     @test CoupledSDEs(sciml_prob(corr_noninv)).noise_type == types
 
     # non-invertible correlated Wiener noise
-    g!(du, u, p, t) = (du .= [1 0.3 1; 0.3 1 1]; return nothing)
+    g! = (du, u, p, t) -> (du .= [1 0.3 1; 0.3 1 1]; return nothing)
     corr_alt = CoupledSDEs(f!, zeros(2); g = g!, noise_prototype = zeros(2, 3))
     types = corr_alt.noise_type
     @test values(types) == (true, true, true, false)
     @test CoupledSDEs(sciml_prob(corr_alt)).noise_type == types
 
     # non-autonomous Wiener noise
-    g!(du, u, p, t) = (du .= 1 / (1 + t); return nothing)
+    g! = (du, u, p, t) -> (du .= 1 / (1 + t); return nothing)
     addit_non_autom = CoupledSDEs(f!, zeros(2); g = g!)
     types = addit_non_autom.noise_type
     @test values(types) == (true, false, true, false)
@@ -54,7 +54,7 @@ end
 
 @testset "Multiplicative noise" begin
     # multiplicative linear Wiener noise
-    g!(du, u, p, t) = (du .= u; return nothing)
+    g! = (du, u, p, t) -> (du .= u; return nothing)
     linear_multipli = CoupledSDEs(f!, rand(2) ./ 10; g = g!)
     types = linear_multipli.noise_type
     @test values(types) == (false, true, true, false)
@@ -67,21 +67,21 @@ end
     @test CoupledSDEs(sciml_prob(lin_multipli_alt)).noise_type == types
 
     # multiplicative nonlinear Wiener noise
-    g!(du, u, p, t) = (du .= u .^ 2; return nothing)
+    g! = (du, u, p, t) -> (du .= u .^ 2; return nothing)
     nonlinear_multiplicative = CoupledSDEs(f!, rand(2); g = g!)
     types = nonlinear_multiplicative.noise_type
     @test values(types) == (false, true, false, false)
     @test CoupledSDEs(sciml_prob(nonlinear_multiplicative)).noise_type == types
 
     # non-autonomous linear multiplicative Wiener noise
-    g!(du, u, p, t) = (du .= u ./ (1 + t); return nothing)
+    g! = (du, u, p, t) -> (du .= u ./ (1 + t); return nothing)
     addit_non_autom = CoupledSDEs(f!, zeros(2); g = g!)
     types = addit_non_autom.noise_type
     @test values(types) == (false, false, true, false)
     @test CoupledSDEs(sciml_prob(addit_non_autom)).noise_type == types
 
     # non-autonomous nonlinear multiplicative Wiener noise
-    g!(du, u, p, t) = (du .= u .^ 2 ./ (1 + t); return nothing)
+    g! = (du, u, p, t) -> (du .= u .^ 2 ./ (1 + t); return nothing)
     addit_non_autom = CoupledSDEs(f!, zeros(2); g = g!)
     types = addit_non_autom.noise_type
     @test values(types) == (false, false, false, false)

From cb158c98ee0af581fbfbb12e73d3f414c1321f7f Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Thu, 3 Apr 2025 15:56:33 +0100
Subject: [PATCH 15/19] remove RODAS usage

---
 CHANGELOG.md       |  1 +
 test/Project.toml  |  1 -
 test/continuous.jl | 11 +++--------
 3 files changed, 4 insertions(+), 9 deletions(-)

diff --git a/CHANGELOG.md b/CHANGELOG.md
index 45153690..82c064ca 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -1,6 +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.
+- Fixed a plethora of bugs and test failures associated with updates in the SciML ecosystem.
 
 # v3.13.0
 
diff --git a/test/Project.toml b/test/Project.toml
index 8a40e594..1a3a01b4 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -3,7 +3,6 @@ DiffEqNoiseProcess = "77a26b50-5914-5dd7-bc55-306e6241c503"
 DynamicalSystemsBase = "6e36e845-645a-534a-86f2-f5d4aa5a06b4"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
-OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
diff --git a/test/continuous.jl b/test/continuous.jl
index 6a99c4a8..507a6671 100644
--- a/test/continuous.jl
+++ b/test/continuous.jl
@@ -1,7 +1,7 @@
 using DynamicalSystemsBase, Test
 
 using OrdinaryDiffEqTsit5: ODEProblem, Tsit5
-using OrdinaryDiffEqRosenbrock: Rodas5
+using OrdinaryDiffEqVerner: Vern9
 
 include("test_system_function.jl")
 
@@ -38,14 +38,9 @@ end
     lorenz_oop = CoupledODEs(lorenz_rule, u0, p0)
     @test lorenz_oop.integ.alg isa Tsit5
 
-    # also test ODEproblem creation
     prob = lorenz_oop.integ.sol.prob
-
-    ds = CoupledODEs(prob, (alg=Rodas5(autodiff=false), abstol=0.0, reltol=1e-6, verbose=false))
-
-    @test ds.integ.alg isa Rodas5
+    ds = CoupledODEs(prob, (alg=Vern9(), abstol=0.0, reltol=1e-6, verbose=false))
+    @test ds.integ.alg isa Vern9
     @test ds.integ.opts.verbose == false
 
-    @test_throws ArgumentError CoupledODEs(prob; diffeq = (alg=Rodas5(autodiff=false), ))
-
 end

From 445bac19f161cd92f9917645c7fa8d8dcd5a7760 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Thu, 3 Apr 2025 16:28:16 +0100
Subject: [PATCH 16/19] fix jacobian (?)

---
 src/core_systems/jacobian.jl | 13 ++++++++++++-
 test/jacobian.jl             |  5 ++---
 2 files changed, 14 insertions(+), 4 deletions(-)

diff --git a/src/core_systems/jacobian.jl b/src/core_systems/jacobian.jl
index 41958750..c0198106 100644
--- a/src/core_systems/jacobian.jl
+++ b/src/core_systems/jacobian.jl
@@ -6,7 +6,7 @@ import ForwardDiff
     jacobian(ds::CoreDynamicalSystem)
 
 Construct the Jacobian rule for the dynamical system `ds`.
-If the system already has a Jacobian rule constructed via ModelingToolkit it returns this,
+If the system already has a Jacobian rule constructed via ModelingToolkit.jl it returns this,
 otherwise it constructs the Jacobian rule with automatic differentiation using module
 [`ForwardDiff`](https://github.com/JuliaDiff/ForwardDiff.jl).
 
@@ -21,6 +21,17 @@ at the state `u`, parameters `p` and time `t` and save the result in `J0`.
 """
 function jacobian(ds::CoreDynamicalSystem{IIP}) where {IIP}
     if ds isa ContinuousTimeDynamicalSystem
+        # TODO: This is the correct API way to obtain the Jacobian,
+        # however it relies on MTK dependency, so we can't do it.
+        # if has_referrenced_model(ds)
+        #     model = referrenced_sciml_model(ds)
+        #     Joop, Jiip = generate_jacobian(model; expression = Val{false})
+        #     if IIP
+        #         jac = Jiip
+        #     else
+        #         jac = Joop
+        #     end
+        # end
         prob = referrenced_sciml_prob(ds)
         if prob.f isa SciMLBase.AbstractDiffEqFunction && !isnothing(prob.f.jac)
             jac = prob.f.jac
diff --git a/test/jacobian.jl b/test/jacobian.jl
index bf1c42ae..fccf1826 100644
--- a/test/jacobian.jl
+++ b/test/jacobian.jl
@@ -42,10 +42,9 @@ end
     sys = structural_simplify(sys)
 
     prob = ODEProblem(sys, [1.0, 1.0], (0.0, 1.0); jac=true)
-    ode = CoupledODEs(prob)
+    ds = CoupledODEs(prob)
 
-    jac = jacobian(ode)
-    @test jac.jac_oop isa RuntimeGeneratedFunction
+    jac = jacobian(ds)
     @test jac([1.0, 1.0], [], 0.0) == [3 0;0 -3]
 
     @testset "CoupledSDEs" begin

From 00c0fd04467f4a820ca6d2dfd4bdf4184f643110 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Thu, 3 Apr 2025 17:19:21 +0100
Subject: [PATCH 17/19] fix jacian?

---
 test/jacobian.jl | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/test/jacobian.jl b/test/jacobian.jl
index fccf1826..a7d114e0 100644
--- a/test/jacobian.jl
+++ b/test/jacobian.jl
@@ -59,7 +59,6 @@ end
         sde = CoupledSDEs(prob)
 
         jac = jacobian(sde)
-        @test jac.jac_oop isa RuntimeGeneratedFunction
-        @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

From 0fca0d513ed71ab720ea8d2147a9faf9dd83c662 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Thu, 3 Apr 2025 18:05:05 +0100
Subject: [PATCH 18/19] add back the type check

only way to ensure the symbolic jacobian was extracted.
---
 test/jacobian.jl | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/test/jacobian.jl b/test/jacobian.jl
index a7d114e0..56503c7c 100644
--- a/test/jacobian.jl
+++ b/test/jacobian.jl
@@ -30,7 +30,7 @@ end
 
 @testset "MTK Jacobian" begin
     using ModelingToolkit
-    using ModelingToolkit: Num, RuntimeGeneratedFunctions.RuntimeGeneratedFunction
+    using ModelingToolkit: Num, GeneratedFunctionWrapper
     using DynamicalSystemsBase: SciMLBase
     @independent_variables t
     @variables u(t)[1:2]
@@ -45,6 +45,7 @@ end
     ds = CoupledODEs(prob)
 
     jac = jacobian(ds)
+    @test jac isa GeneratedFunctionWrapper
     @test jac([1.0, 1.0], [], 0.0) == [3 0;0 -3]
 
     @testset "CoupledSDEs" begin
@@ -59,6 +60,7 @@ end
         sde = CoupledSDEs(prob)
 
         jac = jacobian(sde)
+        @test jac isa GeneratedFunctionWrapper
         @test jac([1.0, 1.0], [], 0.0) == [3 0; 0 -3]
     end
 end

From c2a7d97cfa36885df4275f48e36362d95b46ea98 Mon Sep 17 00:00:00 2001
From: Datseris <datseris.george@gmail.com>
Date: Thu, 3 Apr 2025 19:10:38 +0100
Subject: [PATCH 19/19] correct version

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 35b8d9fa..ebe6a4c3 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.1"
+version = "3.14.0"
 
 [deps]
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"