diff --git a/docs/Project.toml b/docs/Project.toml
index 08c3ed624d..96d9feed0e 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -26,7 +26,12 @@ OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b"
 OptimizationNLopt = "4e6fcdb7-1186-4e1f-a706-475e75c168bb"
 OptimizationOptimJL = "36348300-93cb-4f02-beb5-3c3902f8871e"
 OptimizationOptimisers = "42dfb2eb-d2b4-4451-abcd-913932933ac1"
-OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqBDF = "6ad6398a-0878-4a85-9266-38940aa047c8"
+OrdinaryDiffEqDefault = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
+OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
+OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
+OrdinaryDiffEqVerner = "79d7bb75-1356-48c1-b8c0-6832512096c2"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -64,7 +69,12 @@ OptimizationBBO = "0.4"
 OptimizationNLopt = "0.3"
 OptimizationOptimJL = "0.4"
 OptimizationOptimisers = "0.3"
-OrdinaryDiffEq = "6.80.1"
+OrdinaryDiffEqBDF = "1"
+OrdinaryDiffEqDefault = "1"
+OrdinaryDiffEqRosenbrock = "1"
+OrdinaryDiffEqSDIRK = "1"
+OrdinaryDiffEqTsit5 = "1"
+OrdinaryDiffEqVerner = "1"
 Plots = "1.40"
 QuasiMonteCarlo = "0.3"
 SciMLBase = "2.46"
diff --git a/docs/src/api.md b/docs/src/api.md
index 564af6bb23..573076851e 100644
--- a/docs/src/api.md
+++ b/docs/src/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, OrdinaryDiffEq, StochasticDiffEq, JumpProcesses, Plots
+using Catalyst, OrdinaryDiffEqTsit5, StochasticDiffEq, JumpProcesses, Plots
 t = default_t()
 @parameters β γ
 @species S(t) I(t) R(t)
@@ -393,4 +393,4 @@ Finally, we provide the following helper functions to plot and animate spatial l
 lattice_plot
 lattice_animation
 lattice_kymograph
-```
\ No newline at end of file
+```
diff --git a/docs/src/faqs.md b/docs/src/faqs.md
index aba6e739f0..095031c6c2 100644
--- 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, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqTsit5, Plots
 rn = @reaction_network ABtoC begin
   (k₊,k₋), A + B <--> C
 end
@@ -132,7 +132,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, OrdinaryDiffEq
+using Catalyst, Plots, OrdinaryDiffEqTsit5
 t = default_t()
 @parameters β=1e-4 ν=.01
 @species S(t)=999.0 I(t)=1.0 R(t)=0.0
@@ -175,7 +175,7 @@ Julia `Symbol`s corresponding to each variable/parameter to their values, or
 from ModelingToolkit symbolic variables/parameters to their values. Using
 `Symbol`s we have
 ```@example faq4
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqTsit5
 rn = @reaction_network begin
     α, S + I --> 2I
     β, I --> R
diff --git a/docs/src/index.md b/docs/src/index.md
index 49533f9f43..7c4588fafe 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -32,7 +32,7 @@ etc).
 - [Steady states](@ref homotopy_continuation) (and their [stabilities](@ref steady_state_stability)) can be computed for model ODE representations.
 
 #### [Features of Catalyst composing with other packages](@id doc_index_features_composed)
-- [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) Can be used to numerically solver generated reaction rate equation ODE models.
+- [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) Can be used to numerically solve generated reaction rate equation ODE models.
 - [StochasticDiffEq.jl](https://github.com/SciML/StochasticDiffEq.jl) can be used to numerically solve generated Chemical Langevin Equation SDE models.
 - [JumpProcesses.jl](https://github.com/SciML/JumpProcesses.jl) can be used to numerically sample generated Stochastic Chemical Kinetics Jump Process models.
 - Support for [parallelization of all simulations](@ref ode_simulation_performance_parallelisation), including parallelization of [ODE simulations on GPUs](@ref ode_simulation_performance_parallelisation_GPU) using [DiffEqGPU.jl](https://github.com/SciML/DiffEqGPU.jl).
@@ -92,7 +92,7 @@ Pkg.add("Catalyst")
 
 Many Catalyst features require the installation of additional packages. E.g. for ODE-solving and simulation plotting
 ```julia
-Pkg.add("OrdinaryDiffEq")
+Pkg.add("OrdinaryDiffEqDefault")
 Pkg.add("Plots")
 ```
 is also needed.
@@ -124,7 +124,7 @@ an ordinary differential equation.
 
 ```@example home_simple_example
 # Fetch required packages.
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqDefault, 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
index 872bcc27c7..d7a2493109 100644
--- 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 also wish to install the `OrdinaryDiffEq` and `Plots` packages (for numeric simulation of models, and plotting, respectively).
+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/).
 ```julia
-Pkg.add("OrdinaryDiffEq")
+Pkg.add("OrdinaryDiffEqDefault")
 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 OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 using Plots
 ```
 Here, if we restart Julia, these `using` commands *must be rerun*.
@@ -253,4 +253,4 @@ If you are a new Julia user who has used this tutorial, and there was something
 ---
 ## References
 [^1]: [Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, Samuel A. Isaacson, *Catalyst: Fast and flexible modeling of reaction networks*, PLOS Computational Biology (2023).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530)
-[^2]: [Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah, *Julia: A Fresh Approach to Numerical Computing*, SIAM Review (2017).](https://epubs.siam.org/doi/abs/10.1137/141000671)
\ No newline at end of file
+[^2]: [Jeff Bezanson, Alan Edelman, Stefan Karpinski, Viral B. Shah, *Julia: A Fresh Approach to Numerical Computing*, SIAM Review (2017).](https://epubs.siam.org/doi/abs/10.1137/141000671)
diff --git a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md
index 11addc6057..8f1dc1d5c5 100644
--- 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("OrdinaryDiffEq")
+Pkg.add("OrdinaryDiffEqTsit5")
 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, OrdinaryDiffEq, Plots, Latexify
+using Catalyst, OrdinaryDiffEqTsit5, Plots, Latexify
 ```
 
 Let's start by using the Catalyst [`@reaction_network`](@ref) macro to specify a
@@ -386,4 +386,4 @@ A more detailed summary of the precise mathematical equations Catalyst can gener
 
 ---
 ## References
-1. [Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, Samuel A. Isaacson, *Catalyst: Fast and flexible modeling of reaction networks*, PLOS Computational Biology (2023).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530)
\ No newline at end of file
+1. [Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, Samuel A. Isaacson, *Catalyst: Fast and flexible modeling of reaction networks*, PLOS Computational Biology (2023).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530)
diff --git a/docs/src/inverse_problems/behaviour_optimisation.md b/docs/src/inverse_problems/behaviour_optimisation.md
index 81b30107c0..ef5e427e21 100644
--- a/docs/src/inverse_problems/behaviour_optimisation.md
+++ b/docs/src/inverse_problems/behaviour_optimisation.md
@@ -17,7 +17,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 OrdinaryDiffEq, Plots
+using OrdinaryDiffEqTsit5, 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]
@@ -107,4 +107,4 @@ If you use this functionality in your research, please cite the following paper
 ---
 ## References
 [^1]: [Mykel J. Kochenderfer, Tim A. Wheeler *Algorithms for Optimization*, The MIT Press (2019).](https://algorithmsbook.com/optimization/files/optimization.pdf)
-[^2]: [Lea Goentoro, Oren Shoval, Marc W Kirschner, Uri Alon *The incoherent feedforward loop can provide fold-change detection in gene regulation*, Molecular Cell (2009).](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2896310/)
\ No newline at end of file
+[^2]: [Lea Goentoro, Oren Shoval, Marc W Kirschner, Uri Alon *The incoherent feedforward loop can provide fold-change detection in gene regulation*, Molecular Cell (2009).](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2896310/)
diff --git a/docs/src/inverse_problems/examples/ode_fitting_oscillation.md b/docs/src/inverse_problems/examples/ode_fitting_oscillation.md
index f5eaa4c65f..ff3fb63805 100644
--- a/docs/src/inverse_problems/examples/ode_fitting_oscillation.md
+++ b/docs/src/inverse_problems/examples/ode_fitting_oscillation.md
@@ -4,7 +4,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 OrdinaryDiffEq
+using OrdinaryDiffEqRosenbrock
 using Optimization
 using OptimizationOptimisers # Required for the ADAM optimizer.
 using SciMLSensitivity # Required for `Optimization.AutoZygote()` automatic differentiation option.
diff --git a/docs/src/inverse_problems/global_sensitivity_analysis.md b/docs/src/inverse_problems/global_sensitivity_analysis.md
index e2a10759f9..037ce9fd74 100644
--- a/docs/src/inverse_problems/global_sensitivity_analysis.md
+++ b/docs/src/inverse_problems/global_sensitivity_analysis.md
@@ -22,7 +22,7 @@ end
 ```
 We will study the peak number of infected cases's ($max(I(t))$) sensitivity to the system's three parameters. We create a function which simulates the system from a given initial condition and measures this property:
 ```@example gsa_1
-using OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 
 u0 = [:S => 999.0, :I => 1.0, :E => 0.0, :R => 0.0]
 p_dummy = [:β => 0.0, :a => 0.0, :γ => 0.0]
@@ -157,4 +157,4 @@ If you use this functionality in your research, [in addition to Catalyst](@ref d
 
 ---
 ## References
-[^1]: [Saltelli, A et al. *Global Sensitivity Analysis. The Primer*, Wiley (2008).](http://www.andreasaltelli.eu/file/repository/A_Saltelli_Marco_Ratto_Terry_Andres_Francesca_Campolongo_Jessica_Cariboni_Debora_Gatelli_Michaela_Saisana_Stefano_Tarantola_Global_Sensitivity_Analysis_The_Primer_Wiley_Interscience_2008_.pdf)
\ No newline at end of file
+[^1]: [Saltelli, A et al. *Global Sensitivity Analysis. The Primer*, Wiley (2008).](http://www.andreasaltelli.eu/file/repository/A_Saltelli_Marco_Ratto_Terry_Andres_Francesca_Campolongo_Jessica_Cariboni_Debora_Gatelli_Michaela_Saisana_Stefano_Tarantola_Global_Sensitivity_Analysis_The_Primer_Wiley_Interscience_2008_.pdf)
diff --git a/docs/src/inverse_problems/optimization_ode_param_fitting.md b/docs/src/inverse_problems/optimization_ode_param_fitting.md
index f0226a5a14..82cff5cf1c 100644
--- a/docs/src/inverse_problems/optimization_ode_param_fitting.md
+++ b/docs/src/inverse_problems/optimization_ode_param_fitting.md
@@ -21,7 +21,7 @@ u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0]
 ps_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5]
 
 # Generate synthetic data.
-using OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 oprob_true = ODEProblem(rn, u0, (0.0, 10.0), ps_true)
 true_sol = solve(oprob_true)
 data_sol = solve(oprob_true; saveat=1.0)
@@ -39,7 +39,7 @@ Catalyst.PNG(plot(plt; fmt = :png, dpi = 200)) # hide
 
 Next, we will use DiffEqParamEstim to build a loss function to measure how well our model's solutions fit the data.
 ```@example diffeq_param_estim_1 
-using DiffEqParamEstim, Optimization
+using DiffEqParamEstim, Optimization, OrdinaryDiffEqTsit5
 ps_dummy = [:kB => 0.0, :kD => 0.0, :kP => 0.0]
 oprob = ODEProblem(rn, u0, (0.0, 10.0), ps_dummy)
 loss_function = build_loss_objective(oprob, Tsit5(), L2Loss(data_ts, data_vals), Optimization.AutoForwardDiff(); 
diff --git a/docs/src/model_creation/conservation_laws.md b/docs/src/model_creation/conservation_laws.md
index 383b85812c..c3cffa1840 100644
--- a/docs/src/model_creation/conservation_laws.md
+++ b/docs/src/model_creation/conservation_laws.md
@@ -10,7 +10,7 @@ end
 ```
 If we simulate it, we note that while the concentrations of $X₁$ and $X₂$ change throughout the simulation, the total concentration of $X$ ($= X₁ + X₂$) is constant:
 ```@example conservation_laws
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:X₁ => 80.0, :X₂ => 20.0]
 ps = [:k₁ => 10.0, :k₂ => 2.0]
 oprob = ODEProblem(rs, u0, (0.0, 1.0), ps)
@@ -47,7 +47,7 @@ Here, Catalyst encodes all conserved quantities in a single, [vector-valued](@re
 
 Practically, the `remove_conserved = true` argument can be provided when a `ReactionSystem` is converted to an `ODEProblem`:
 ```@example conservation_laws
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:X₁ => 80.0, :X₂ => 20.0]
 ps = [:k₁ => 10.0, :k₂ => 2.0]
 oprob = ODEProblem(rs, u0, (0.0, 1.0), ps; remove_conserved = true)
@@ -88,4 +88,4 @@ Finally, the `conservationlaws` function yields a $m \times n$ matrix, where $n$
 ```@example conservation_laws
 conservationlaws(rs)
 ```
-I.e. in this case we have a single conserved quantity, which contains a single copy each of the system's two species.
\ No newline at end of file
+I.e. in this case we have a single conserved quantity, which contains a single copy each of the system's two species.
diff --git a/docs/src/model_creation/constraint_equations.md b/docs/src/model_creation/constraint_equations.md
index e9a66a31b8..8df1e51391 100644
--- a/docs/src/model_creation/constraint_equations.md
+++ b/docs/src/model_creation/constraint_equations.md
@@ -28,7 +28,7 @@ creating these two systems.
 Here, to create differentials with respect to time (for our differential equations), we must import the time differential operator from Catalyst. We do this through `D = default_time_deriv()`. Here, `D(V)` denotes the differential of the variable `V` with respect to time.
 
 ```@example ceq1
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqTsit5, Plots
 t = default_t()
 D = default_time_deriv()
 
@@ -73,7 +73,7 @@ plot(sol)
 As an alternative to the previous approach, we could have constructed our
 `ReactionSystem` all at once by directly using the symbolic interface:
 ```@example ceq2
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqTsit5, Plots
 t = default_t()
 D = default_time_deriv()
 
@@ -106,7 +106,7 @@ callback interface is illustrated [here](https://docs.sciml.ai/DiffEqDocs/stable
 
 Let's first create our equations and unknowns/species again
 ```@example ceq3
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqTsit5, Plots
 t = default_t()
 D = default_time_deriv()
 
@@ -172,4 +172,4 @@ plot(sol)
 For a detailed discussion on how to directly use the lower-level but more
 flexible DifferentialEquations.jl event/callback interface, see the
 [tutorial](https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Event-handling-using-callbacks)
-on event handling using callbacks.
\ No newline at end of file
+on event handling using callbacks.
diff --git a/docs/src/model_creation/dsl_advanced.md b/docs/src/model_creation/dsl_advanced.md
index a04a632aa9..ce4cb43f7f 100644
--- a/docs/src/model_creation/dsl_advanced.md
+++ b/docs/src/model_creation/dsl_advanced.md
@@ -100,7 +100,7 @@ end
 ```
 Next, if we simulate the model, we do not need to provide values for species or parameters that have default values. In this case all have default values, so both `u0` and `ps` can be empty vectors:
 ```@example dsl_advanced_defaults
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = []
 tspan = (0.0, 10.0)
 p = []
@@ -142,6 +142,7 @@ end
 ```
 Please note that as the parameter `X₀` does not occur as part of any reactions, Catalyst's DSL cannot infer whether it is a species or a parameter. This must hence be explicitly declared. We can now simulate our model while providing `X`'s value through the `X₀` parameter:
 ```@example dsl_advanced_defaults
+using OrdinaryDiffEqTsit5
 u0 = []
 p = [:X₀ => 1.0, :p => 1.0, :d => 0.5]
 oprob = ODEProblem(rn, u0, tspan, p)
@@ -262,7 +263,7 @@ end
 ```
 Now, we can also declare our initial conditions and parameter values as vectors as well:
 ```@example dsl_advanced_vector_variables
-using OrdinaryDiffEq, Plots # hide
+using OrdinaryDiffEqDefault, Plots # hide
 u0 = [:X => [0.0, 2.0]]
 tspan = (0.0, 1.0)
 ps = [:k => [1.0, 2.0]]
@@ -332,7 +333,7 @@ end
 ```
 We can now simulate our model using normal syntax (initial condition values for observables should not, and can not, be provided):
 ```@example dsl_advanced_observables
-using OrdinaryDiffEq
+using OrdinaryDiffEqTsit5
 u0 = [:X => 1.0, :Y => 2.0, :XY => 0.0]
 tspan = (0.0, 10.0)
 ps = [:kB => 1.0, :kD => 1.5]
@@ -489,7 +490,7 @@ X
 ```
 Next, we can now use these to e.g. designate initial conditions and parameter values for model simulations:
 ```@example dsl_advanced_programmatic_unpack
-using OrdinaryDiffEq, Plots # hide
+using OrdinaryDiffEqDefault, Plots # hide
 u0 = [X => 0.1]
 tspan = (0.0, 10.0)
 ps = [p => 1.0, d => 0.2]
@@ -538,4 +539,4 @@ nothing # hide
 ```
 
 !!! note
-    When using interpolation, expressions like `2$spec` won't work; the multiplication symbol must be explicitly included like `2*$spec`.
\ No newline at end of file
+    When using interpolation, expressions like `2$spec` won't work; the multiplication symbol must be explicitly included like `2*$spec`.
diff --git a/docs/src/model_creation/examples/basic_CRN_library.md b/docs/src/model_creation/examples/basic_CRN_library.md
index 618e195b1c..4d995b59af 100644
--- a/docs/src/model_creation/examples/basic_CRN_library.md
+++ b/docs/src/model_creation/examples/basic_CRN_library.md
@@ -18,7 +18,7 @@ nothing # hide
 ```
 We can now simulate our model using all three interpretations. First, we perform a reaction rate equation-based ODE simulation:
 ```@example crn_library_birth_death
-using OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 oprob = ODEProblem(bd_process, u0, tspan, ps)
 osol = solve(oprob)
 nothing # hide
@@ -52,7 +52,7 @@ plot(oplt, splt, jplt; lw = 3, size=(800,700), layout = (3,1))
 ## [Two-state model](@id basic_CRN_library_two_states)
 The two-state model describes a component (here called $X$) which can exist in two different forms (here called $X₁$ and $X₂$). It switches between these forms at linear rates. First, we simulate the model using both ODEs and SDEs:
 ```@example crn_library_two_states
-using Catalyst, OrdinaryDiffEq, StochasticDiffEq, Plots
+using Catalyst, OrdinaryDiffEqDefault, StochasticDiffEq, Plots
 two_state_model = @reaction_network begin
     (k₁,k₂), X₁ <--> X₂
 end
@@ -96,7 +96,7 @@ u0 = [:S => 301, :E => 100, :SE => 0, :P => 0]
 tspan = (0., 100.)
 ps = [:kB => 0.00166, :kD => 0.0001, :kP => 0.1]
 
-using OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 oprob = ODEProblem(mm_system, u0, tspan, ps)
 osol  = solve(oprob)
 
@@ -132,7 +132,7 @@ end
 ```
 First, we perform a deterministic ODE simulation:
 ```@example crn_library_sir
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:S => 99, :I => 1, :R => 0]
 tspan = (0.0, 500.0)
 ps = [:α => 0.001, :β => 0.01]
@@ -179,7 +179,7 @@ Below, we perform a simple deterministic ODE simulation of the system. Next, we
 
 In two separate plots.
 ```@example crn_library_cc
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:S₁ => 1.0, :C => 0.05, :S₂ => 1.2, :S₁C => 0.0, :CP => 0.0, :P => 0.0]
 tspan = (0., 15.)
 ps = [:k₁ => 5.0, :k₂ => 5.0, :k₃ => 100.0]
@@ -206,7 +206,7 @@ end
 ```
 We can simulate the model for two different initial conditions, demonstrating the existence of two different stable steady states.
 ```@example crn_library_wilhelm
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0_1 = [:X => 1.5, :Y => 0.5]
 u0_2 = [:X => 2.5, :Y => 0.5]
 tspan = (0., 10.)
@@ -234,7 +234,7 @@ A simple example of such a loop is a transcription factor which activates its ow
 
 We simulate the self-activation loop from a single initial condition using both deterministic (ODE) and stochastic (jump) simulations. We note that while the deterministic simulation reaches a single steady state, the stochastic one switches between two different states.
 ```@example crn_library_self_activation
-using JumpProcesses, OrdinaryDiffEq, Plots
+using JumpProcesses, OrdinaryDiffEqDefault, Plots
 u0 = [:X => 4]
 tspan = (0.0, 1000.0)
 ps = [:v₀ => 0.1, :v => 2.0, :K => 10.0, :n => 2, :d => 0.1]
@@ -267,7 +267,7 @@ end
 ```
 It is generally known to (for reaction rate equation-based ODE simulations) produce oscillations when $B > 1 + A^2$. However, this result is based on models generated when *combinatorial adjustment of rates is not performed*. Since Catalyst [automatically perform these adjustments](@ref introduction_to_catalyst_ratelaws), and one reaction contains a stoichiometric constant $>1$, the threshold will be different. Here, we trial two different values of $B$. In both cases, $B < 1 + A^2$, however, in the second case the system can generate oscillations.
 ```@example crn_library_brusselator
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:X => 1.0, :Y => 1.0]
 tspan = (0., 50.)
 ps1 = [:A => 1.0, :B => 1.0]
@@ -296,7 +296,7 @@ end
 ```
 Whether the Repressilator oscillates or not depends on its parameter values. Here, we will perform deterministic (ODE) simulations for two different values of $K$, showing that it oscillates for one value and not the other one. Next, we will perform stochastic (SDE) simulations for both $K$ values, showing that the stochastic model can sustain oscillations in both cases. This is an example of the phenomena of *noise-induced oscillation*.
 ```@example crn_library_brusselator
-using OrdinaryDiffEq, StochasticDiffEq, Plots
+using OrdinaryDiffEqDefault, StochasticDiffEq, Plots
 u0 = [:X => 50.0, :Y => 15.0, :Z => 15.0]
 tspan = (0., 200.)
 ps1 = [:v => 10.0, :K => 20.0, :n => 3, :d => 0.1]
@@ -337,7 +337,7 @@ end
 ```
 Here we simulate the model for a single initial condition, showing both time-state space and phase space how it reaches a [*strange attractor*](https://www.dynamicmath.xyz/strange-attractors/).
 ```@example crn_library_chaos
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5]
 tspan = (0.0, 50.0)
 p = [:k1 => 2.1, :k2 => 0.7, :k3 => 2.9, :k4 => 1.1, :k5 => 1.0, :k6 => 0.5, :k7 => 2.7]
@@ -348,4 +348,4 @@ plt1 = plot(sol; title = "Time-state space")
 plt2 = plot(sol; idxs = (:X, :Y, :Z), title = "Phase space")
 plot(plt1, plt2; layout = (1,2), size = (800,400))
 plot!(bottom_margin = 3Plots.Measures.mm) # hide
-```
\ No newline at end of file
+```
diff --git a/docs/src/model_creation/examples/hodgkin_huxley_equation.md b/docs/src/model_creation/examples/hodgkin_huxley_equation.md
index fd6313647a..acc913cab7 100644
--- a/docs/src/model_creation/examples/hodgkin_huxley_equation.md
+++ b/docs/src/model_creation/examples/hodgkin_huxley_equation.md
@@ -18,7 +18,7 @@ complete model.
 
 We begin by importing some necessary packages:
 ```@example hh1
-using ModelingToolkit, Catalyst, NonlinearSolve, Plots, OrdinaryDiffEq
+using ModelingToolkit, Catalyst, NonlinearSolve, Plots, OrdinaryDiffEqRosenbrock
 ```
 
 ## Building the model via the Catalyst DSL
@@ -192,4 +192,4 @@ sol = solve(oprob)
 plot(sol, idxs = V, legend = :outerright)
 ```
 
-We observe the same solutions as from our original model.
\ No newline at end of file
+We observe the same solutions as from our original model.
diff --git a/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md b/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md
index eb5180c468..5826071c50 100644
--- a/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md
+++ b/docs/src/model_creation/examples/programmatic_generative_linear_pathway.md
@@ -56,7 +56,7 @@ nothing # hide
 ```
 Next, we prepare an ODE for each model (scaling the initial concentration of $X_0$ and the value of $\tau$ appropriately for each model).
 ```@example programmatic_generative_linear_pathway_dsl
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0_n3 = [:X0 => 3*1.0, :X1 => 0.0, :X2 => 0.0, :X3 => 0.0]
 ps_n3 = [:τ => 1.0/3]
 oprob_n3 = ODEProblem(lp_n3, u0_n3, (0.0, 5.0), ps_n3)
@@ -117,7 +117,7 @@ nothing # hide
 ```
 We can now simulate linear pathways of arbitrary lengths using a simple syntax. We use this to recreate our previous result from the DSL:
 ```@example programmatic_generative_linear_pathway_generative
-using OrdinaryDiffEq, Plots # hide
+using OrdinaryDiffEqDefault, Plots # hide
 sol_n3 = solve(generate_oprob(3))
 sol_n10 = solve(generate_oprob(10))
 plot(sol_n3; idxs = :Xend, label = "n = 3")
@@ -134,4 +134,4 @@ plot!(sol_n20; idxs = :Xend, label = "n = 20")
 ## References
 [^1]: [N. Korsbo, H. Jönsson *It’s about time: Analysing simplifying assumptions for modelling multi-step pathways in systems biology*, PLoS Computational Biology (2020).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1007982)
 [^2]: [J. Metz, O. Diekmann *The Abstract Foundations of Linear Chain Trickery* (1991).](https://ir.cwi.nl/pub/1559/1559D.pdf)
-[^3]: D. Fargue *Reductibilite des systemes hereditaires a des systemes dynamiques (regis par des equations differentielles aux derivees partielles)*, Comptes rendus de l'Académie des Sciences (1973).
\ No newline at end of file
+[^3]: D. Fargue *Reductibilite des systemes hereditaires a des systemes dynamiques (regis par des equations differentielles aux derivees partielles)*, Comptes rendus de l'Académie des Sciences (1973).
diff --git a/docs/src/model_creation/model_file_loading_and_export.md b/docs/src/model_creation/model_file_loading_and_export.md
index 9f13fd1035..9884dd2afc 100644
--- a/docs/src/model_creation/model_file_loading_and_export.md
+++ b/docs/src/model_creation/model_file_loading_and_export.md
@@ -80,7 +80,7 @@ Here, .net files not only contain information regarding the reaction network its
 
 A parsed reaction network's content can then be provided to various problem types for simulation. E.g. here we perform an ODE simulation of our repressilator model:
 ```julia
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqDefault, Plots
 tspan = (0.0, 10000.0)
 oprob = ODEProblem(prn.rn, Float64[], tspan, Float64[])
 sol = solve(oprob)
@@ -106,7 +106,7 @@ prn, cbs = load_SBML("brusselator.xml", massaction = true)
 ```
 Here, while [ReactionNetworkImporters generates a `ParsedReactionSystem` only](@ref model_file_import_export_sbml_rni_net), SBMLImporter generates a `ParsedReactionSystem` (here stored in `prn`) and a [so-called `CallbackSet`](https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#CallbackSet) (here stored in `cbs`). While `prn` can be used to create various problems, when we simulate them, we must also supply `cbs`. E.g. to simulate our brusselator we use:
 ```julia
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqDefault, Plots
 tspan = (0.0, 50.0)
 oprob = ODEProblem(prn.rn, prn.u0, tspan, prn.p)
 sol = solve(oprob; callback = cbs)
@@ -165,4 +165,4 @@ If you use any of this functionality in your research, [in addition to Catalyst]
     year = {2024},
     lastchecked = {2024-06-02}
 }
-```
\ No newline at end of file
+```
diff --git a/docs/src/model_creation/parametric_stoichiometry.md b/docs/src/model_creation/parametric_stoichiometry.md
index 0eac502aae..65e6cde31b 100644
--- a/docs/src/model_creation/parametric_stoichiometry.md
+++ b/docs/src/model_creation/parametric_stoichiometry.md
@@ -7,7 +7,7 @@ use symbolic stoichiometries, and discuss several caveats to be aware of.
 Let's first consider a simple reversible reaction where the number of reactants
 is a parameter, and the number of products is the product of two parameters.
 ```@example s1
-using Catalyst, Latexify, OrdinaryDiffEq, ModelingToolkit, Plots
+using Catalyst, Latexify, OrdinaryDiffEqTsit5, ModelingToolkit, Plots
 revsys = @reaction_network revsys begin
     @parameters m::Int64 n::Int64
     k₊, m*A --> (m*n)*B
diff --git a/docs/src/model_creation/reactionsystem_content_accessing.md b/docs/src/model_creation/reactionsystem_content_accessing.md
index 5318293b99..b54666015d 100644
--- a/docs/src/model_creation/reactionsystem_content_accessing.md
+++ b/docs/src/model_creation/reactionsystem_content_accessing.md
@@ -20,7 +20,7 @@ rs.X1
 ```
 to access e.g. `X1`. This symbolic variable can be used just like those [declared using `@parameters` and `@species`](@ref programmatic_CRN_construction):
 ```@example model_accessing_symbolic_variables
-using OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 u0 = [rs.X1 => 1.0, rs.X2 => 2.0]
 ps = [rs.k1 => 2.0, rs.k2 => 4.0]
 oprob = ODEProblem(rs, u0, (0.0, 10.0), ps)
@@ -204,7 +204,7 @@ rs.kₜ
 
 Practically, when we simulate our hierarchical model, we use all of this to designate initial conditions and parameters. I.e. below we designate values for the two `Xᵢ` species and `d` parameters separately:
 ```@example model_accessing_hierarchical
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [rs.nucleus.Xᵢ => 0.0, rs.cytoplasm.Xᵢ => 0.0, rs.cytoplasm.Xₐ => 0.0]
 ps = [rs.nucleus.p => 1.0, rs.nucleus.d => 0.2, rs.cytoplasm.kₐ => 5.0, rs.cytoplasm.d => 0.2, rs.kₜ => 0.1]
 oprob = ODEProblem(rs, u0, (0.0, 10.0), ps)
@@ -247,4 +247,4 @@ Catalyst.get_rxs(rs)
 Other examples of similar pairs of functions are `get_unknowns` and `unknowns`, and `get_observed` and `observed`.
 
 !!! note
-    Due to the large number of available accessor functions, most `get_` functions are not directly exported by Catalyst. Hence, these must be used as `Catalyst.get_rxs`, rather than `get_rxs` directly. Again, a full list of accessor functions and instructions on how to use them can be found in [Catalyst's API](@ref api).
\ No newline at end of file
+    Due to the large number of available accessor functions, most `get_` functions are not directly exported by Catalyst. Hence, these must be used as `Catalyst.get_rxs`, rather than `get_rxs` directly. Again, a full list of accessor functions and instructions on how to use them can be found in [Catalyst's API](@ref api).
diff --git a/docs/src/model_simulation/ensemble_simulations.md b/docs/src/model_simulation/ensemble_simulations.md
index 79a8dce1ec..d897a0b5d5 100644
--- a/docs/src/model_simulation/ensemble_simulations.md
+++ b/docs/src/model_simulation/ensemble_simulations.md
@@ -50,7 +50,7 @@ Previously, we assumed that each simulation used the same initial conditions and
 
 Here, we first create an `ODEProblem` of our previous self-activation loop:
 ```@example ensemble
-using OrdinaryDiffEq
+using OrdinaryDiffEqTsit5
 oprob = ODEProblem(sa_model, u0, tspan, ps)
 nothing # hide
 ```
@@ -74,4 +74,4 @@ plot(sols)
 ```
 Here we see that the deterministic model (unlike the stochastic one), only activates for some initial conditions (while other tends to an inactive state).
 
-The `EnsembleProblem` constructor accepts a few additional optional options (`output_func`, `reduction`, `u_init`, and `safetycopy`), which are described in more detail [here](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#Building-a-Problem). These can be used to e.g. rerun an individual simulation which does not fulfil some condition, or extract and save only some selected information from a simulation (rather than the full simulation).
\ No newline at end of file
+The `EnsembleProblem` constructor accepts a few additional optional options (`output_func`, `reduction`, `u_init`, and `safetycopy`), which are described in more detail [here](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#Building-a-Problem). These can be used to e.g. rerun an individual simulation which does not fulfil some condition, or extract and save only some selected information from a simulation (rather than the full simulation).
diff --git a/docs/src/model_simulation/examples/interactive_brusselator_simulation.md b/docs/src/model_simulation/examples/interactive_brusselator_simulation.md
index 053ead2540..e50c9b14bd 100644
--- a/docs/src/model_simulation/examples/interactive_brusselator_simulation.md
+++ b/docs/src/model_simulation/examples/interactive_brusselator_simulation.md
@@ -9,7 +9,7 @@ Let's again use the oscillating Brusselator model, extending the basic simulatio
 
 ```@example interactive_brusselator; continued = true
 using Catalyst
-using OrdinaryDiffEq
+using OrdinaryDiffEqTsit5
 using GLMakie
 GLMakie.activate!(inline = true) # hide
 
diff --git a/docs/src/model_simulation/examples/periodic_events_simulation.md b/docs/src/model_simulation/examples/periodic_events_simulation.md
index c4cf4ce2cf..37cd3bddbb 100644
--- a/docs/src/model_simulation/examples/periodic_events_simulation.md
+++ b/docs/src/model_simulation/examples/periodic_events_simulation.md
@@ -14,7 +14,7 @@ end
 ```
 We can now simulate this model, observing how a 24-hour cycle is reached
 ```@example periodic_event_example
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [:X => 150.0, :Xᴾ => 50.0]
 ps = [:kₚ => 0.1, :kᵢ => 0.1, :l => 1.0]
 tspan = (0.0, 100.0)
diff --git a/docs/src/model_simulation/ode_simulation_performance.md b/docs/src/model_simulation/ode_simulation_performance.md
index 0dbe9877a0..9e546a5732 100644
--- a/docs/src/model_simulation/ode_simulation_performance.md
+++ b/docs/src/model_simulation/ode_simulation_performance.md
@@ -15,7 +15,7 @@ Generally, ODE problems can be categorised into [*stiff ODEs* and *non-stiff ODE
 
 Here we simulate the (stiff) [Brusselator](@ref basic_CRN_library_brusselator) model using the `Tsit5` solver (which is designed for non-stiff ODEs):
 ```@example ode_simulation_performance_1
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqTsit5, Plots
 
 brusselator = @reaction_network begin
     A, ∅ --> X
@@ -40,6 +40,7 @@ sol1.retcode
 ```
 Next, we instead try the `Rodas5P` solver (which is designed for stiff problems):
 ```@example ode_simulation_performance_1
+using OrdinaryDiffEqRosenbrock
 sol2 = solve(oprob, Rodas5P())
 plot(sol2)
 ```
@@ -56,7 +57,7 @@ Finally, we should note that stiffness is not tied to the model equations only.
 ## [ODE solver selection](@id ode_simulation_performance_solvers)
 OrdinaryDiffEq implements an unusually large number of ODE solvers, with the performance of the simulation heavily depending on which one is chosen. These are provided as the second argument to the `solve` command, e.g. here we use the `Tsit5` solver to simulate a simple [birth-death process](@ref basic_CRN_library_bd):
 ```@example ode_simulation_performance_2
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqTsit5
 
 bd_model = @reaction_network begin
     (p,d), 0 <--> X
@@ -69,16 +70,19 @@ oprob = ODEProblem(bd_model, u0, tspan, ps)
 solve(oprob, Tsit5())
 nothing # hide
 ```
-If no solver argument is provided to `solve`, one is automatically selected:
+If no solver argument is provided to `solve`, and the `OrdinaryDiffEqDefault` sub-library or meta `OrdinaryDiffEq` library are loaded, then one is automatically selected:
 ```@example ode_simulation_performance_2
+using OrdinaryDiffEqDefault
 solve(oprob)
 nothing # hide
 ```
+
 While the default choice is typically enough for most single simulations, if performance is important, it can be worthwhile exploring the available solvers to find one that is especially suited for the given problem. A complete list of possible ODE solvers, with advice on optimal selection, can be found [here](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/). This section will give some general advice.
 
 The most important part of solver selection is to select one appropriate for [the problem's stiffness](@ref ode_simulation_performance_stiffness). Generally, the `Tsit5` solver is good for non-stiff problems, and `Rodas5P` for stiff problems. For large stiff problems (with many species), `FBDF` can be a good choice. We can illustrate the impact of these choices by simulating our birth-death process using the `Tsit5`, `Vern7` (an explicit solver yielding [low error in the solution](@ref ode_simulation_performance_error)), `Rodas5P`, and `FBDF` solvers (benchmarking their respective performance using [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl)):
 ```julia
-using BenchmarkTools
+using BenchmarkTools 
+using OrdinaryDiffEqTsit5, OrdinaryDiffEqRosenbrock, OrdinaryDiffEqVerner, OrdinaryDiffEqBDF
 @btime solve(oprob, Tsit5())
 @btime solve(oprob, Vern7())
 @btime solve(oprob, Rodas5P())
@@ -106,7 +110,7 @@ By default, OrdinaryDiffEq computes the Jacobian using [*automatic differentiati
 
 To use this option, simply set `jac = true` when constructing an `ODEProblem`:
 ```@example ode_simulation_performance_3
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqDefault
 
 brusselator = @reaction_network begin
     A, ∅ --> X
@@ -132,7 +136,7 @@ nothing # hide
 ### [Linear solver selection](@id ode_simulation_performance_symbolic_jacobian_linear_solver)
 When implicit solvers use e.g. the Newton-Raphson method to (at each simulation time step) solve a (typically non-linear) equation, they actually solve a linearised version of this equation. For this, they use a linear solver, the choice of which can impact performance. To specify one, we use the `linsolve` option (given to the solver function, *not* the `solve` command). E.g. to use the `KLUFactorization` linear solver (which requires loading the [LinearSolve.jl](https://github.com/SciML/LinearSolve.jl) package) we run
 ```@example ode_simulation_performance_3
-using LinearSolve
+using LinearSolve, OrdinaryDiffEqRosenbrock
 solve(oprob, Rodas5P(linsolve = KLUFactorization()))
 nothing # hide
 ```
@@ -174,7 +178,7 @@ Generally, the use of preconditioners is only recommended for advanced users who
 ## [Elimination of system conservation laws](@id ode_simulation_performance_conservation_laws)
 Previously, we have described how Catalyst, when it generates ODEs, is able to [detect and eliminate conserved quantities](@ref conservation_laws). In certain cases, doing this can improve performance. E.g. in the following example we will eliminate the single conserved quantity in a [two-state model](@ref basic_CRN_library_two_states). This results in a differential algebraic equation with a single differential equation and a single algebraic equation (as opposed to two differential equations). However, as the algebraic equation is fully determined by the ODE solution, Catalyst moves it to be an observable and our new system therefore only contains one ODE that must be solved numerically. Conservation laws can be eliminated by providing the `remove_conserved = true` option to `ODEProblem`:
 ```@example ode_simulation_performance_conservation_laws
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqTsit5
 
 # Declare model.
 rs = @reaction_network begin
@@ -209,7 +213,7 @@ end
 ```
 The model can be simulated, showing how $P$ is produced from $S$:
 ```@example ode_simulation_performance_4
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqTsit5, Plots
 u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0]
 tspan = (0.0, 50.0)
 ps = [:kB => 1.0, :kD => 0.1, :kP => 0.5, :d => 0.1]
@@ -296,7 +300,7 @@ Which backend package you should use depends on your available hardware, with th
 
 Next, we declare our model and `ODEProblem`. However, we make all values `Float64` (by appending `f0` to them) and all vectors static (by adding `@SVector` before their declaration, something which requires the [StaticArrays](https://github.com/JuliaArrays/StaticArrays.jl) package).
 ```@example ode_simulation_performance_5
-using Catalyst, OrdinaryDiffEq, StaticArrays
+using Catalyst, OrdinaryDiffEqDefault, StaticArrays
 
 mm_model = @reaction_network begin
     kB, S + E --> SE
@@ -306,7 +310,7 @@ mm_model = @reaction_network begin
 end
 @unpack S, E, SE, P, kB, kD, kP, d = mm_model
 
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = @SVector [S => 1.0f0, E => 1.0f0, SE => 0.0f0, P => 0.0f0]
 tspan = (0.0f0, 50.0f0)
 p = @SVector [kB => 1.0f0, kD => 0.1f0, kP => 0.5f0, d => 0.1f0]
@@ -358,4 +362,4 @@ If you use GPU simulations in your research, please cite the following paper to
 
 ---
 ## References
-[^1]: [E. Hairer, G. Wanner, *Solving Ordinary Differential Equations II*, Springer (1996).](https://link.springer.com/book/10.1007/978-3-642-05221-7)
\ No newline at end of file
+[^1]: [E. Hairer, G. Wanner, *Solving Ordinary Differential Equations II*, Springer (1996).](https://link.springer.com/book/10.1007/978-3-642-05221-7)
diff --git a/docs/src/model_simulation/simulation_introduction.md b/docs/src/model_simulation/simulation_introduction.md
index e2c453708b..2a70d542b4 100644
--- a/docs/src/model_simulation/simulation_introduction.md
+++ b/docs/src/model_simulation/simulation_introduction.md
@@ -99,7 +99,7 @@ nothing # hide
 ```
 Next, we can simulate the model (requires loading the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package). Simulations are performed using the `solve` function.
 ```@example simulation_intro_ode
-using OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 sol = solve(oprob)
 nothing # hide
 ```
@@ -118,6 +118,7 @@ Some additional considerations:
 ### [Designating solvers and solver options](@id simulation_intro_solver_options)
 While good defaults are generally selected, OrdinaryDiffEq enables the user to customise simulations through a long range of options that can be provided to the `solve` function. This includes specifying a [solver algorithm](https://en.wikipedia.org/wiki/Numerical_methods_for_ordinary_differential_equations), which can be provided as a second argument to `solve` (if none is provided, a suitable choice is automatically made). E.g. here we specify that the `Rodas5P` method should be used:
 ```@example simulation_intro_ode
+using OrdinaryDiffEqRosenbrock
 sol = solve(oprob, Rodas5P())
 nothing # hide
 ```
@@ -348,7 +349,7 @@ end
 ```
 This type of model will generate so called *variable rate jumps* (`VariableRateJump`s in JumpProcesses.jl). Such models can be simulated in Catalyst too, but note that now a method for time-stepping the solver must be provided to `solve`. Here ODE solvers should be given as they are used to handle integrating the explicitly time-dependent propensities for problems with variable rates, i.e. the proceeding example can be solved like
 ```@example simulation_intro_jumps
-using OrdinaryDiffEq
+using OrdinaryDiffEqTsit5
 u0map = [:P => 0]
 pmap = [:f => 1.0, :A => 2.0, :ϕ => 0.0, :d => 1.0]
 tspan = (0.0, 24.0)
@@ -405,4 +406,4 @@ For jump simulations:
 	author = {Zagatti, Guilherme Augusto and Isaacson, Samuel A. and Rackauckas, Christopher and Ilin, Vasily and Ng, See-Kiong and Bressan, Stéphane},
 	year = {2023},
 }
-```
\ No newline at end of file
+```
diff --git a/docs/src/model_simulation/simulation_plotting.md b/docs/src/model_simulation/simulation_plotting.md
index f8c40d8684..e48e10bac1 100644
--- a/docs/src/model_simulation/simulation_plotting.md
+++ b/docs/src/model_simulation/simulation_plotting.md
@@ -7,7 +7,7 @@ Catalyst uses the [Plots.jl](https://github.com/JuliaPlots/Plots.jl) package for
 ## [Common plotting options](@id simulation_plotting_options)
 Let us consider the oscillating [Brusselator](@ref basic_CRN_library_brusselator) model. We have previously shown how model simulation solutions can be plotted using the `plot` function. Here we plot an ODE simulation from the Brusselator:
 ```@example simulation_plotting
-using Catalyst, OrdinaryDiffEq, Plots
+using Catalyst, OrdinaryDiffEqDefault, Plots
 
 brusselator = @reaction_network begin
     A, ∅ --> X
@@ -84,4 +84,4 @@ The plot file type is automatically determined from the extension (if none is gi
 By default, simulations are plotted as species concentrations over time. However, [phase space](https://en.wikipedia.org/wiki/Phase_space#:~:text=In%20dynamical%20systems%20theory%20and,point%20in%20the%20phase%20space.) plots are also possible. This is done by designating the axis arguments using the `idxs` option, but providing them as a tuple. E.g. here we plot our simulation in $X-Y$ space:
 ```@example simulation_plotting
 plot(sol; idxs = (:X, :Y))
-```
\ No newline at end of file
+```
diff --git a/docs/src/model_simulation/simulation_structure_interfacing.md b/docs/src/model_simulation/simulation_structure_interfacing.md
index 1c37df3672..063f9f61ed 100644
--- a/docs/src/model_simulation/simulation_structure_interfacing.md
+++ b/docs/src/model_simulation/simulation_structure_interfacing.md
@@ -52,7 +52,7 @@ To modify a problem, the `remake` function should be used. It takes an already c
 
 Here we modify our problem to increase the initial condition concentrations of the two substrates ($S₁$ and $S₂$), and also confirm that the new problem is different from the old (unchanged) one:
 ```@example structure_indexing
-using OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 oprob_new = remake(oprob; u0 = [:S₁ => 5.0, :S₂ => 2.5])
 oprob_new != oprob
 ```
@@ -150,7 +150,7 @@ get_S(oprob)
 ## [Interfacing using symbolic representations](@id simulation_structure_interfacing_symbolic_representation)
 When e.g. [programmatic modelling is used](@ref programmatic_CRN_construction), species and parameters can be represented as *symbolic variables*. These can be used to index a problem, just like symbol-based representations can. Here we create a simple [two-state model](@ref basic_CRN_library_two_states) programmatically, and use its symbolic variables to check, and update, an initial condition:
 ```@example structure_indexing_symbolic_variables
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqDefault
 t = default_t()
 @species X1(t) X2(t)
 @parameters k1 k2
diff --git a/docs/src/spatial_modelling/lattice_reaction_systems.md b/docs/src/spatial_modelling/lattice_reaction_systems.md
index 32df70adcb..6182c2b59c 100644
--- a/docs/src/spatial_modelling/lattice_reaction_systems.md
+++ b/docs/src/spatial_modelling/lattice_reaction_systems.md
@@ -47,7 +47,7 @@ In this example we used non-uniform values for $X1$'s initial condition, but uni
 
 We can now simulate our model:
 ```@example spatial_intro_basics
-using OrdinaryDiffEq
+using OrdinaryDiffEqDefault
 sol = solve(oprob)
 nothing # hide
 ```
diff --git a/docs/src/spatial_modelling/lattice_simulation_plotting.md b/docs/src/spatial_modelling/lattice_simulation_plotting.md
index 9d7622e3d6..e135a637fb 100644
--- a/docs/src/spatial_modelling/lattice_simulation_plotting.md
+++ b/docs/src/spatial_modelling/lattice_simulation_plotting.md
@@ -19,7 +19,7 @@ The first two functions can be applied to [graph](@ref spatial_lattice_modelling
 ## [Animation and plotting of 1d Cartesian or masked lattice simulations](@id lattice_simulation_plotting_1d_grids)
 Let us consider a spatial simulation on a 1d Cartesian grid lattice:
 ```@example lattice_plotting_1d
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqDefault
 two_state_model = @reaction_network begin
     (k1,k2), X1 <--> X2
 end
@@ -62,7 +62,7 @@ For more information of either function, and additional optional arguments, plea
 ## [Animation and plotting of 2d Cartesian or masked lattice simulations](@id lattice_simulation_plotting_2d_grids)
 Two-dimensional lattice simulations can be plotted in the same manner as one-dimensional ones. However, instead of displaying a species's value as a line plot, it is displayed as a heatmap. E.g. here we simulate a spatial [Brusselator](@ref basic_CRN_library_brusselator) model and display the value of $X$ at a designated time point.
 ```@example lattice_plotting_2d
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqBDF
 brusselator = @reaction_network begin
     A, ∅ --> X
     1, 2X + Y --> 3X
@@ -95,7 +95,7 @@ Again, please check the API pages for the [`lattice_plot`](@ref) and [`lattice_a
 ## [Animation and plotting of graph lattice simulations](@id lattice_simulation_plotting_graphs)
 Finally, we consider lattice simulations on graph lattices. We first simulate a simple [birth-death process](@ref basic_CRN_library_bd) on a (6-node cyclic) graph lattice.
 ```@example lattice_plotting_graphs
-using Catalyst, Graphs, OrdinaryDiffEq
+using Catalyst, Graphs, OrdinaryDiffEqDefault
 rs = @reaction_network begin
     (p,d), 0 <--> X
 end
diff --git a/docs/src/spatial_modelling/lattice_simulation_structure_ interaction.md b/docs/src/spatial_modelling/lattice_simulation_structure_ interaction.md
index 36a95c7952..58f6d784d1 100644
--- a/docs/src/spatial_modelling/lattice_simulation_structure_ interaction.md	
+++ b/docs/src/spatial_modelling/lattice_simulation_structure_ interaction.md	
@@ -17,7 +17,7 @@ Furthermore, there are some cases of interfacing which are currently not support
 ## [Retrieving values from lattice simulations](@id lattice_simulation_structure_interaction_simulation_species)
 Let us consider a simulation of a [`LatticeReactionSystem`](@ref):
 ```@example lattice_struct_interaction_sims
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqDefault
 two_state_model = @reaction_network begin
     (k1,k2), X1 <--> X2
 end
@@ -57,7 +57,7 @@ lat_getu(sol, :X1, lrs; t = [0.5, 0.75])
 ## [Retrieving and updating species values in problems and integrators](@id lattice_simulation_structure_interaction_prob_int_species)
 Let us consider a spatial `ODEProblem`
 ```@example lattice_struct_interaction_prob_ints
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqDefault
 two_state_model = @reaction_network begin
     (k1,k2), X1 <--> X2
 end
diff --git a/docs/src/spatial_modelling/spatial_ode_simulations.md b/docs/src/spatial_modelling/spatial_ode_simulations.md
index 0b7253c9ba..a052f7ed08 100644
--- a/docs/src/spatial_modelling/spatial_ode_simulations.md
+++ b/docs/src/spatial_modelling/spatial_ode_simulations.md
@@ -7,7 +7,7 @@ Previously we have described [how to select ODE solvers, and how this can impact
 ## [Jacobian options for spatial ODE simulations](@id spatial_lattice_ode_simulations_jacobians)
 We have previously described how, when [implicit solvers are used to solve stiff ODEs](@ref ode_simulation_performance_stiffness), the [strategy for computing the system Jacobian](@ref ode_simulation_performance_jacobian) is important. This is especially the case for spatial simulations, where the Jacobian often is large and highly sparse. Catalyst implements special methods for spatial Jacobians. To utilise these, provide the `jac = true` argument to your `ODEProblem` when it is created (if `jac = false`, which is the default, [*automatic differentiation*](https://en.wikipedia.org/wiki/Automatic_differentiation) will be used for Jacobian computation). Here we simulate a [Brusselator](@ref basic_CRN_library_brusselator) while designating to use Catalyst's computed Jacobian:
 ```@example spatial_ode
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqBDF
 brusselator = @reaction_network begin
     A, ∅ --> X
     1, 2X + Y --> 3X
@@ -32,4 +32,4 @@ sol = solve(oprob, FBDF())
 nothing # hide
 ``` 
 
-It is possible to use `sparse = true` while `jac = false`, in which case a sparse Jacobian is computed using automatic differentiation. 
\ No newline at end of file
+It is possible to use `sparse = true` while `jac = false`, in which case a sparse Jacobian is computed using automatic differentiation. 
diff --git a/docs/src/steady_state_functionality/dynamical_systems.md b/docs/src/steady_state_functionality/dynamical_systems.md
index ab48f99bac..8be5065afc 100644
--- a/docs/src/steady_state_functionality/dynamical_systems.md
+++ b/docs/src/steady_state_functionality/dynamical_systems.md
@@ -22,7 +22,7 @@ nothing # hide
 ```
 Next, for any application of DynamicalSystems.jl, our `ODEProblem` must be converted into a so-called `CoupledODEs` structure. This is done by combining the ODE with the solver (and potential solver options) with which we wish to simulate it (just like when it is simulated using `solve`). Here, we will simply designate the `Tsit5` numeric solver (but provide no other options).
 ```@example dynamical_systems_basins
-using DynamicalSystems, OrdinaryDiffEq
+using DynamicalSystems, OrdinaryDiffEqTsit5
 ds = CoupledODEs(oprob, (alg = Tsit5(),))
 ```
 We can now compute the basins of attraction. This is done by first creating a grid that designates which subspace of phase-space we wish to investigate (here, the corresponding basin of attraction is found for every point on the grid). Next, we create a `AttractorsViaRecurrences` struct, that maps initial conditions to attractors, and then use that as input to the `basins_of_attraction` function.
@@ -69,7 +69,7 @@ end
 ```
 We can simulate the model, noting that its behaviour seems chaotic.
 ```@example dynamical_systems_lyapunov
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqRosenbrock, Plots
 
 u0 = [:X => 1.5, :Y => 1.5, :Z => 1.5]
 tspan = (0.0, 100.0)
@@ -93,6 +93,7 @@ Here, the largest exponent is positive, suggesting that the model is chaotic (or
 
 Next, we consider the [Brusselator] model. First we simulate the model for two similar initial conditions, confirming that they converge to the same limit cycle:
 ```@example dynamical_systems_lyapunov
+using OrdinaryDiffEqTsit5
 brusselator = @reaction_network begin
     A, ∅ --> X
     1, 2X + Y --> 3X
@@ -151,4 +152,4 @@ If you want to learn more about analysing dynamical systems, including chaotic b
 ## References
 [^1]: [S. H. Strogatz, *Nonlinear Dynamics and Chaos*, Westview Press (1994).](http://users.uoa.gr/~pjioannou/nonlin/Strogatz,%20S.%20H.%20-%20Nonlinear%20Dynamics%20And%20Chaos.pdf)
 [^2]: [A. M. Lyapunov, *The general problem of the stability of motion*, International Journal of Control (1992).](https://www.tandfonline.com/doi/abs/10.1080/00207179208934253)
-[^3]: [G. Datseris, U. Parlitz, *Nonlinear dynamics: A concise introduction interlaced with code*, Springer (2022).](https://link.springer.com/book/10.1007/978-3-030-91032-7)
\ No newline at end of file
+[^3]: [G. Datseris, U. Parlitz, *Nonlinear dynamics: A concise introduction interlaced with code*, Springer (2022).](https://link.springer.com/book/10.1007/978-3-030-91032-7)
diff --git a/docs/src/steady_state_functionality/nonlinear_solve.md b/docs/src/steady_state_functionality/nonlinear_solve.md
index fb4ae52f0a..4a84198ee9 100644
--- a/docs/src/steady_state_functionality/nonlinear_solve.md
+++ b/docs/src/steady_state_functionality/nonlinear_solve.md
@@ -102,11 +102,11 @@ Next, we provide these as an input to a `SteadyStateProblem`
 ssprob = SteadyStateProblem(dimer_production, u0, p)
 nothing # hide
 ```
-Finally, we can find the steady states using the `solver` command. Since `SteadyStateProblem`s are solved through forward ODE simulation, we must load the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package, and [select an ODE solver](@ref simulation_intro_solver_options). Any available ODE solver can be used, however, it has to be encapsulated by the `DynamicSS()` function. E.g. here we designate the `Rodas5P` solver:
+Finally, we can find the steady states using the `solver` command. Since `SteadyStateProblem`s are solved through forward ODE simulation, we must load the sublibrary of the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package that corresponds to the [selected ODE solver](@ref simulation_intro_solver_options). Any available ODE solver can be used, however, it has to be encapsulated by the `DynamicSS()` function. E.g. here we use the `Rodas5P` solver which is loaded from the `OrdinaryDiffEqRosenbrock` sublibrary:
 
 (which requires loading the SteadyStateDiffEq package).
 ```@example steady_state_solving_simulation
-using SteadyStateDiffEq, OrdinaryDiffEq
+using SteadyStateDiffEq, OrdinaryDiffEqRosenbrock
 solve(ssprob, DynamicSS(Rodas5P()))
 ```
 Note that, unlike for nonlinear system solving, `u0` is not just an initial guess of the solution, but the initial conditions from which the steady state simulation is carried out. This means that, for a system with multiple steady states, we can determine the steady states associated with specific initial conditions (which is not possible when the nonlinear solving approach is used). This also permits us to easily [handle the presence of conservation laws](@ref steady_state_solving_nonlinear_conservation_laws). The forward ODE simulation approach (unlike homotopy continuation and nonlinear solving) cannot find unstable steady states.
@@ -146,4 +146,4 @@ If you use this functionality in your research, [in addition to Catalyst](@ref d
 
 ---
 ## References
-[^1]: [J. Nocedal, S. J. Wright, *Numerical Optimization*, Springer (2006).](https://www.math.uci.edu/~qnie/Publications/NumericalOptimization.pdf)
\ No newline at end of file
+[^1]: [J. Nocedal, S. J. Wright, *Numerical Optimization*, Springer (2006).](https://www.math.uci.edu/~qnie/Publications/NumericalOptimization.pdf)
diff --git a/docs/src/v14_migration_guide.md b/docs/src/v14_migration_guide.md
index 070f6e5a6d..2603472fb0 100644
--- a/docs/src/v14_migration_guide.md
+++ b/docs/src/v14_migration_guide.md
@@ -40,7 +40,7 @@ Catalyst.iscomplete(rs)
 ```
 We can now go on and use our model for e.g. simulations:
 ```@example v14_migration_1
-using OrdinaryDiffEq, Plots
+using OrdinaryDiffEqDefault, Plots
 u0 = [X => 0.1]
 tspan = (0.0, 10.0)
 ps = [p => 1.0, d => 0.2]
@@ -186,7 +186,7 @@ For more details on how to query various structures for parameter and species va
 #### Modification of problems with conservation laws broken
 While it is possible to update e.g. `ODEProblem`s using the [`remake`](@ref simulation_structure_interfacing_problems_remake) function, this is currently not possible if the `remove_conserved = true` option was used. E.g. while
 ```@example v14_migration_5
-using Catalyst, OrdinaryDiffEq
+using Catalyst, OrdinaryDiffEqDefault
 rn = @reaction_network begin
     (k1,k2), X1 <--> X2
 end
@@ -212,4 +212,4 @@ In early versions of Catalyst, parameters and species were provided as vectors (
 *Users should never use vector-forms to represent parameter and species values*
 
 #### Additional deprecated functions
-The `reactionparams`, `numreactionparams`, and `reactionparamsmap` functions have been deprecated.
\ No newline at end of file
+The `reactionparams`, `numreactionparams`, and `reactionparamsmap` functions have been deprecated.
diff --git a/docs/unpublished/pdes.md b/docs/unpublished/pdes.md
index eb4fe09495..46e758bd6e 100644
--- a/docs/unpublished/pdes.md
+++ b/docs/unpublished/pdes.md
@@ -11,7 +11,7 @@ parameters in (Kim et al., J. Chem. Phys., 146, 2017).
 
 First we load the packages we'll use
 ```julia
-using Catalyst, MethodOfLines, DomainSets, OrdinaryDiffEq, Plots, Random, Distributions
+using Catalyst, MethodOfLines, DomainSets, OrdinaryDiffEqSDIRK, Plots, Random, Distributions
 using ModelingToolkit: scalarize, unwrap, operation, defaults
 ```
 
diff --git a/docs/unpublished/petab_ode_param_fitting.md b/docs/unpublished/petab_ode_param_fitting.md
index 64b61df5bf..a57e431912 100644
--- a/docs/unpublished/petab_ode_param_fitting.md
+++ b/docs/unpublished/petab_ode_param_fitting.md
@@ -20,7 +20,7 @@ u0 = [:S => 1.0, :E => 1.0, :SE => 0.0, :P => 0.0]
 p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5]
 
 # Generate synthetic data.
-using OrdinaryDiffEq
+using OrdinaryDiffEqTsit5
 oprob_true = ODEProblem(rn, u0, (0.0, 10.0), p_true)
 true_sol = solve(oprob_true, Tsit5())
 data_sol = solve(oprob_true, Tsit5(); saveat=1.0)
@@ -238,7 +238,7 @@ u0 = [:E => 1.0, :SE => 0.0, :P => 0.0]
 p_true = [:kB => 1.0, :kD => 0.1, :kP => 0.5]
 
 # Simulate data.
-using OrdinaryDiffEq
+using OrdinaryDiffEqTsit5
 t1, d1 = let 
     oprob_true = ODEProblem(rn, [:S=>1.0; u0], (0.0, 10.0), p_true)
     data_sol = solve(oprob_true, Tsit5(); saveat=1.0)
@@ -384,6 +384,7 @@ While in our basic example, we do not provide any additional information to our
 
 Here is an example, taken from the [more detailed PEtab.jl documentation](https://sebapersson.github.io/PEtab.jl/dev/Boehm/#Creating-a-PEtabODEProblem)
 ```@example petab1
+using OrdinaryDiffEqRosenbrock
 PEtabODEProblem(petab_model, ode_solver=ODESolver(Rodas5P(), abstol=1e-8, reltol=1e-8), gradient_method=:ForwardDiff, hessian_method=:ForwardDiff, verbose=false); nothing # hide
 ```
 ```julia