From 86928e0a468961b8e158046009e8f89c0745cc77 Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Sat, 20 Jun 2026 07:20:57 -0400
Subject: [PATCH 1/2] docs: floor Attractors at 1.38.2 to fix doc build

The Documentation workflow on master fails in the `dynamical_systems_basins`
`@example` block of `docs/src/steady_state_functionality/dynamical_systems.md`,
where `basins_of_attraction` errors with:

    UndefVarError: `referenced_sciml_model` not defined in `Attractors`
    @ Attractors .../src/mapping/attractor_mapping.jl:220

This is a version-skew between two DynamicalSystems-ecosystem packages.
DynamicalSystemsBase renamed `referrenced_sciml_model` -> `referenced_sciml_model`
(fixing the misspelling). Attractors 1.37.0 (the version the docs env resolved to)
calls `DynamicalSystemsBase.referenced_sciml_model` directly with no fallback, so
when paired with a DynamicalSystemsBase that only had the renamed/old symbol it threw.

Attractors 1.38.2 added an `isdefined(DynamicalSystemsBase, :referenced_sciml_model)`
guard that works with both the old (`referrenced_sciml_model`) and new
(`referenced_sciml_model`) names. Add `Attractors` as a direct docs dependency with
a `1.38.2` compat floor so the docs environment can no longer resolve to the broken
`Attractors < 1.38.2`. The page already uses `AttractorsViaRecurrences`,
`basins_of_attraction`, and `heatmap_basins_attractors` directly, so making
Attractors an explicit dependency is natural.

Verified locally on Julia 1.12: the failing `@example` block runs successfully with
Attractors 1.38.4 (resolved under the new floor), producing the two expected
fixed-point attractors with no UndefVarError.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 docs/Project.toml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/Project.toml b/docs/Project.toml
index 0547d99a4b..9bec3b0198 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
@@ -58,6 +59,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 
 [compat]
+Attractors = "1.38.2"
 BenchmarkTools = "1.5"
 BifurcationKit = "0.5, 0.8"
 CairoMakie = "0.15"

From 69276c13bd2933fc44c40a2645866c85b857fa7e Mon Sep 17 00:00:00 2001
From: ChrisRackauckas-Claude <accounts@chrisrackauckas.com>
Date: Sat, 20 Jun 2026 12:24:33 -0400
Subject: [PATCH 2/2] docs: adapt example code to latest SciML release APIs

The Documentation `build` job had many `@example` blocks failing after the
latest SciML stack (OrdinaryDiffEq v2 / SciMLBase 3 / BifurcationKit 0.8 /
SciMLLogging 2) introduced breaking API changes. Adapt the doc code:

- `verbose = false` (Bool) is rejected by `solve` (DiffEqBase
  `_process_verbose_param`). Use `verbose = SciMLLogging.None()`; add
  `SciMLLogging` as a docs dependency and import it in the affected pages.
- `Rodas5P(autodiff = false)` (Bool) is rejected; use
  `autodiff = AutoFiniteDiff()` (add `ADTypes` docs dep + import).
- `EnsembleProblem` `prob_func` is now called as `prob_func(prob, ctx)`
  (where `ctx::EnsembleContext` exposes `ctx.sim_id`), not
  `prob_func(prob, i, repeat)`.
- The OrdinaryDiffEq solver-level `precs` keyword moved to the LinearSolve
  linsolve: `Rodas5P(linsolve = KrylovJL_GMRES(precs = ...))` with the new
  2-argument `precs(A, p) -> (Pl, Pr)` signature.
- BifurcationKit renamed `PeriodicOrbitOCollProblem` to `Collocation`.
- `LBFGSB` is now exported by both `Optim` and `OptimizationLBFGSB`, causing
  an ambiguity `UndefVarError`; qualify it as `OptimizationLBFGSB.LBFGSB()`.

Each changed block was verified locally on Julia 1.12 (the Documentation
workflow's version) by running it against the resolved docs environment.

Co-Authored-By: Chris Rackauckas <accounts@chrisrackauckas.com>
Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
---
 docs/Project.toml                             |  4 +++
 .../behaviour_optimisation.md                 |  6 ++--
 .../examples/ode_fitting_oscillation.md       |  3 +-
 .../global_sensitivity_analysis.md            | 12 +++----
 .../inverse_problems/likelihood_profiler.md   |  4 +--
 .../optimization_ode_param_fitting.md         | 16 +++++-----
 .../turing_ode_param_fitting.md               | 12 +++----
 .../examples/noise_modelling_approaches.md    |  4 +--
 .../model_simulation/ensemble_simulations.md  |  7 ++---
 .../ode_simulation_performance.md             | 31 ++++++++-----------
 .../dynamical_systems.md                      |  8 ++---
 .../bifurcationkit_periodic_orbits.md         |  4 +--
 12 files changed, 55 insertions(+), 56 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 9bec3b0198..a772b36c52 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,4 +1,5 @@
 [deps]
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 Attractors = "f3fd9213-ca85-4dba-9dfd-7fc91308fec7"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
@@ -47,6 +48,7 @@ QuasiMonteCarlo = "8a4e6c94-4038-4cdc-81c3-7e6ffdb2a71b"
 ReactionNetworkImporters = "b4db0fb7-de2a-5028-82bf-5021f5cfa881"
 SBMLImporter = "210efffb-c3c8-456d-a807-6f55560b12fe"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLLogging = "a6db7da4-7206-11f0-1eab-35f2a5dbe1d1"
 SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -59,6 +61,7 @@ Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Turing = "fce5fe82-541a-59a6-adf8-730c64b5f9a0"
 
 [compat]
+ADTypes = "1.22"
 Attractors = "1.38.2"
 BenchmarkTools = "1.5"
 BifurcationKit = "0.5, 0.8"
@@ -106,6 +109,7 @@ QuasiMonteCarlo = "0.3"
 ReactionNetworkImporters = "1.1.0"
 SBMLImporter = "4.0.0"
 SciMLBase = "2.46, 3"
+SciMLLogging = "1, 2"
 SciMLSensitivity = "7.60"
 SpecialFunctions = "2.4"
 StaticArrays = "1.9"
diff --git a/docs/src/inverse_problems/behaviour_optimisation.md b/docs/src/inverse_problems/behaviour_optimisation.md
index 1f2bfa6c21..389a8ad6eb 100644
--- a/docs/src/inverse_problems/behaviour_optimisation.md
+++ b/docs/src/inverse_problems/behaviour_optimisation.md
@@ -24,7 +24,7 @@ Incoherent feedforward loops (network motifs where a single component both activ
 
 Our model consists of 3 species: $X$ (the input node), $Y$ (an intermediary), and $Z$ (the output node). In it, $X$ activates the production of both $Y$ and $Z$, with $Y$ also deactivating $Z$. When $X$ is activated, there will be a brief time window where $Y$ is still inactive, and $Z$ is activated. However, as $Y$ becomes active, it will turn $Z$ off. This creates a pulse of $Z$ activity. To trigger the system, we create [an event](@ref events), which increases the production rate of $X$ ($pX$) by a factor of $10$ at time $t = 10$.
 ```@example behaviour_optimization
-using Catalyst
+using Catalyst, SciMLLogging
 incoherent_feed_forward = @reaction_network begin
     @discrete_events [10.0] => [pX => 10*pX]
     pX, 0 --> X
@@ -52,7 +52,7 @@ 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
@@ -60,7 +60,7 @@ 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
index 3f06866f16..ea002e3e8b 100644
--- a/docs/src/inverse_problems/examples/ode_fitting_oscillation.md
+++ b/docs/src/inverse_problems/examples/ode_fitting_oscillation.md
@@ -24,6 +24,7 @@ First, we fetch the required packages.
 ```@example pe_osc_example
 using Catalyst
 using OrdinaryDiffEqRosenbrock
+using SciMLLogging
 using OptimizationBase
 using OptimizationOptimisers # Required for the ADAM optimizer.
 using SciMLSensitivity # Required for the `AutoZygote()` automatic differentiation option.
@@ -78,7 +79,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
index ea98f82263..2ed8e0b208 100644
--- a/docs/src/inverse_problems/global_sensitivity_analysis.md
+++ b/docs/src/inverse_problems/global_sensitivity_analysis.md
@@ -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, OrdinaryDiffEqDefault, SymbolicIndexingInterface, SciMLLogging
 
 # 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
 
@@ -65,7 +65,7 @@ While local sensitivities are primarily used as a subroutine of other methodolog
 ## [Basic example](@id global_sensitivity_analysis_basic_example)
 We will consider a simple [SEIR model of an infectious disease](https://en.wikipedia.org/wiki/Compartmental_models_in_epidemiology). This is an expansion of the classic [SIR model](@ref basic_CRN_library_sir) with an additional *exposed* state, $E$, denoting individuals who are latently infected but currently unable to transmit their infection to others.
 ```@example gsa_1
-using Catalyst
+using Catalyst, SciMLLogging
 seir_model = @reaction_network begin
     10^β, S + I --> E + I
     10^a, E --> I
@@ -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
diff --git a/docs/src/inverse_problems/likelihood_profiler.md b/docs/src/inverse_problems/likelihood_profiler.md
index 549ef9ba72..b3b84efcc3 100644
--- a/docs/src/inverse_problems/likelihood_profiler.md
+++ b/docs/src/inverse_problems/likelihood_profiler.md
@@ -57,7 +57,7 @@ petab_sol = calibrate_multistart(petab_prob, Optim.IPNewton(), 3)
 # Compute and plot the likelihood profiles.
 using LikelihoodProfiler, OptimizationLBFGSB
 pl_prob = ProfileLikelihoodProblem(petab_sol, petab_prob)
-pl_method = OptimizationProfiler(optimizer = LBFGSB(), stepper = FixedStep(; initial_step = 5e-3))
+pl_method = OptimizationProfiler(optimizer = OptimizationLBFGSB.LBFGSB(), stepper = FixedStep(; initial_step = 5e-3))
 pl_sol = LikelihoodProfiler.solve(pl_prob, pl_method)
 plot(pl_sol) # Parameters are plotted on a log scale by default.
 ```
@@ -193,7 +193,7 @@ Profile likelihood computation with LikelihoodProfiler consists of three steps:
 ```@example likelihood_profiler_basics
 using LikelihoodProfiler, OptimizationLBFGSB
 pl_prob = ProfileLikelihoodProblem(petab_sol, petab_prob)
-pl_method = OptimizationProfiler(optimizer = LBFGSB(), stepper = FixedStep(; initial_step = 1e-2))
+pl_method = OptimizationProfiler(optimizer = OptimizationLBFGSB.LBFGSB(), stepper = FixedStep(; initial_step = 1e-2))
 pl_sol = LikelihoodProfiler.solve(pl_prob, pl_method)
 nothing # hide
 ```
diff --git a/docs/src/inverse_problems/optimization_ode_param_fitting.md b/docs/src/inverse_problems/optimization_ode_param_fitting.md
index e3b2d87cbf..bd1cf0aeb8 100644
--- a/docs/src/inverse_problems/optimization_ode_param_fitting.md
+++ b/docs/src/inverse_problems/optimization_ode_param_fitting.md
@@ -22,7 +22,7 @@ Pkg.add("Plots")
 ```
 The following code provides a brief example of how to perform parameter fitting using the [Optimization.jl](https://github.com/SciML/Optimization.jl) package.
 ```julia
-using Catalyst, OrdinaryDiffEqDefault, OptimizationBase, OptimizationNLopt, SymbolicIndexingInterface
+using Catalyst, OrdinaryDiffEqDefault, OptimizationBase, OptimizationNLopt, SymbolicIndexingInterface, SciMLLogging
 
 # What we know: A model, an initial condition, and sampled datapoints.
 rn = @reaction_network begin
@@ -43,7 +43,7 @@ function loss(p, (oprob_base, p_setter, t_samples, X_samples))
     oprob = remake(oprob_base; p)
 
     # Simulate the model. If sucesfull, return sum-of-squares distance as loss.
-    sol = solve(oprob; saveat = t_samples, verbose = false, maxiters = 10000)
+    sol = solve(oprob; saveat = t_samples, verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return sum(abs2, sol[:X] .- X_samples)
 end
@@ -76,7 +76,7 @@ For simple parameter fitting problems (such as the one outlined below), [PEtab.j
 
 Let us consider a [Michaelis-Menten enzyme kinetics model](@ref basic_CRN_library_mm), where an enzyme ($E$) converts a substrate ($S$) into a product ($P$):
 ```@example optimization_paramfit_1
-using Catalyst
+using Catalyst, SciMLLogging
 rn = @reaction_network begin
     kB, S + E --> SE
     kD, SE --> S + E
@@ -113,14 +113,14 @@ oprob_base = ODEProblem(rn, u0, (0.0, 10.0), ps_init)
 function objective_function(p, _)
     p = Pair.([:kB, :kD, :kP], p)
     oprob = remake(oprob_base; p)
-    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = false, maxiters = 10000)
+    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return sum((sol .- data_vals) .^2)
 end
 ```
 When our optimisation algorithm searches parameter space it will likely consider many highly non-plausible parameter sets. To better handle this we:
 1. Add `maxiters = 10000` to our `solve` command. As most well-behaved ODEs can be solved in relatively few timesteps, this speeds up the optimisation procedure by preventing us from spending too much time trying to simulate (for the model) unsuitable parameter sets.
-2. Add `verbose = false` to our `solve` command. This prevents (potentially a very large number of) warnings from being printed to our output as unsuitable parameter sets are simulated.
+2. Add `verbose = SciMLLogging.None()` to our `solve` command. This prevents (potentially a very large number of) warnings from being printed to our output as unsuitable parameter sets are simulated.
 3. Add the line `SciMLBase.successful_retcode(sol) || return Inf`, which returns an infinite value for parameter sets which does not lead to successful simulations.
 
 To improve optimisation performance, rather than creating a new `ODEProblem` in each iteration, we pre-declare one which we [apply `remake` to](@ref simulation_structure_interfacing_problems_remake). We also use the `saveat = data_ts, save_idxs = :P` arguments to only save the values of the measured species at the measured time points.
@@ -194,7 +194,7 @@ In this case we simply modify our objective function to take this into account:
 function objective_function_S_P(p, _)
     p = Pair.([:kB, :kD, :kP], p)
     oprob = remake(oprob_base; p)
-    sol = solve(oprob; saveat = data_ts, save_idxs = [:S, :P], verbose = false, maxiters = 10000)
+    sol = solve(oprob; saveat = data_ts, save_idxs = [:S, :P], verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return sum((sol[:S] .- data_vals_S) .^2 + (sol[:P] .- data_vals_P) .^2)
 end
@@ -228,7 +228,7 @@ If we from previous knowledge know that $kD = 0.1$, and only want to fit the val
 function objective_function_known_kD(p, _)
     p = Pair.([:kB, :kD, :kP], [p[1], 0.1, p[2]])
     oprob = remake(oprob_base; p)
-    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = false, maxiters = 10000)
+    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return sum((sol .- data_vals) .^2)
 end
@@ -268,7 +268,7 @@ corresponds to the same true parameter values as used previously (`[:kB => 1.0,
 function objective_function_logtransformed(p, _)
     p = Pair.([:kB, :kD, :kP], 10.0 .^ p)
     oprob = remake(oprob_base; p)
-    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = false, maxiters = 10000)
+    sol = solve(oprob; saveat = data_ts, save_idxs = :P, verbose = SciMLLogging.None(), maxiters = 10000)
     SciMLBase.successful_retcode(sol) || return Inf
     return sum((sol .- data_vals) .^2)
 end
diff --git a/docs/src/inverse_problems/turing_ode_param_fitting.md b/docs/src/inverse_problems/turing_ode_param_fitting.md
index 2c4be60422..d3cc9e4c7a 100644
--- a/docs/src/inverse_problems/turing_ode_param_fitting.md
+++ b/docs/src/inverse_problems/turing_ode_param_fitting.md
@@ -23,7 +23,7 @@ Pkg.add("Turing")
 The following code provides a minimal example of how to infer parameter posteriors from data using the [Turing.jl](https://github.com/TuringLang/Turing.jl) package.
 ```julia
 # Create reaction network model (an SIR model).
-using Catalyst
+using Catalyst, SciMLLogging
 sir = @reaction_network begin
     γ, S + I --> 2I
     ν, I --> R
@@ -56,7 +56,7 @@ using Turing, MCMCChains
     # Simulate the model for parameter values γ, ν. Saves the solution at the measurement times.
     p = setp_oop(oprob_base, [γ, ν])
     oprob_base = remake(oprob_base; p)
-    sol = solve(oprob_base; saveat, verbose = false, maxiters = 10000)
+    sol = solve(oprob_base; saveat, verbose = SciMLLogging.None(), maxiters = 10000)
 
     # If simulation was unsuccessful, the likelihood is -Inf.
     if !SciMLBase.successful_retcode(sol)
@@ -96,7 +96,7 @@ In frequentist parameter fitting, we can use a cost function based on likelihood
 ## [Inferring parameter posterior distributions for an ODE model using Turing](@id turing_parameter_fitting_basic_example)
 For this example, we will consider a simple [SIR model of an infectious disease](@ref basic_CRN_library_sir).
 ```@example turing_paramfit
-using Catalyst
+using Catalyst, SciMLLogging
 sir = @reaction_network begin
     γ, S + I --> 2I
     ν, I --> R
@@ -142,7 +142,7 @@ setp_oop = SymbolicIndexingInterface.setp_oop(oprob_true, [:γ, :ν])
     # Simulate the model for parameter values γ, ν. Saves the solution at the measurement times.
     p = setp_oop(oprob_base, [γ, ν])
     oprob_base = remake(oprob_base; p)
-    sol = solve(oprob_base; saveat, verbose = false, maxiters = 10000)
+    sol = solve(oprob_base; saveat, verbose = SciMLLogging.None(), maxiters = 10000)
 
     # If simulation was unsuccessful, the likelihood is -Inf.
     if !SciMLBase.successful_retcode(sol)
@@ -159,7 +159,7 @@ nothing # hide
 ```
 
 Some specific comments regarding how we have declared the model above:
-- Like for [normal parameter fitting](@ref optimization_parameter_fitting_basics), we use the `maxiters = 10000` (to prevent spending a long time simulating unfeasible parameter sets) and `verbose = false` (to prevent unnecessary printing of warning messages) arguments to `solve`.
+- Like for [normal parameter fitting](@ref optimization_parameter_fitting_basics), we use the `maxiters = 10000` (to prevent spending a long time simulating unfeasible parameter sets) and `verbose = SciMLLogging.None()` (to prevent unnecessary printing of warning messages) arguments to `solve`.
 - Again, we need to handle parameter sets where the model cannot be successfully simulated. Here, we use `Turing.@addlogprob! -Inf` to set a non-existent likelihood, and `return nothing` to stop further evaluation of the specific parameter set.
 - Just like for normal parameter fitting, we wish to [fit parameters on a log scale](@ref optimization_parameter_fitting_log_scale). Here we do so by declaring log-scaled prior distributions.
 - Here we assume that we (correctly) know that the noise is normally distributed. However, we assume that we *do not know the standard deviation*. Instead, we make the standard deviation a third parameter whose value we infer as part of the inference process. More complicated noise formulas can be used (and are sometimes even advisable[^2]).
@@ -209,7 +209,7 @@ collect(chain.value[rand(1:n_steps), 1:3, rand(1:n_chains)])
 
 We can use this to e.g. draw $10$ random parameter sets from the posterior distribution, simulate the model for these parameter sets, and plot the resulting ensemble simulation. For this, we will create an `EnsembleProblem` from our `ODEProblem` using the approach described [here](@ref ensemble_simulations_varying_conditions).
 ```@example turing_paramfit
-function prob_func(prob, _, _)
+function prob_func(prob, ctx)
     γ, ν = collect(chain.value[rand(1:n_steps), 1:2, rand(1:n_chains)])
     remake(prob; p = [:γ => γ, :ν => ν])
 end
diff --git a/docs/src/model_creation/examples/noise_modelling_approaches.md b/docs/src/model_creation/examples/noise_modelling_approaches.md
index 1a2bd66993..2dcfbf2036 100644
--- a/docs/src/model_creation/examples/noise_modelling_approaches.md
+++ b/docs/src/model_creation/examples/noise_modelling_approaches.md
@@ -66,7 +66,7 @@ Again we will perform ensemble simulation. Instead of creating an `SDEProblem`,
 ```@example noise_modelling_approaches
 using Distributions
 p_dists = Dict([:v => Normal(10.0, 2.0), :K => Normal(20.0, 5.0), :n => Normal(3, 0.2), :d => Normal(0.1, 0.02)])
-function prob_func(prob, i, repeat)
+function prob_func(prob, ctx)
     p = [par => rand(p_dists[par]) for par in keys(p_dists)]
     return remake(prob; p)
 end
@@ -112,7 +112,7 @@ nothing # hide
 ```
 Finally, we will again perform ensemble simulations of our model. This time, at the beginning of each simulation, we will use `make_K_series` to generate a new $K$, and set this as the `K_in` parameter's value.
 ```@example noise_modelling_approaches
-function prob_func_Kin(prob, i, repeat)
+function prob_func_Kin(prob, ctx)
     p = [ps; :K_in => make_K_series()]    
     return remake(prob; p)
 end
diff --git a/docs/src/model_simulation/ensemble_simulations.md b/docs/src/model_simulation/ensemble_simulations.md
index 974412e0df..a09e9cc73d 100644
--- a/docs/src/model_simulation/ensemble_simulations.md
+++ b/docs/src/model_simulation/ensemble_simulations.md
@@ -96,13 +96,12 @@ nothing # hide
 ```
 Next, we wish to simulate the model for a range of initial conditions of $X$`. To do this we create a problem function, which takes the following arguments:
 - `prob`: The problem given to our `EnsembleProblem` (which is the problem that `prob_func` modifies in each iteration).
-- `i`: The number of this specific Monte Carlo iteration in the interval `1:trajectories`.
-- `repeat`: The iteration of the repeat of the simulation. Typically `1`, but potentially higher if [the simulation re-running option](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#Building-a-Problem) is used.
+- `ctx`: An `EnsembleContext` carrying information about the current simulation. Its `ctx.sim_id` field gives the number of this specific Monte Carlo iteration in the interval `1:trajectories`. Further fields (such as `ctx.repeat`, the repeat iteration of the simulation) are described [here](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#Building-a-Problem).
 
 Here we will use the following problem function (utilising [remake](@ref simulation_structure_interfacing_problems_remake)), which will provide a uniform range of initial concentrations of $X$:
 ```@example ensemble
-function prob_func(prob, i, repeat)
-    remake(prob; u0 = [:X => i * 5.0])
+function prob_func(prob, ctx)
+    remake(prob; u0 = [:X => ctx.sim_id * 5.0])
 end
 nothing # hide
 ```
diff --git a/docs/src/model_simulation/ode_simulation_performance.md b/docs/src/model_simulation/ode_simulation_performance.md
index 9c4f3fcd81..a00fe85654 100644
--- a/docs/src/model_simulation/ode_simulation_performance.md
+++ b/docs/src/model_simulation/ode_simulation_performance.md
@@ -178,20 +178,16 @@ When an implicit method solves a linear equation through an (iterative) matrix-f
 
 In practice, preconditioners are implemented as functions with a specific set of arguments. How to implement these is non-trivial, and we recommend reading OrdinaryDiffEq's documentation pages [here](https://docs.sciml.ai/DiffEqDocs/stable/features/linear_nonlinear/#Preconditioners:-precs-Specification) and [here](https://docs.sciml.ai/DiffEqDocs/stable/tutorials/advanced_ode_example/#Adding-a-Preconditioner). In this example, we will define an [Incomplete LU](https://en.wikipedia.org/wiki/Incomplete_LU_factorization) preconditioner (which requires the [IncompleteLU.jl](https://github.com/haampie/IncompleteLU.jl) package):
 ```@example ode_simulation_performance_3
-using IncompleteLU
-function incompletelu(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
-    if newW === nothing || newW
-        Pl = ilu(convert(AbstractMatrix, W), τ = 50.0)
-    else
-        Pl = Plprev
-    end
-    Pl, nothing
+using IncompleteLU, LinearAlgebra
+function incompletelu(A, p)
+    Pl = ilu(convert(AbstractMatrix, A), τ = 50.0)
+    return (Pl, I)
 end
 nothing # hide
 ```
-Next, `incompletelu` can be supplied to our solver using the `precs` argument:
+A preconditioner function takes the linear system's matrix (`A`) and the parameters (`p`), and returns a tuple with the left and right preconditioners (here we use no right preconditioner, hence we return the identity operator `I`). Next, `incompletelu` can be supplied to our matrix-free linear solver using its `precs` argument:
 ```@example ode_simulation_performance_3
-solve(oprob, Rodas5P(linsolve = KrylovJL_GMRES(), precs = incompletelu, concrete_jac = true))
+solve(oprob, Rodas5P(linsolve = KrylovJL_GMRES(precs = incompletelu), concrete_jac = true))
 nothing # hide
 ```
 Finally, we note that when using preconditioners with a matrix-free method (like `KrylovJL_GMRES`, which is also the only case when these are relevant), the `concrete_jac = true` argument is required.
@@ -251,17 +247,16 @@ To parallelise our simulations, we first need to create an `EnsembleProblem`. Th
 - The `ODEProblem` corresponds to the model simulation (`SDEProblem` and `JumpProblem`s can also be supplied, enabling the parallelisation of these problem types).
 - A function, `prob_func`, describing how to modify the problem for each simulation. If we wish to simulate the same, unmodified problem, in each simulation (primarily relevant for stochastic simulations), this argument is not required.
 
-Here, `prob_func` takes 3 arguments:
+Here, `prob_func` takes 2 arguments:
 - `prob`: The problem that it modifies at the start of each individual run (which will be the same as `EnsembleProblem`'s first argument).
-- `i`: The index of the specific simulation (in the array of all simulations that are performed).
-- `repeat`: The repeat of a specific simulation in the array. We will not use this option in this example, however, it is discussed in more detail [here](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#Building-a-Problem).
+- `ctx`: An `EnsembleContext` carrying information about the current simulation. Its `ctx.sim_id` field gives the index of the specific simulation (in the array of all simulations that are performed). Further fields (such as `ctx.repeat`) are described [here](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#Building-a-Problem).
 
-and output the `ODEProblem` simulated in the i'th simulation.
+and output the `ODEProblem` simulated in the `ctx.sim_id`'th simulation.
 
 Let us assume that we wish to simulate our model 100 times, for $kP = 0.01, 0.02, ..., 0.99, 1.0$. We define our `prob_func` using [`remake`](@ref simulation_structure_interfacing_problems_remake):
 ```@example ode_simulation_performance_4
-function prob_func(prob, i, repeat)
-    return remake(prob; p = [:kP => 0.01*i])
+function prob_func(prob, ctx)
+    return remake(prob; p = [:kP => 0.01*ctx.sim_id])
 end
 nothing # hide
 ```
@@ -343,8 +338,8 @@ nothing # hide
 ```
 When we declare our `prob_func` and `EnsembleProblem` we need to ensure that the updated `ODEProblem` uses `Float32`:
 ```@example ode_simulation_performance_5
-function prob_func(prob, i, repeat)
-    return remake(prob; p = [:kP => 0.0001f0*i])
+function prob_func(prob, ctx)
+    return remake(prob; p = [:kP => 0.0001f0*ctx.sim_id])
 end
 eprob = EnsembleProblem(oprob; prob_func = prob_func)
 nothing # hide
diff --git a/docs/src/steady_state_functionality/dynamical_systems.md b/docs/src/steady_state_functionality/dynamical_systems.md
index 31af340b9e..7e8c1114ca 100644
--- a/docs/src/steady_state_functionality/dynamical_systems.md
+++ b/docs/src/steady_state_functionality/dynamical_systems.md
@@ -100,11 +100,11 @@ plot(sol; idxs=(:X, :Y, :Z))
 ```
 Next, like when we [computed basins of attraction](@ref dynamical_systems_basins_of_attraction), we create a `CoupledODEs` corresponding to the model and state for which we wish to compute our Lyapunov spectrum. Like previously, `tspan` must provide some small interval (at least `(0.0, 1.0)` is recommended), but else have no impact on the computed Lyapunov spectrum.
 ```@example dynamical_systems_lyapunov
-using DynamicalSystems
-ds = CoupledODEs(oprob, (alg = Rodas5P(autodiff = false),))
+using DynamicalSystems, ADTypes
+ds = CoupledODEs(oprob, (alg = Rodas5P(autodiff = AutoFiniteDiff()),))
 nothing # hide
 ```
-Here, the `autodiff = false` argument is required when Lyapunov spectrums are computed. We can now provide our `CoupledODEs` (`ds`) to `lyapunovspectrum` to compute the lyapunov spectrum. This function requires a second argument (here set to `100`). Generally setting this to a higher value will increase accuracy, but also increase runtime (since `lyapunovspectrum` is fast for most systems, setting this to a large value is recommended).
+Here, the `autodiff = AutoFiniteDiff()` argument is required when Lyapunov spectrums are computed. We can now provide our `CoupledODEs` (`ds`) to `lyapunovspectrum` to compute the lyapunov spectrum. This function requires a second argument (here set to `100`). Generally setting this to a higher value will increase accuracy, but also increase runtime (since `lyapunovspectrum` is fast for most systems, setting this to a large value is recommended).
 ```@example dynamical_systems_lyapunov
 lyapunovspectrum(ds, 100)
 ```
@@ -134,7 +134,7 @@ plot!(osol2; idxs = (:X, :Y))
 ```
 Next, we compute the Lyapunov spectrum at one of the initial conditions:
 ```@example dynamical_systems_lyapunov
-ds = CoupledODEs(oprob1, (alg = Rodas5P(autodiff = false),))
+ds = CoupledODEs(oprob1, (alg = Rodas5P(autodiff = AutoFiniteDiff()),))
 lyapunovspectrum(ds, 100)
 ```
 Here, all Lyapunov exponents are negative, confirming that the brusselator is non-chaotic.
diff --git a/docs/src/steady_state_functionality/examples/bifurcationkit_periodic_orbits.md b/docs/src/steady_state_functionality/examples/bifurcationkit_periodic_orbits.md
index de4a094323..95ee1d7986 100644
--- a/docs/src/steady_state_functionality/examples/bifurcationkit_periodic_orbits.md
+++ b/docs/src/steady_state_functionality/examples/bifurcationkit_periodic_orbits.md
@@ -72,9 +72,9 @@ First, we set the options for the continuation. Just like for bifurcation diagra
 opts_po = ContinuationPar(opts_br)
 nothing # hide
 ```
-During the continuation we will compute the periodic orbit using the [Collocation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/#Collocation-method) method (however, the [Trapezoid](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/#Trapezoid-method) and [Shooting](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/#Shooting-method) method also exist, with each having their advantages and disadvantages). For this, we create a [`PeriodicOrbitOCollProblem`](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/library/#BifurcationKit.PeriodicOrbitOCollProblem) which we will supply to our continuation computation.
+During the continuation we will compute the periodic orbit using the [Collocation](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/#Collocation-method) method (however, the [Trapezoid](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/#Trapezoid-method) and [Shooting](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/periodicOrbit/#Shooting-method) method also exist, with each having their advantages and disadvantages). For this, we create a [`Collocation`](https://bifurcationkit.github.io/BifurcationKitDocs.jl/stable/library/#BifurcationKit.Collocation) problem which we will supply to our continuation computation.
 ```@example bifurcationkit_periodic_orbits
-poprob = PeriodicOrbitOCollProblem(50, 4)
+poprob = Collocation(50, 4)
 nothing # hide
 ```
 Finally, we will also create a `record_from_solution` function. This is a function which records information of the solution at each step of the continuation (which we can later plot or investigate using other means).