SciML · TorkelE · Jun 15, 2026 · Jun 15, 2026 · Jun 22, 2026 · Jun 22, 2026
diff --git a/Project.toml b/Project.toml
@@ -63,14 +63,10 @@ Latexify = "0.16.6"
 LinearAlgebra = "1.10"
 MacroTools = "0.5.5"
 Makie = "0.22.1, 0.23, 0.24"
-ModelingToolkitBase = "1.17"
+ModelingToolkitBase = "1.46"
 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 @@ DomainSets = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
 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",
-    "Logging", "NonlinearSolve", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault",
-    "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqVerner",
+    "Logging", "NonlinearSolve", "OrdinaryDiffEq", "OrdinaryDiffEqCore",
     "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
@@ -29,17 +29,12 @@ NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
 Optim = "429524aa-4258-5aef-a3af-852621145aeb"
 OptimizationBBO = "3e6eede4-6085-4f62-9a71-46d9bc1eb92b"
 OptimizationBase = "bca83a33-5cc9-4baa-983d-23429ab6bcbb"
+OptimizationEvolutionary = "cb963754-43f6-435e-8d4b-99009ff27753"
 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 +64,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"
@@ -85,19 +80,14 @@ NetworkLayout = "0.4"
 NonlinearSolve = "4"
 NonlinearSolveFirstOrder = "1, 2.1"
 Optim = "2"
-OptimizationBBO = "0.4"
+OptimizationBBO = "0.4.7"
 OptimizationBase = "4, 5.0"
+OptimizationEvolutionary = "0.4.9"
 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/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

diff --git a/docs/src/inverse_problems/global_sensitivity_analysis.md b/docs/src/inverse_problems/global_sensitivity_analysis.md
@@ -9,7 +9,7 @@ Pkg.activate(; temp = true) # Creates a temporary environment, which is deleted
 Pkg.add("Catalyst")
 Pkg.add("GlobalSensitivity")
 Pkg.add("Plots")
-Pkg.add("OrdinaryDiffEqDefault")
+Pkg.add("OrdinaryDiffEq")
 ```
 ```@raw html
 </details>
@@ -19,7 +19,7 @@ Pkg.add("OrdinaryDiffEqDefault")
 ```
 The following code provides a brief example of how to perform global sensitivity analysis using the [GlobalSensitivity.jl](https://github.com/SciML/GlobalSensitivity.jl) package.
 ```julia
-using Catalyst, GlobalSensitivity, OrdinaryDiffEqDefault, SymbolicIndexingInterface
+using Catalyst, GlobalSensitivity, OrdinaryDiffEq, SymbolicIndexingInterface
 
 # Designates a model, parameter set, and set of initial conditions.
 # For this infectious disease model we will determine the peak number of cases's sensitivity to the parameters.
@@ -38,7 +38,7 @@ function peak_cases(p)
     # Updates the ODEProblem with teh proposed parameter set.
     p = p_setter(oprob_base, p)
     oprob = remake(oprob_base; p)
-    sol = solve(oprob; maxiters = 100000, verbose = false)
+    sol = solve(oprob; maxiters = 100000, verbose = SciMLLogging.None())
     return maximum(sol[:I])
 end
 
@@ -74,7 +74,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 OrdinaryDiffEqDefault
+using OrdinaryDiffEq
 
 u0 = [:S => 999.0, :I => 1.0, :E => 0.0, :R => 0.0]
 p_dummy = [:β => 0.0, :a => 0.0, :γ => 0.0]
@@ -83,7 +83,7 @@ oprob_base = ODEProblem(seir_model, u0, (0.0, 10000.0), p_dummy)
 function peak_cases(p)
     ps = [:β => p[1], :a => p[2], :γ => p[3]]
     oprob = remake(oprob_base; p = ps)
-    sol = solve(oprob; maxiters = 100000, verbose = false)
+    sol = solve(oprob; maxiters = 100000, verbose = SciMLLogging.None())
     SciMLBase.successful_retcode(sol) || return Inf
     return maximum(sol[:I])
 end
@@ -106,7 +106,7 @@ on the domain $10^β ∈ (-3.0,-1.0)$, $10^a ∈ (-2.0,0.0)$, $10^γ ∈ (-2.0,0
     We should make a couple of notes about the example above:
     - Here, we write our parameters on the forms $10^β$, $10^a$, and $10^γ$, which transforms them into log-space. As [previously described](@ref optimization_parameter_fitting_log_scale), this is advantageous in the context of inverse problems such as this one.
     - For GSA, where a function is evaluated a large number of times, it is ideal to write it as performant as possible. Hence, we initially create a base `ODEProblem`, and then apply the [`remake`](@ref simulation_structure_interfacing_problems_remake) function to it in each evaluation of `peak_cases` to generate a problem which is solved for that specific parameter set.
-    - Again, as [previously described in other inverse problem tutorials](@ref optimization_parameter_fitting_basics), when exploring a function over large parameter spaces, we will likely simulate our model for unsuitable parameter sets. To reduce time spent on these, and to avoid excessive warning messages, we provide the `maxiters = 100000` and `verbose = false` arguments to `solve`.
+    - Again, as [previously described in other inverse problem tutorials](@ref optimization_parameter_fitting_basics), when exploring a function over large parameter spaces, we will likely simulate our model for unsuitable parameter sets. To reduce time spent on these, and to avoid excessive warning messages, we provide the `maxiters = 100000` and `verbose = SciMLLogging.None()` arguments to `solve`.
     - As we have encountered in [a few other cases](@ref optimization_parameter_fitting_basics), the `gsa` function is not able to take parameter inputs of the map form usually used for Catalyst. Hence, as a first step in `peak_cases` we convert the parameter vector to this form. Next, we remember that the order of the parameters when we e.g. evaluate the GSA output, or set the parameter bounds, corresponds to the order used in `ps = [:β => p[1], :a => p[2], :γ => p[3]]`.
 
 
@@ -174,7 +174,7 @@ Previously, we have demonstrated GSA on functions with scalar outputs. However,
 function peak_cases_2(p)
     ps = [:β => p[1], :a => p[2], :γ => p[3]]
     oprob = remake(oprob_base; p = ps)
-    sol = solve(oprob; maxiters = 100000, verbose = false)
+    sol = solve(oprob; maxiters = 100000, verbose = SciMLLogging.None())
     SciMLBase.successful_retcode(sol) || return Inf
     return [maximum(sol[:E]), maximum(sol[:I])]
 end