SciML · TorkelE · Jun 12, 2026 · Jun 15, 2026 · Jun 15, 2026 · Jun 16, 2026
diff --git a/Project.toml b/Project.toml
@@ -65,12 +65,8 @@ MacroTools = "0.5.5"
 Makie = "0.22.1, 0.23, 0.24"
 ModelingToolkitBase = "1.17"
 NetworkLayout = "0.4.7"
-OrdinaryDiffEqBDF = "1, 2"
 OrdinaryDiffEqCore = "3.22, 4"
-OrdinaryDiffEqDefault = "1, 2"
-OrdinaryDiffEqRosenbrock = "1, 2"
-OrdinaryDiffEqTsit5 = "1, 2"
-OrdinaryDiffEqVerner = "1, 2"
+OrdinaryDiffEq = "6, 7"
 Parameters = "0.12, 0.13"
 Reexport = "1.0"
 RuntimeGeneratedFunctions = "0.5.12"
@@ -93,12 +89,8 @@ ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
-OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
-OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
-OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
-OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
@@ -112,7 +104,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["DataInterpolations", "DiffEqCallbacks", "DiffEqNoiseProcess", "DomainSets",
-    "ExplicitImports", "Logging", "NonlinearSolve", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault",
-    "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqVerner",
+    "ExplicitImports", "Logging", "NonlinearSolve", "OrdinaryDiffEqCore", "OrdinaryDiffEq",
     "Pkg", "Plots", "Random", "SafeTestsets", "StableRNGs",
     "StaticArrays", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test"]
diff --git a/README.md b/README.md
@@ -70,7 +70,7 @@ an ordinary differential equation.
 
 ```julia
 # Fetch required packages.
-using Catalyst, OrdinaryDiffEqDefault, Plots
+using Catalyst, OrdinaryDiffEq, Plots
 
 # Create model.
 model = @reaction_network begin

diff --git a/docs/Project.toml b/docs/Project.toml
@@ -21,25 +21,22 @@ Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LikelihoodProfiler = "93acb638-a083-5915-8dce-d129bc6a3f59"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
+Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
 MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d"
 ModelingToolkitBase = "7771a370-6774-4173-bd38-47e70ca0b839"
+ModelingToolkitNeuralNets = "f162e290-f571-43a6-83d9-22ecc16da15f"
 NetworkLayout = "46757867-2c16-5918-afeb-47bfcb05e46a"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
+Optimisers = "3bd65402-5787-11e9-1adc-39752487f4e2"
 OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b"
 OptimizationBase = "bca83a33-5cc9-4baa-983d-23429ab6bcbb"
 OptimizationLBFGSB = "22f7324a-a79d-40f2-bebe-3af60c77bd15"
 OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
-OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
-OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
-OrdinaryDiffEqNonlinearSolve = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
-OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
-OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
-OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
-OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 PEtab = "48d54b35-e43e-4a66-a5a1-dde6b987cf69"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
@@ -69,7 +66,7 @@ DiffEqNoiseProcess = "5.31.1"
 Distributions = "0.25"
 Documenter = "1.11.1"
 DynamicPolynomials = "0.6"
-DynamicalSystems = "3.6.7"
+DynamicalSystems = "3.6.8"
 GlobalSensitivity = "2.6"
 GraphMakie = "0.6"
 Graphs = "1.11.1"
@@ -79,25 +76,22 @@ JumpProcesses = "9.23"
 Latexify = "0.16.5"
 LikelihoodProfiler = "1.5.1"
 LinearSolve = "2.30, 3"
+Lux = "1.31.4"
 MCMCChains = "6, 7"
 ModelingToolkitBase = "1.17"
+ModelingToolkitNeuralNets = "2.5.1"
 NetworkLayout = "0.4"
 NonlinearSolve = "4"
 NonlinearSolveFirstOrder = "1, 2.1"
 Optim = "2"
-OptimizationBBO = "0.4"
+Optimisers = "0.4.7"
+OptimizationBBO = "0.4.7"
 OptimizationBase = "4, 5.0"
 OptimizationLBFGSB = "1"
 OptimizationNLopt = "0.3"
 OptimizationOptimJL = "0.4"
 OptimizationOptimisers = "0.3"
-OrdinaryDiffEqBDF = "1, 2"
-OrdinaryDiffEqDefault = "1, 2"
-OrdinaryDiffEqNonlinearSolve = "1, 2"
-OrdinaryDiffEqRosenbrock = "1, 2"
-OrdinaryDiffEqSDIRK = "1, 2"
-OrdinaryDiffEqTsit5 = "1, 2"
-OrdinaryDiffEqVerner = "1, 2"
+OrdinaryDiffEq = "7.1.0"
 PEtab = "5"
 Plots = "1.40"
 QuasiMonteCarlo = "0.3"

diff --git a/docs/old_files/unused_tutorials/finite_state_projection_simulation.md b/docs/old_files/unused_tutorials/finite_state_projection_simulation.md
@@ -9,8 +9,7 @@ Pkg.activate(; temp = true) # Creates a temporary environment, which is deleted
 Pkg.add("Catalyst")
 Pkg.add("FiniteStateProjection")
 Pkg.add("JumpProcesses")
-Pkg.add("OrdinaryDiffEqDefault")
-Pkg.add("OrdinaryDiffEqRosenbrock")
+Pkg.add("OrdinaryDiffEq")
 Pkg.add("Plots")
 Pkg.add("SteadyStateDiffEq")
 ```
@@ -44,7 +43,7 @@ oprob = ODEProblem(fsp_sys, u0, tspan, ps)
 
 # Simulate ODE (it can be quite large, so consider performance options).
 # Plot solution as a heatmap at a specific time point.
-using OrdinaryDiffEqRosenbrock, Plots
+using OrdinaryDiffEq, Plots
 osol = solve(oprob, Rodas5P())
 heatmap(0:19, 0:19, osol(50.0); xguide = "Y", yguide = "X")
 ```
@@ -110,7 +109,7 @@ We also plot the full distribution using the `bar` function. Finally, the initia
 
 Now, we can finally create an `ODEProblem` using our `FSPSystem`, initial conditions, and the parameters declared previously. We can simulate this `ODEProblem` like any other ODE.
 ```@example state_projection_one_species
-using OrdinaryDiffEqDefault
+using OrdinaryDiffEq
 oprob = ODEProblem(fsp_sys, u0, tspan, ps)
 osol = solve(oprob)
 nothing # hide
@@ -151,7 +150,7 @@ nothing # hide
 Finally, we can simulate the model just like in the 1-dimensional case. As we are simulating an ODE with $25⋅25 = 625$ states, we need to make some considerations regarding performance. In this case, we will simply specify the `Rodas5P()` ODE solver (more extensive advice on performance can be found [here](@ref ode_simulation_performance)). Here, we perform a simulation with a long time span ($t = 100.0$), aiming to find the system's steady state distribution. Next, we plot it using the `heatmap` function.
 ```@example state_projection_multi_species
 using Plots # hide
-using OrdinaryDiffEqRosenbrock
+using OrdinaryDiffEq
 oprob = ODEProblem(fsp_sys, u0, 100.0, ps)
 osol = solve(oprob, Rodas5P())
 heatmap(0:24, 0:24, osol[end]; xguide = "X₂", yguide = "X")
@@ -163,7 +162,7 @@ heatmap(0:24, 0:24, osol[end]; xguide = "X₂", yguide = "X")
 ## [Finite state projection steady state simulations](@id state_projection_steady_state_sim)
 Previously, we have shown how the [SteadyStateDiffEq.jl](https://github.com/SciML/SteadyStateDiffEq.jl) package can be used to [find an ODE's steady state through forward simulation](@ref steady_state_stability). The same interface can be used for ODEs generated through FiniteStateProjection. Below, we use this to find the steady state of the dimerisation example studied in the last example.
 ```@example state_projection_multi_species
-using SteadyStateDiffEq, OrdinaryDiffEqRosenbrock
+using SteadyStateDiffEq, OrdinaryDiffEq
 ssprob = SteadyStateProblem(fsp_sys, u0, ps)
 sssol = solve(ssprob, DynamicSS(Rodas5P()))
 heatmap(0:24, 0:24, sssol; xguide = "X₂", yguide = "X")

diff --git a/docs/pages.jl b/docs/pages.jl
@@ -67,6 +67,7 @@ pages = Any[
     "inverse_problems/likelihood_profiler.md",
     "inverse_problems/global_sensitivity_analysis.md",
     "inverse_problems/turing_ode_param_fitting.md",
+    "inverse_problems/udes.md",
     "Examples" => Any[
         "inverse_problems/examples/ode_fitting_oscillation.md"
     ]

diff --git a/docs/src/api/core_api.md b/docs/src/api/core_api.md
@@ -35,7 +35,7 @@ corresponding chemical reaction ODE models, chemical Langevin equation SDE
 models, and stochastic chemical kinetics jump process models.
 
 ```@example ex1
-using Catalyst, OrdinaryDiffEqTsit5, StochasticDiffEq, JumpProcesses, Plots
+using Catalyst, OrdinaryDiffEq, StochasticDiffEq, JumpProcesses, Plots
 t = default_t()
 @parameters β γ
 @species S(t) I(t) R(t)

diff --git a/docs/src/faqs.md b/docs/src/faqs.md
@@ -5,7 +5,7 @@ One can directly use symbolic variables to index into SciML solution objects.
 Moreover, observables can also be evaluated in this way. For example,
 consider the system
 ```@example faq1
-using Catalyst, OrdinaryDiffEqNonlinearSolve, OrdinaryDiffEqTsit5, Plots
+using Catalyst, OrdinaryDiffEq, Plots
 rn = @reaction_network ABtoC begin
   (k₊,k₋), A + B <--> C
 end
@@ -134,7 +134,7 @@ When directly constructing a `ReactionSystem`, we can set the symbolic values to
 have the desired default values, and this will automatically be propagated
 through to the equation solvers:
 ```@example faq3
-using Catalyst, Plots, OrdinaryDiffEqTsit5
+using Catalyst, Plots, OrdinaryDiffEq
 t = default_t()
 @parameters β=1e-4 ν=.01
 @species S(t)=999.0 I(t)=1.0 R(t)=0.0
@@ -166,7 +166,7 @@ Julia `Symbol`s corresponding to each variable/parameter to their values, or
 from ModelingToolkitBase symbolic variables/parameters to their values. Using
 `Symbol`s we have
 ```@example faq4
-using Catalyst, OrdinaryDiffEqTsit5
+using Catalyst, OrdinaryDiffEq
 rn = @reaction_network begin
     α, S + I --> 2I
     β, I --> R

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -73,7 +73,7 @@ Pkg.add("Catalyst")
 
 Many Catalyst features require the installation of additional packages. E.g. for ODE-solving and simulation plotting
 ```julia
-Pkg.add("OrdinaryDiffEqDefault")
+Pkg.add("OrdinaryDiffEq")
 Pkg.add("Plots")
 ```
 is also needed.
@@ -116,7 +116,7 @@ an ordinary differential equation.
 
 ```@example home_simple_example
 # Fetch required packages.
-using Catalyst, OrdinaryDiffEqDefault, Plots
+using Catalyst, OrdinaryDiffEq, Plots
 
 # Create model.
 model = @reaction_network begin

diff --git a/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md b/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md
@@ -55,15 +55,15 @@ To import a Julia package into a session, you can use the `using PackageName` co
 using Pkg
 Pkg.add("Catalyst")
 ```
-Here, the Julia package manager package (`Pkg`) is by default installed on your computer when Julia is installed, and can be activated directly. Next, we install an ODE solver from a sub-library of the larger `OrdinaryDiffEq` package, and install the `Plots` package for making graphs. We will import the recommended default solver from the `OrdinaryDiffEqDefault` sub-library. A full list of `OrdinaryDiffEq` solver sublibraries can be found on the sidebar of [this page](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
+Here, the Julia package manager package (`Pkg`) is by default installed on your computer when Julia is installed, and can be activated directly. Next, we install the `OrdinaryDiffEq` package, which provides ODE solvers, and the `Plots` package for making graphs. A full list of available `OrdinaryDiffEq` solvers can be found in [its documentation](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
 ```julia
-Pkg.add("OrdinaryDiffEqDefault")
+Pkg.add("OrdinaryDiffEq")
 Pkg.add("Plots")
 ```
 Once a package has been installed through the `Pkg.add` command, this command does not have to be repeated if we restart our Julia session. We can now import all three packages into our current session with:
 ```@example ex2
 using Catalyst
-using OrdinaryDiffEqDefault
+using OrdinaryDiffEq
 using Plots
 ```
 Here, if we restart Julia, these `using` commands *must be rerun*.

diff --git a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
@@ -20,7 +20,7 @@ Pkg.activate("catalyst_introduction")
 
 # packages we will use in this tutorial
 Pkg.add("Catalyst")
-Pkg.add("OrdinaryDiffEqTsit5")
+Pkg.add("OrdinaryDiffEq")
 Pkg.add("Plots")
 Pkg.add("Latexify")
 Pkg.add("JumpProcesses")
@@ -29,7 +29,7 @@ Pkg.add("StochasticDiffEq")
 
 We next load the basic packages we'll need for our first example:
 ```@example tut1
-using Catalyst, OrdinaryDiffEqTsit5, Plots, Latexify
+using Catalyst, OrdinaryDiffEq, Plots, Latexify
 ```
 
 Let's start by using the Catalyst [`@reaction_network`](@ref) macro to specify a

diff --git a/docs/src/inverse_problems/behaviour_optimisation.md b/docs/src/inverse_problems/behaviour_optimisation.md
@@ -9,7 +9,7 @@ Pkg.activate(; temp = true) # Creates a temporary environment, which is deleted
 Pkg.add("Catalyst")
 Pkg.add("OptimizationBase")
 Pkg.add("OptimizationBBO")
-Pkg.add("OrdinaryDiffEqDefault")
+Pkg.add("OrdinaryDiffEq")
 Pkg.add("Plots")
 ```
 ```@raw html
@@ -35,7 +35,7 @@ end
 ```
 To demonstrate this pulsing behaviour we will simulate the system for an example parameter set. We select an initial condition (`u0`) so the system begins in a steady state.
 ```@example behaviour_optimization
-using OrdinaryDiffEqDefault, Plots
+using OrdinaryDiffEq, Plots
 example_p = [:pX => 0.1, :pY => 1.0, :pZ => 1.0]
 tspan = (0.0, 50.0)
 example_u0 = [:X => 0.1, :Y => 0.1, :Z => 1.0]
@@ -52,15 +52,15 @@ function pulse_amplitude(p, _)
     p = Dict([:pX => p[1], :pY => p[2], :pZ => p[2]])
     u0 = [:X => p[:pX], :Y => p[:pX]*p[:pY], :Z => p[:pZ]/p[:pY]^2]
     oprob_local = remake(oprob; u0, p)
-    sol = solve(oprob_local; verbose = false, maxiters = 10000)
+    sol = solve(oprob_local; verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return -(maximum(sol[:Z]) - sol[:Z][1])
 end
 nothing # hide
 ```
 This objective function takes two arguments (a parameter value `p`, and an additional one which we will ignore but is discussed in a note [here](@ref optimization_parameter_fitting_basics)). It first calculates the new initial steady state concentration for the given parameter set. Next, it creates an updated `ODEProblem` using the steady state as initial conditions and the, to the objective function provided, input parameter set.  Finally, Optimization.jl finds the function's *minimum value*, so to find the *maximum* relative pulse amplitude, we make our objective function return the negative pulse amplitude.
 
-As described [in our tutorial on parameter fitting using Optimization.jl](@ref optimization_parameter_fitting_basics) we use `remake`, `verbose = false`, `maxiters = 10000`, and a check on the simulations return code, all providing various advantages to the optimisation procedure (as explained in that tutorial).
+As described [in our tutorial on parameter fitting using Optimization.jl](@ref optimization_parameter_fitting_basics) we use `remake`, `verbose = SciMLLogging.None()`, `maxiters = 10000`, and a check on the simulations return code, all providing various advantages to the optimisation procedure (as explained in that tutorial).
 
 Just like for [parameter fitting](@ref optimization_parameter_fitting_basics), we create an `OptimizationProblem` using our objective function, and some initial guess of the parameter values. We also [set upper and lower bounds](@ref optimization_parameter_fitting_constraints) for each parameter using the `lb` and `ub` optional arguments (in this case limiting each parameter's value to the interval $(0.1,10.0)$).
 ```@example behaviour_optimization

diff --git a/docs/src/inverse_problems/examples/ode_fitting_oscillation.md b/docs/src/inverse_problems/examples/ode_fitting_oscillation.md
@@ -9,7 +9,7 @@ Pkg.activate(; temp = true) # Creates a temporary environment, which is deleted
 Pkg.add("Catalyst")
 Pkg.add("OptimizationBase")
 Pkg.add("OptimizationOptimisers")
-Pkg.add("OrdinaryDiffEqRosenbrock")
+Pkg.add("OrdinaryDiffEq")
 Pkg.add("Plots")
 Pkg.add("SciMLSensitivity")
 ```
@@ -23,7 +23,7 @@ In this example we will use [Optimization.jl](https://github.com/SciML/Optimizat
 First, we fetch the required packages.
 ```@example pe_osc_example
 using Catalyst
-using OrdinaryDiffEqRosenbrock
+using OrdinaryDiffEq
 using OptimizationBase
 using OptimizationOptimisers # Required for the ADAM optimizer.
 using SciMLSensitivity # Required for the `AutoZygote()` automatic differentiation option.
@@ -78,7 +78,7 @@ function optimize_p(pinit, tend,
         p = set_p(prob, p)
         newtimes = filter(<=(tend), sample_times)
         newprob = remake(prob; p)
-        sol = Array(solve(newprob, Rosenbrock23(); saveat = newtimes, verbose = false, maxiters = 10000))
+        sol = Array(solve(newprob, Rosenbrock23(); saveat = newtimes, verbose = SciMLLogging.None(), maxiters = 10000))
         loss = sum(abs2, sol .- sample_vals[:, 1:size(sol,2)])
         return loss
     end