ChrisRackauckas · Aug 6, 2026 · Aug 6, 2026 · Aug 6, 2026
diff --git a/Project.toml b/Project.toml
@@ -11,6 +11,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
@@ -21,6 +22,7 @@ LinearAlgebra = "1"
 PreallocationTools = "1.2"
 PrecompileTools = "1.2.1"
 SciMLBase = "3.1"
+SciMLOperators = "1.22.0"
 SparseArrays = "1"
 Test = "1"
 julia = "1"

diff --git a/docs/make.jl b/docs/make.jl
@@ -207,7 +207,6 @@ makedocs(;
     ),
     draft = IS_LIVESERVER,
     pages = _PAGES,
-    warnonly = true
 )
 
 deploydocs(;

diff --git a/docs/src/literate_wyos/diffusion_equations.jl b/docs/src/literate_wyos/diffusion_equations.jl
@@ -320,7 +320,7 @@ fvm_prob = FVMProblem(
     final_time
 )
 
-using LinearSolve #src
+using LinearSolve
 
 # ````julia
 # using BenchmarkTools

diff --git a/docs/src/math.md b/docs/src/math.md
@@ -133,8 +133,8 @@ fig, ax, sc = poly(control_volume_verts[cv_idx],
 triplot!(tri,
     show_convex_hull=false,
     show_constrained_edges=false,
-    show_all_points=true,
-    point_color=:black,
+    show_points=true,
+    markercolor=:black,
     markersize=11)
 scatter!(ax, [a, b, c], color=:cyan, markersize=11, strokecolor=:black, strokewidth=1)
 linesegments!(ax, linesegments, color=:blue, linewidth=3)

diff --git a/docs/src/tutorials/diffusion_equation_in_a_wedge_with_mixed_boundary_conditions.md b/docs/src/tutorials/diffusion_equation_in_a_wedge_with_mixed_boundary_conditions.md
@@ -109,7 +109,7 @@ In my experience, I've found that `TRBDF2(linsolve=KLUFactorization())` typicall
 has the best performance for these problems.
 
 ````@example diffusion_equation_in_a_wedge_with_mixed_boundary_conditions
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.01, parallel = Val(false))
 ind = findall(DelaunayTriangulation.each_point_index(tri)) do i #hide
     !DelaunayTriangulation.has_vertex(tri, i) #hide
@@ -172,7 +172,7 @@ prob = FVMProblem(mesh, BCs; diffusion_function = D, initial_condition, final_ti
 
 flux = (x, y, t, α, β, γ, p) -> (-α, -β)
 
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.01, parallel = Val(false))
 
 using CairoMakie

diff --git a/docs/src/tutorials/diffusion_equation_on_an_annulus.md b/docs/src/tutorials/diffusion_equation_on_an_annulus.md
@@ -99,7 +99,7 @@ prob = FVMProblem(mesh, BCs;
 ````
 
 ````@example diffusion_equation_on_an_annulus
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.2)
 sol |> tc #hide
 ````
@@ -246,7 +246,7 @@ prob = FVMProblem(mesh, BCs;
     final_time,
     initial_condition)
 
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.2)
 
 fig = Figure(fontsize = 38)

diff --git a/...ls/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.md b/...ls/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.md
@@ -78,7 +78,7 @@ prob = FVMSystem(u_prob, v_prob)
 Now that we have our system, we can solve.
 
 ````@example gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 10.0, parallel = Val(false))
 sol |> tc #hide
 ````

diff --git a/docs/src/tutorials/helmholtz_equation_with_inhomogeneous_boundary_conditions.md b/docs/src/tutorials/helmholtz_equation_with_inhomogeneous_boundary_conditions.md
@@ -91,7 +91,7 @@ algorithm above.
 using NonlinearSolve
 sol = solve(steady_prob, NewtonRaphson())
 copyto!(prob.initial_condition, sol.u) # this also changes steady_prob's initial condition
-using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
+using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization())))
 sol |> tc #hide
 ````
@@ -132,7 +132,7 @@ steady_prob = SteadyFVMProblem(prob)
 using NonlinearSolve
 sol = solve(steady_prob, NewtonRaphson())
 copyto!(prob.initial_condition, sol.u) # this also changes steady_prob's initial condition
-using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
+using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization())))
 
 using CairoMakie

diff --git a/docs/src/tutorials/laplaces_equation_with_internal_dirichlet_conditions.md b/docs/src/tutorials/laplaces_equation_with_internal_dirichlet_conditions.md
@@ -150,7 +150,7 @@ steady_prob = SteadyFVMProblem(prob)
 Now let's solve the problem.
 
 ````@example laplaces_equation_with_internal_dirichlet_conditions
-using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
+using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization())))
 sol |> tc #hide
 ````
@@ -226,7 +226,7 @@ prob = FVMProblem(mesh, BCs, ICs;
 
 steady_prob = SteadyFVMProblem(prob)
 
-using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
+using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization())))
 
 fig, ax, sc = tricontourf(tri, sol.u, levels = LinRange(0, 100, 28))

diff --git a/...near_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.md b/...near_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.md
@@ -126,7 +126,7 @@ prob = FVMProblem(mesh, BCs;
 Now we can solve and visualise the solution.
 
 ````@example piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 times = [0, 10, 25, 50, 100, 200, 250]
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = times)
 sol |> tc #hide
@@ -351,7 +351,7 @@ prob = FVMProblem(mesh, BCs;
     flux_parameters,
     final_time)
 
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 times = [0, 10, 25, 50, 100, 200, 250]
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = times)
 

diff --git a/docs/src/tutorials/porous_fisher_equation_and_travelling_waves.md b/docs/src/tutorials/porous_fisher_equation_and_travelling_waves.md
@@ -64,7 +64,7 @@ takes the solution some time to evolve towards the travelling wave solution.
 Now with this preamble out of the way, let us solve this problem.
 
 ````@example porous_fisher_equation_and_travelling_waves
-using DelaunayTriangulation, FiniteVolumeMethod, OrdinaryDiffEq, LinearSolve
+using DelaunayTriangulation, FiniteVolumeMethod, OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 a, b, c, d, nx, ny = 0.0, 3.0, 0.0, 40.0, 60, 80
 tri = triangulate_rectangle(a, b, c, d, nx, ny; single_boundary = false)
 mesh = FVMGeometry(tri)
@@ -144,7 +144,7 @@ An uncommented version of this example is given below.
 You can view the source code for this file [here](https://github.com/SciML/FiniteVolumeMethod.jl/tree/main/docs/src/literate_tutorials/porous_fisher_equation_and_travelling_waves.jl).
 
 ```julia
-using DelaunayTriangulation, FiniteVolumeMethod, OrdinaryDiffEq, LinearSolve
+using DelaunayTriangulation, FiniteVolumeMethod, OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 a, b, c, d, nx, ny = 0.0, 3.0, 0.0, 40.0, 60, 80
 tri = triangulate_rectangle(a, b, c, d, nx, ny; single_boundary = false)
 mesh = FVMGeometry(tri)

diff --git a/docs/src/tutorials/porous_medium_equation.md b/docs/src/tutorials/porous_medium_equation.md
@@ -78,7 +78,7 @@ prob = FVMProblem(mesh, BCs;
 
 ````@example porous_medium_equation
 # Step 4: Solve
-using LinearSolve, OrdinaryDiffEq
+using LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()); saveat = 3.0)
 sol |> tc #hide
 ````
@@ -219,7 +219,7 @@ prob = FVMProblem(mesh, BCs;
     final_time)
 
 # Step 4: Solve
-using LinearSolve, OrdinaryDiffEq
+using LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()); saveat = 3.0)
 
 # Step 5: Visualise

diff --git a/docs/src/tutorials/reaction_diffusion_brusselator_system_of_pdes.md b/docs/src/tutorials/reaction_diffusion_brusselator_system_of_pdes.md
@@ -183,7 +183,7 @@ system = FVMSystem(Φ_prob, Ψ_prob)
 We can now solve the problem just as we've done previously.
 
 ````@example reaction_diffusion_brusselator_system_of_pdes
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(system, TRBDF2(linsolve = KLUFactorization()), saveat = 1.0)
 sol |> tc #hide
 ````
@@ -266,7 +266,7 @@ mesh = FVMGeometry(tri)
 
 system = FVMSystem(Φ_prob, Ψ_prob)
 
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(system, TRBDF2(linsolve = KLUFactorization()), saveat = 1.0)
 
 sol.u[3]

diff --git a/docs/src/tutorials/solving_mazes_with_laplaces_equation.md b/docs/src/tutorials/solving_mazes_with_laplaces_equation.md
@@ -90,7 +90,7 @@ steady_prob = SteadyFVMProblem(prob)
 ````
 
 ````@example solving_mazes_with_laplaces_equation
-using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
+using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization(), autodiff = AutoFiniteDiff())))
 sol |> tc #hide
 ````
@@ -200,7 +200,7 @@ prob = FVMProblem(mesh, BCs;
     final_time = final_time)
 steady_prob = SteadyFVMProblem(prob)
 
-using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
+using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization(), autodiff = AutoFiniteDiff())))
 
 tricontourf(tri, sol.u, colormap = :matter)

diff --git a/docs/src/wyos/diffusion_equations.md b/docs/src/wyos/diffusion_equations.md
@@ -277,7 +277,7 @@ end
 Let's now test the function. We use the same problem as in [this tutorial](../tutorials/diffusion_equation_on_a_square_plate.md).
 
 ````@example diffusion_equations
-using DelaunayTriangulation, OrdinaryDiffEq, LinearAlgebra, SparseArrays
+using DelaunayTriangulation, OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearAlgebra, SparseArrays, SciMLOperators
 tri = triangulate_rectangle(0, 2, 0, 2, 50, 50, single_boundary = true)
 mesh = FVMGeometry(tri)
 BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> zero(x), Dirichlet)
@@ -351,6 +351,10 @@ fvm_prob = FVMProblem(mesh, BCs;
     final_time)
 ````
 
+````@example diffusion_equations
+using LinearSolve
+````
+
 ````julia
 using BenchmarkTools
 @btime solve($diff_eq, $Tsit5(), saveat = $0.05);
@@ -590,7 +594,7 @@ function diffusion_equation(mesh::FVMGeometry,
     return prob
 end
 
-using DelaunayTriangulation, OrdinaryDiffEq, LinearAlgebra, SparseArrays
+using DelaunayTriangulation, OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearAlgebra, SparseArrays, SciMLOperators
 tri = triangulate_rectangle(0, 2, 0, 2, 50, 50, single_boundary = true)
 mesh = FVMGeometry(tri)
 BCs = BoundaryConditions(mesh, (x, y, t, u, p) -> zero(x), Dirichlet)

diff --git a/docs/src/wyos/laplaces_equation.md b/docs/src/wyos/laplaces_equation.md
@@ -148,7 +148,7 @@ fvm_prob = SteadyFVMProblem(FVMProblem(mesh, BCs;
 ````
 
 ````@example laplaces_equation
-using SteadyStateDiffEq, OrdinaryDiffEq
+using SteadyStateDiffEq, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 fvm_sol = solve(fvm_prob, DynamicSS(TRBDF2()))
 fvm_sol |> tc #hide
 ````
@@ -315,7 +315,7 @@ fvm_prob = SteadyFVMProblem(FVMProblem(mesh, BCs;
     initial_condition,
     final_time = Inf))
 
-using SteadyStateDiffEq, OrdinaryDiffEq
+using SteadyStateDiffEq, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 fvm_sol = solve(fvm_prob, DynamicSS(TRBDF2()))
 
 ax = Axis(fig[1, 2], xlabel = "x", ylabel = "y", width = 600, height = 600)

diff --git a/docs/src/wyos/linear_reaction_diffusion_equations.md b/docs/src/wyos/linear_reaction_diffusion_equations.md
@@ -38,7 +38,7 @@ We can also reuse `apply_dirichlet_conditions!`, `apply_dudt_conditions`, and
 `boundary_edge_contributions!` from the diffusion equation example. Here is our implementation.
 
 ````@example linear_reaction_diffusion_equations
-using FiniteVolumeMethod, SparseArrays, OrdinaryDiffEq, LinearAlgebra
+using FiniteVolumeMethod, SparseArrays, OrdinaryDiffEq, LinearAlgebra, SciMLOperators
 const FVM = FiniteVolumeMethod
 function linear_source_contributions!(
         A, mesh, conditions, source_function, source_parameters)
@@ -191,7 +191,7 @@ An uncommented version of this example is given below.
 You can view the source code for this file [here](https://github.com/SciML/FiniteVolumeMethod.jl/tree/main/docs/src/literate_wyos/linear_reaction_diffusion_equations.jl).
 
 ```julia
-using FiniteVolumeMethod, SparseArrays, OrdinaryDiffEq, LinearAlgebra
+using FiniteVolumeMethod, SparseArrays, OrdinaryDiffEq, LinearAlgebra, SciMLOperators
 const FVM = FiniteVolumeMethod
 function linear_source_contributions!(
         A, mesh, conditions, source_function, source_parameters)

diff --git a/docs/src/wyos/mean_exit_time.md b/docs/src/wyos/mean_exit_time.md
@@ -203,7 +203,7 @@ fvm_prob = SteadyFVMProblem(FVMProblem(mesh, BCs, ICs;
 Let's compare the two solutions.
 
 ````@example mean_exit_time
-using SteadyStateDiffEq, OrdinaryDiffEq
+using SteadyStateDiffEq, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 fvm_sol = solve(fvm_prob, DynamicSS(TRBDF2()))
 fvm_sol |> tc #hide
 ````
@@ -364,7 +364,7 @@ fvm_prob = SteadyFVMProblem(FVMProblem(mesh, BCs, ICs;
     final_time = Inf,
     initial_condition))
 
-using SteadyStateDiffEq, OrdinaryDiffEq
+using SteadyStateDiffEq, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 fvm_sol = solve(fvm_prob, DynamicSS(TRBDF2()))
 
 ax = Axis(fig[1, 2], width = 600, height = 600, title = "Template")

diff --git a/docs/src/wyos/poissons_equation.md b/docs/src/wyos/poissons_equation.md
@@ -151,7 +151,7 @@ fvm_prob = SteadyFVMProblem(FVMProblem(mesh, BCs;
 ````
 
 ````@example poissons_equation
-using SteadyStateDiffEq, OrdinaryDiffEq
+using SteadyStateDiffEq, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 fvm_sol = solve(fvm_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization())))
 fvm_sol |> tc #hide
 ````
@@ -491,7 +491,7 @@ fvm_prob = SteadyFVMProblem(FVMProblem(mesh, BCs;
     initial_condition,
     final_time = Inf))
 
-using SteadyStateDiffEq, OrdinaryDiffEq
+using SteadyStateDiffEq, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 fvm_sol = solve(fvm_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization())))
 
 prob = PoissonsEquation(mesh, BCs; source_function = source_function)

diff --git a/src/FiniteVolumeMethod.jl b/src/FiniteVolumeMethod.jl
@@ -16,7 +16,8 @@ using DelaunayTriangulation: DelaunayTriangulation, Triangulation,
 using LinearAlgebra: LinearAlgebra, norm
 using PreallocationTools: PreallocationTools, DiffCache, get_tmp
 using SciMLBase: SciMLBase, CallbackSet, DiscreteCallback, LinearProblem,
-    MatrixOperator, ODEFunction, ODEProblem, SteadyStateProblem
+    ODEFunction, ODEProblem, SteadyStateProblem
+using SciMLOperators: SciMLOperators, MatrixOperator
 using SparseArrays: SparseArrays, sparse
 using Base.Threads