SciML · ChrisRackauckas · Jul 29, 2025
diff --git a/docs/make.jl b/docs/make.jl
@@ -1,4 +1,5 @@
-using Documenter, DiffEqBase, SciMLBase, OrdinaryDiffEq, OrdinaryDiffEqBDF, OrdinaryDiffEqCore, StochasticDiffEq, DelayDiffEq, SteadyStateDiffEq, DiffEqCallbacks
+using Documenter, DiffEqBase, SciMLBase, OrdinaryDiffEq, OrdinaryDiffEqBDF,
+      OrdinaryDiffEqCore, StochasticDiffEq, DelayDiffEq, SteadyStateDiffEq, DiffEqCallbacks
 import ODEProblemLibrary,
        SDEProblemLibrary, DDEProblemLibrary, DAEProblemLibrary, BVProblemLibrary
 using Sundials, DASKR
@@ -18,19 +19,19 @@ if isdir(ordinartdiffeq_docs_path)
     # Create the OrdinaryDiffEq API directory in the docs
     ordinary_diffeq_dest = joinpath(@__DIR__, "src", "api", "ordinarydiffeq")
     mkpath(dirname(ordinary_diffeq_dest))
-    
+
     # Copy all the docs from OrdinaryDiffEq.jl
-    cp(ordinartdiffeq_docs_path, ordinary_diffeq_dest, force=true)
-    
+    cp(ordinartdiffeq_docs_path, ordinary_diffeq_dest, force = true)
+
     # Copy the pages.jl file from OrdinaryDiffEq.jl
     ordinary_diffeq_pages_dest = joinpath(@__DIR__, "ordinarydiffeq_pages.jl")
     ordinary_diffeq_pages_file = joinpath(ordinartdiffeq_docs_root, "pages.jl")
-    cp(ordinary_diffeq_pages_file, ordinary_diffeq_pages_dest, force=true)
-    
+    cp(ordinary_diffeq_pages_file, ordinary_diffeq_pages_dest, force = true)
+
     # Copy the common_first_steps.jl file from OrdinaryDiffEq.jl
     common_first_steps_dest = joinpath(@__DIR__, "common_first_steps.jl")
     common_first_steps_file = joinpath(ordinartdiffeq_docs_root, "common_first_steps.jl")
-    cp(common_first_steps_file, common_first_steps_dest, force=true)
+    cp(common_first_steps_file, common_first_steps_dest, force = true)
 end
 
 # Copy StochasticDiffEq.jl documentation
@@ -40,14 +41,14 @@ if isdir(stochasticdiffeq_docs_path)
     # Create the StochasticDiffEq API directory in the docs
     stochastic_diffeq_dest = joinpath(@__DIR__, "src", "api", "stochasticdiffeq")
     mkpath(dirname(stochastic_diffeq_dest))
-    
+
     # Copy all the docs from StochasticDiffEq.jl
-    cp(stochasticdiffeq_docs_path, stochastic_diffeq_dest, force=true)
-    
+    cp(stochasticdiffeq_docs_path, stochastic_diffeq_dest, force = true)
+
     # Copy the pages.jl file from StochasticDiffEq.jl
     stochastic_diffeq_pages_dest = joinpath(@__DIR__, "stochasticdiffeq_pages.jl")
     stochastic_diffeq_pages_file = joinpath(stochasticdiffeq_docs_root, "pages.jl")
-    cp(stochastic_diffeq_pages_file, stochastic_diffeq_pages_dest, force=true)
+    cp(stochastic_diffeq_pages_file, stochastic_diffeq_pages_dest, force = true)
 end
 
 ENV["PLOTS_TEST"] = "true"
@@ -118,13 +119,13 @@ makedocs(
         "https://github.com/SciML/ColPrac/blob/master/README.md",
         "https://github.com/SciML/DiffEqDevTools.jl/blob/master/src/ode_tableaus.jl",
         "https://github.com/SciML/DiffEqProblemLibrary.jl/blob/master/lib/BVProblemLibrary/src/BVProblemLibrary.jl",
-        "https://github.com/SciML/DiffEqProblemLibrary.jl/blob/master/lib/DDEProblemLibrary/src/DDEProblemLibrary.jl",
+        "https://github.com/SciML/DiffEqProblemLibrary.jl/blob/master/lib/DDEProblemLibrary/src/DDEProblemLibrary.jl"
     ],
     doctest = false, clean = true,
     warnonly = [:missing_docs, :docs_block],
     format = Documenter.HTML(assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/DiffEqDocs/stable/",
-              size_threshold = 500 * 2^10),
+        size_threshold = 500 * 2^10),
     sitename = "DifferentialEquations.jl",
     authors = "Chris Rackauckas",
     pages = pages)

diff --git a/docs/pages.jl b/docs/pages.jl
@@ -31,7 +31,7 @@ stochastic_diffeq_pages_file = joinpath(@__DIR__, "stochasticdiffeq_pages.jl")
 stochastic_diffeq_pages = []
 if isfile(stochastic_diffeq_pages_file)
     include(stochastic_diffeq_pages_file)
-    
+
     # Transform StochasticDiffEq pages to have the api/stochasticdiffeq prefix
     function transform_stochasticdiffeq_pages(pages_array)
         transformed = []
@@ -49,7 +49,7 @@ if isfile(stochastic_diffeq_pages_file)
         end
         return transformed
     end
-    
+
     stochastic_diffeq_pages = transform_stochasticdiffeq_pages(pages)
 end
 

diff --git a/docs/src/basics/faq.md b/docs/src/basics/faq.md
@@ -12,9 +12,9 @@ For guidelines on debugging ODE solve issues, see
 First of all, don't panic. You may have experienced one of the following warnings:
 
 > dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
->
+> 
 > NaN dt detected. Likely a NaN value in the state, parameters, or derivative value caused this outcome.
->
+> 
 > Instability detected. Aborting
 
 These are all pointing to a similar behavior: for some reason or another, the
@@ -572,7 +572,8 @@ though, an `Error: SingularException` is also possible if the linear solver fail
 
 ```julia
 import DifferentialEquations as DE, OrdinaryDiffEq as ODE, LinearSolve
-DE.solve(prob, ODE.Rodas4(linsolve = LinearSolve.KLUFactorization(; reuse_symbolic = false)))
+DE.solve(
+    prob, ODE.Rodas4(linsolve = LinearSolve.KLUFactorization(; reuse_symbolic = false)))
 ```
 
 For more details about possible linear solvers, consult the [LinearSolve.jl documentation](https://docs.sciml.ai/LinearSolve/stable/)

diff --git a/docs/src/basics/integrator.md b/docs/src/basics/integrator.md
@@ -250,7 +250,8 @@ For example, if one wants to iterate but only stop at specific values, one can
 choose:
 
 ```julia
-integrator = DE.init(prob, DE.Tsit5(); dt = 1 // 2^(4), tstops = [0.5], advance_to_tstop = true)
+integrator = DE.init(
+    prob, DE.Tsit5(); dt = 1 // 2^(4), tstops = [0.5], advance_to_tstop = true)
 for (u, t) in tuples(integrator)
     @test t ∈ [0.5, 1.0]
 end

diff --git a/docs/src/basics/plot.md b/docs/src/basics/plot.md
@@ -160,7 +160,8 @@ xy = Plots.plot(sol, plotdensity = 10000, idxs = (1, 2))
 xz = Plots.plot(sol, plotdensity = 10000, idxs = (1, 3))
 yz = Plots.plot(sol, plotdensity = 10000, idxs = (2, 3))
 xyz = Plots.plot(sol, plotdensity = 10000, idxs = (1, 2, 3))
-Plots.plot(Plots.plot(xyzt, xyz), Plots.plot(xy, xz, yz, layout = (1, 3), w = 1), layout = (2, 1))
+Plots.plot(
+    Plots.plot(xyzt, xyz), Plots.plot(xy, xz, yz, layout = (1, 3), w = 1), layout = (2, 1))
 ```
 
 An example using the functions:

diff --git a/docs/src/basics/solution.md b/docs/src/basics/solution.md
@@ -1,6 +1,6 @@
 # [Solution Handling](@id solution)
 
-The solution is an `RecursiveArrayTools.AbstractDiffEqArray`. 
+The solution is an `RecursiveArrayTools.AbstractDiffEqArray`.
 [See RecursiveArrayTools.jl for more information on the interface](https://docs.sciml.ai/RecursiveArrayTools/stable/).
 The following is a more DiffEq-centric explanation of the interface.
 
@@ -19,7 +19,7 @@ derivative at each timestep `du` or the spatial discretization `x`, `y`, etc.
 ## Array Interface
 
 !!! note
-
+    
     In 2023 the linear indexing `sol[i]` was deprecated. It previously had the behavior that
     `sol[i] = sol.u[i]`. However, this is incompatible with standard `AbstractArray` interfaces,
     Since if `A = VectorOfArray([[1,2],[3,4]])` and `A` is supposed to act like `[1 3; 2 4]`,
@@ -50,7 +50,7 @@ will address first by component and lastly by time, and thus
 sol[i, j]
 ```
 
-will be the `i`th component at timestep `j`. Hence, `sol[j][i] == sol[i, j]`. This is done because Julia is column-major, 
+will be the `i`th component at timestep `j`. Hence, `sol[j][i] == sol[i, j]`. This is done because Julia is column-major,
 so the leading dimension should be contiguous in memory. If the independent variables had shape
 (for example, was a matrix), then `i` is the linear index. We can also access
 solutions with shape:
@@ -186,12 +186,14 @@ error state of the solution. Return codes are now implemented as an enum using E
 rather than symbols.
 
 To check if a solution was successful, use:
+
 ```julia
 SciMLBase.successful_retcode(sol)
 ```
 
 !!! warning
-    Previous iterations of the interface suggested using `sol.retcode == :Success`, 
+
+    Previous iterations of the interface suggested using `sol.retcode == :Success`,
     however, that is now not advised because there are more than one return code that can be interpreted
     as successful. For example, `Terminated` is a successful run to a manual termination, and would be missed
     if only checking for Success. Therefore we highly recommend you use `SciMLBase.successful_retcode(sol)` instead.
@@ -214,7 +216,7 @@ following are major return codes to know:
   - `ConvergenceFailure`: The internal implicit solvers failed to converge.
   - `Failure`: General uncategorized failures or errors.
 
-For a complete list of return codes and their properties, see the 
+For a complete list of return codes and their properties, see the
 [SciMLBase ReturnCode documentation](https://docs.sciml.ai/SciMLBase/stable/interfaces/Solutions/#retcodes).
 
 ## Problem-Specific Features

diff --git a/docs/src/examples/beeler_reuter.md b/docs/src/examples/beeler_reuter.md
@@ -93,8 +93,9 @@ end
 The finite-difference Laplacian is calculated in-place by a 5-point stencil. The Neumann boundary condition is enforced.
 
 !!! note
-    For more complex PDE discretizations, consider using [MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/) 
-    which can automatically generate finite difference discretizations, or [SciMLOperators.jl](https://docs.sciml.ai/SciMLOperators/stable/) 
+
+    For more complex PDE discretizations, consider using [MethodOfLines.jl](https://docs.sciml.ai/MethodOfLines/stable/)
+    which can automatically generate finite difference discretizations, or [SciMLOperators.jl](https://docs.sciml.ai/SciMLOperators/stable/)
     for defining matrix-free linear operators.
 
 ```julia

diff --git a/docs/src/examples/classical_physics.md b/docs/src/examples/classical_physics.md
@@ -88,7 +88,8 @@ sol = ODE.solve(prob, ODE.DPRKN6())
 Plots.plot(sol, idxs = [2, 1], linewidth = 2, title = "Simple Harmonic Oscillator",
     xaxis = "Time", yaxis = "Elongation", label = ["x" "dx"])
 Plots.plot!(t -> A * cos(ω * t - ϕ), lw = 3, ls = :dash, label = "Analytical Solution x")
-Plots.plot!(t -> -A * ω * sin(ω * t - ϕ), lw = 3, ls = :dash, label = "Analytical Solution dx")
+Plots.plot!(
+    t -> -A * ω * sin(ω * t - ϕ), lw = 3, ls = :dash, label = "Analytical Solution dx")
 ```
 
 Note that the order of the variables (and initial conditions) is `dx`, `x`.
@@ -374,15 +375,17 @@ function HH_acceleration!(dv, v, u, p, t)
 end
 initial_positions = [0.0, 0.1]
 initial_velocities = [0.5, 0.0]
-prob = ODE.SecondOrderODEProblem(HH_acceleration!, initial_velocities, initial_positions, tspan)
+prob = ODE.SecondOrderODEProblem(
+    HH_acceleration!, initial_velocities, initial_positions, tspan)
 sol2 = ODE.solve(prob, ODE.KahanLi8(), dt = 1 / 10);
 ```
 
 Notice that we get the same results:
 
 ```@example physics
 # Plot the orbit
-Plots.plot(sol2, idxs = (3, 4), title = "The orbit of the Hénon-Heiles system", xaxis = "x",
+Plots.plot(
+    sol2, idxs = (3, 4), title = "The orbit of the Hénon-Heiles system", xaxis = "x",
     yaxis = "y", leg = false)
 ```
 

diff --git a/docs/src/examples/kepler_problem.md b/docs/src/examples/kepler_problem.md
@@ -33,7 +33,9 @@ sol = ODE.solve(prob, ODE.KahanLi6(), dt = 1 // 10);
 Let's plot the orbit and check the energy and angular momentum variation. We know that energy and angular momentum should be constant, and they are also called first integrals.
 
 ```@example kepler
-plot_orbit(sol) = Plots.plot(sol, idxs = (3, 4), lab = "Orbit", title = "Kepler Problem Solution")
+function plot_orbit(sol)
+    Plots.plot(sol, idxs = (3, 4), lab = "Orbit", title = "Kepler Problem Solution")
+end
 
 function plot_first_integrals(sol, H, L)
     Plots.plot(initial_first_integrals[1] .- map(u -> H(u.x[2], u.x[1]), sol.u),

diff --git a/docs/src/features/callback_functions.md b/docs/src/features/callback_functions.md
@@ -437,7 +437,8 @@ value of `u[1] = -1.2647055847076505e-15`. You can see this by changing the
 `rootfind` argument of the callback:
 
 ```@example callback4
-floor_event = DE.ContinuousCallback(condition, floor_aff!, rootfind = DE.SciMLBase.RightRootFind)
+floor_event = DE.ContinuousCallback(
+    condition, floor_aff!, rootfind = DE.SciMLBase.RightRootFind)
 u0 = [1.0, 0.0]
 p = [1.0]
 prob = DE.ODEProblem{true}(dynamics!, u0, (0.0, 1.75), p)

diff --git a/docs/src/getting_started.md b/docs/src/getting_started.md
@@ -404,7 +404,8 @@ M = t -> 0.1sin(t)                    # external torque [Nm]
 prob = DE.ODEProblem(pendulum!, u₀, tspan, M)
 sol = DE.solve(prob)
 
-Plots.plot(sol, linewidth = 2, xaxis = "t", label = ["θ [rad]" "ω [rad/s]"], layout = (2, 1))
+Plots.plot(
+    sol, linewidth = 2, xaxis = "t", label = ["θ [rad]" "ω [rad/s]"], layout = (2, 1))
 ```
 
 Note how the external **time-varying** torque `M` is introduced as a **parameter** in the `pendulum!` function. Indeed, as a general principle the parameters can be any type; here we specify `M` as time-varying by representing it by a function, which is expressed by appending the dependence on time `(t)` to the name of the parameter.
@@ -458,9 +459,9 @@ above, with the only change being the type for the initial condition and constan
 ```@example ODE4
 import StaticArrays
 A = StaticArrays.@SMatrix [1.0 0.0 0.0 -5.0
-              4.0 -2.0 4.0 -3.0
-              -4.0 0.0 0.0 1.0
-              5.0 -2.0 2.0 3.0]
+                           4.0 -2.0 4.0 -3.0
+                           -4.0 0.0 0.0 1.0
+                           5.0 -2.0 2.0 3.0]
 u0 = StaticArrays.@SMatrix rand(4, 2)
 tspan = (0.0, 1.0)
 f2(u, p, t) = A * u

diff --git a/docs/src/index.md b/docs/src/index.md
@@ -67,7 +67,7 @@ Additionally, DifferentialEquations.jl comes with built-in analysis features, in
     for guidance on PRs, issues, and other matters relating to contributing to SciML.
   - See the [SciML Style Guide](https://github.com/SciML/SciMLStyle) for common coding practices and other style decisions.
   - There are a few community forums:
-
+    
       + The #diffeq-bridged and #sciml-bridged channels in the
         [Julia Slack](https://julialang.org/slack/)
       + The #diffeq-bridged and #sciml-bridged channels in the

diff --git a/docs/src/solvers/ode_solve.md b/docs/src/solvers/ode_solve.md
@@ -960,22 +960,22 @@ import IRKGaussLegendre
 `IRKGL16(;kwargs...)` has the following arguments:
 
   - second_order_ode (boolean):
-
+    
       + =false (default): for a ODEProblem type
       + =true: for a second-order differential equation
 
   - simd (boolean):
-
+    
       + =true: SIMD-vectorized implementation only available for Float32 or Float64 computations
       + =false (default):  generic implementation that can use with arbitrary Julia-defined number systems
   - mstep: output saved at every 'mstep' steps. Default 1.
   - initial_extrapolation: initialization method for stages.
-
+    
       + =false: simplest initialization
       + =true (default): extrapolation from the stage values of previous step
   - maxtrials: maximum number of attempts to accept adaptive step size
   - threading
-
+    
       + =false (default): sequential execution of the numerical integration
       + =true: computations using threads (shared memory multi-threading) for stage-wise parallelization
 

diff --git a/docs/src/tutorials/advanced_ode_example.md b/docs/src/tutorials/advanced_ode_example.md
@@ -186,7 +186,8 @@ Now let's see how the version with sparsity compares to the version without:
 import BenchmarkTools as BT # for @btime
 BT.@btime DE.solve(prob_ode_brusselator_2d, DE.TRBDF2(); save_everystep = false);
 BT.@btime DE.solve(prob_ode_brusselator_2d_sparse, DE.TRBDF2(); save_everystep = false);
-BT.@btime DE.solve(prob_ode_brusselator_2d_sparse, DE.KenCarp47(; linsolve = DE.KLUFactorization());
+BT.@btime DE.solve(
+    prob_ode_brusselator_2d_sparse, DE.KenCarp47(; linsolve = DE.KLUFactorization());
     save_everystep = false);
 nothing # hide
 ```
@@ -273,7 +274,8 @@ which is more automatic. The setup is very similar to before:
 import AlgebraicMultigrid
 function algebraicmultigrid(W, du, u, p, t, newW, Plprev, Prprev, solverdata)
     if newW === nothing || newW
-        Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(convert(AbstractMatrix, W)))
+        Pl = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(convert(
+            AbstractMatrix, W)))
     else
         Pl = Plprev
     end
@@ -332,7 +334,8 @@ Newton Krylov (with numerical differentiation). Thus, on this problem we could d
 import Sundials
 BT.@btime DE.solve(prob_ode_brusselator_2d, Sundials.CVODE_BDF(); save_everystep = false);
 # Simplest speedup: use :LapackDense
-BT.@btime DE.solve(prob_ode_brusselator_2d, Sundials.CVODE_BDF(; linear_solver = :LapackDense);
+BT.@btime DE.solve(
+    prob_ode_brusselator_2d, Sundials.CVODE_BDF(; linear_solver = :LapackDense);
     save_everystep = false);
 # GMRES Version: Doesn't require any extra stuff!
 BT.@btime DE.solve(prob_ode_brusselator_2d, Sundials.CVODE_BDF(; linear_solver = :GMRES);
@@ -428,7 +431,8 @@ function psetupamg(p, t, u, du, jok, jcurPtr, gamma)
         @. @view(W[idxs]) = @view(W[idxs]) + 1
 
         # Build preconditioner on W
-        preccache2[] = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(W,
+        preccache2[] = AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(
+            W,
             presmoother = AlgebraicMultigrid.Jacobi(rand(size(W,
                 1))),
             postsmoother = AlgebraicMultigrid.Jacobi(rand(size(W,

diff --git a/docs/src/tutorials/bvp_example.md b/docs/src/tutorials/bvp_example.md
@@ -132,7 +132,9 @@ end
 u0 = [1.0, 1.0, 1.0]
 tspan = (0.0, 1.0)
 prob = BVP.SecondOrderBVProblem(f!, bc!, u0, tspan)
-sol = BVP.solve(prob, BVP.MIRKN4(; jac_alg = BVP.BVPJacobianAlgorithm(BVP.AutoForwardDiff())); dt = 0.01)
+sol = BVP.solve(
+    prob, BVP.MIRKN4(; jac_alg = BVP.BVPJacobianAlgorithm(BVP.AutoForwardDiff()));
+    dt = 0.01)
 ```
 
 ## Example 3: Semi-Explicit Boundary Value Differential-Algebraic Equations
@@ -174,6 +176,7 @@ tspan = (0.0, 1.0)
 fun = BVP.BVPFunction(f!, bc!, mass_matrix = [1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 0])
 prob = BVP.BVProblem(fun, u0, tspan)
 sol = BVP.solve(prob,
-    BVP.Ascher4(; zeta = [0.0, 0.0, 1.0], jac_alg = BVP.BVPJacobianAlgorithm(BVP.AutoForwardDiff()));
+    BVP.Ascher4(;
+        zeta = [0.0, 0.0, 1.0], jac_alg = BVP.BVPJacobianAlgorithm(BVP.AutoForwardDiff()));
     dt = 0.01)
 ```
diff --git a/docs/src/tutorials/faster_ode_example.md b/docs/src/tutorials/faster_ode_example.md
@@ -350,7 +350,8 @@ function rober_static(u, p, t)
     du3 = k₂ * y₂^2
     StaticArrays.SA[du1, du2, du3]
 end
-prob = DE.ODEProblem(rober_static, StaticArrays.SA[1.0, 0.0, 0.0], (0.0, 1e5), StaticArrays.SA[0.04, 3e7, 1e4])
+prob = DE.ODEProblem(rober_static, StaticArrays.SA[1.0, 0.0, 0.0],
+    (0.0, 1e5), StaticArrays.SA[0.04, 3e7, 1e4])
 sol = DE.solve(prob, DE.Rosenbrock23())
 ```
 
@@ -717,13 +718,15 @@ nothing # hide
 ```
 
 ```@example faster_ode3
-BT.@btime DE.solve(prob, Sundials.CVODE_BDF(; linear_solver = :GMRES); save_everystep = false);
+BT.@btime DE.solve(
+    prob, Sundials.CVODE_BDF(; linear_solver = :GMRES); save_everystep = false);
 nothing # hide
 ```
 
 ```@example faster_ode3
 prob = DE.ODEProblem(fast_gm!, r0, (0.0, 500.0), p)
-BT.@btime DE.solve(prob, Sundials.CVODE_BDF(; linear_solver = :GMRES); save_everystep = false);
+BT.@btime DE.solve(
+    prob, Sundials.CVODE_BDF(; linear_solver = :GMRES); save_everystep = false);
 nothing # hide
 ```
 

diff --git a/docs/src/tutorials/sde_example.md b/docs/src/tutorials/sde_example.md
@@ -164,7 +164,7 @@ In general, a system of SDEs
 du = f(u,p,t)dt + g(u,p,t)dW,
 ```
 
-where `u` is now a vector of variables, `f` is a vector, and `g` is a matrix, is numerically integrated in the same way as ODEs. A common scenario, which is the default for DifferentialEquations.jl, is that every variable in the system gets a different random kick. This is the case when `g` is a diagonal matrix. Correspondingly, we say that we have a diagonal noise. 
+where `u` is now a vector of variables, `f` is a vector, and `g` is a matrix, is numerically integrated in the same way as ODEs. A common scenario, which is the default for DifferentialEquations.jl, is that every variable in the system gets a different random kick. This is the case when `g` is a diagonal matrix. Correspondingly, we say that we have a diagonal noise.
 
 We handle this in a simple manner by defining the deterministic part `f!(du,u,p,t)` and the stochastic part
 `g!(du2,u,p,t)` as in-place functions, but note that our convention is that the function `g!` only defines and modifies the diagonal entries of `g` matrix.