diff --git a/Project.toml b/Project.toml
index e4140ff2..6619fb36 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "FiniteVolumeMethod"
 uuid = "d4f04ab7-4f65-4d72-8a28-7087bc7f46f4"
 authors = ["Daniel VandenHeuvel <danj.vandenheuvel@gmail.com>"]
-version = "1.2.1"
+version = "1.2.2"
 
 [deps]
 ChunkSplitters = "ae650224-84b6-46f8-82ea-d812ca08434e"
@@ -14,12 +14,14 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-ChunkSplitters = "0.1, 1.0, 2.0, 3.2"
-CommonSolve = "0.2"
-DelaunayTriangulation = "1.0"
-PreallocationTools = "0.4, 1.2"
-PrecompileTools = "1.2"
-SciMLBase = "2.34, 3.1"
+ChunkSplitters = "3.2"
+CommonSolve = "0.2.4"
+DelaunayTriangulation = "1.6.6"
+LinearAlgebra = "1"
+PreallocationTools = "1.2"
+PrecompileTools = "1.2.1"
+SciMLBase = "3.1"
+SparseArrays = "1"
 Test = "1"
 julia = "1"
 
diff --git a/README.md b/README.md
index f17a304c..4b059392 100644
--- a/README.md
+++ b/README.md
@@ -42,7 +42,7 @@ sol = solve(prob, Tsit5(), saveat = 0.001)
 u = Observable(sol.u[1])
 fig, ax, sc = tricontourf(tri, u, levels = 0:5:50, colormap = :matter)
 tightlimits!(ax)
-record(fig, "anim.gif", eachindex(sol)) do i
+record(fig, "anim.gif", eachindex(sol.u)) do i
     u[] = sol.u[i]
 end
 ```
diff --git a/docs/Project.toml b/docs/Project.toml
index c115f5db..7952a116 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -25,11 +25,13 @@ LiveServer = "16fef848-5104-11e9-1b77-fb7a48bbb589"
 NaturalNeighbours = "f16ad982-4edb-46b1-8125-78e5a8b5a9e6"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
 PProf = "e4faabce-9ead-11e9-39d9-4379958e3056"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SimpleGraphs = "55797a34-41de-5266-9ec1-32ac4eb504d3"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
@@ -64,11 +66,13 @@ LiveServer = "1.5.0"
 NaturalNeighbours = "1.3.6"
 NonlinearSolve = "4.19.1"
 OrdinaryDiffEq = "7.0.0"
+OrdinaryDiffEqSDIRK = "2.7.0"
 PProf = "3.2.0"
 ReferenceTests = "0.10.6"
 Roots = "3.0.0"
 SciMLBase = "3.18.0"
-SimpleGraphs = "0.8.10"
+SciMLOperators = "1.22.0"
+SimpleGraphs = "0.8.8 - 0.8"
 Sparspak = "0.3.15"
 SpecialFunctions = "2.8.0"
 StableRNGs = "1.0.4"
diff --git a/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.png b/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.png
index 38aefd7e..dbebcf91 100644
Binary files a/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.png and b/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.png differ
diff --git a/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_exact_comparisons.png b/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_exact_comparisons.png
index 54741dfc..41d6e72f 100644
Binary files a/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_exact_comparisons.png and b/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_exact_comparisons.png differ
diff --git a/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_natural_neighbour_interpolation.png b/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_natural_neighbour_interpolation.png
index a5da960f..f2c095bd 100644
Binary files a/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_natural_neighbour_interpolation.png and b/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_natural_neighbour_interpolation.png differ
diff --git a/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_piecewise_linear_interpolation.png b/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_piecewise_linear_interpolation.png
index 339ba640..e5b1230c 100644
Binary files a/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_piecewise_linear_interpolation.png and b/docs/src/figures/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation_piecewise_linear_interpolation.png differ
diff --git a/docs/src/literate_tutorials/diffusion_equation_in_a_wedge_with_mixed_boundary_conditions.jl b/docs/src/literate_tutorials/diffusion_equation_in_a_wedge_with_mixed_boundary_conditions.jl
index b4a0e5f6..54906c5e 100644
--- a/docs/src/literate_tutorials/diffusion_equation_in_a_wedge_with_mixed_boundary_conditions.jl
+++ b/docs/src/literate_tutorials/diffusion_equation_in_a_wedge_with_mixed_boundary_conditions.jl
@@ -79,13 +79,13 @@ flux = (x, y, t, α, β, γ, p) -> (-α, -β)
 # We now solve the problem. We provide the solver for this problem.
 # In my experience, I've found that `TRBDF2(linsolve=KLUFactorization())` typically
 # has the best performance for these problems.
-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
 end #hide
 using Test #hide
-@test sol[ind, :] ≈ reshape(repeat(initial_condition, length(sol)), :, length(sol))[ind, :] # make sure that missing vertices don't change #hide
+@test sol[ind, :] ≈ reshape(repeat(initial_condition, length(sol.u)), :, length(sol.u))[ind, :] # make sure that missing vertices don't change #hide
 sol |> tc #hide
 
 #-
@@ -151,12 +151,12 @@ function exact_solution(x, y, t, A, ζ, f, α) #src
 end #src
 function compare_solutions(sol, tri, α, f) #src
     n = DelaunayTriangulation.num_points(tri) #src
-    x = zeros(n, length(sol)) #src
-    y = zeros(n, length(sol)) #src
-    u = zeros(n, length(sol)) #src
+    x = zeros(n, length(sol.u)) #src
+    y = zeros(n, length(sol.u)) #src
+    u = zeros(n, length(sol.u)) #src
     ζ = get_ζ_terms(20, 20, α) #src
     A = get_sum_coefficients(20, 20, α, ζ) #src
-    for i in eachindex(sol) #src
+    for i in eachindex(sol.u) #src
         for j in each_solid_vertex(tri) #src
             x[j, i], y[j, i] = get_point(tri, j) #src
             u[j, i] = exact_solution(x[j, i], y[j, i], sol.t[i], A, ζ, f, α) #src
@@ -166,7 +166,7 @@ function compare_solutions(sol, tri, α, f) #src
 end #src
 x, y, u = compare_solutions(sol, tri, α, f) #src
 fig = Figure(fontsize = 64) #src
-for i in eachindex(sol) #src
+for i in eachindex(sol.u) #src
     local ax #src
     ax = Axis(fig[1, i], width = 600, height = 600) #src
     tricontourf!(ax, tri, sol.u[i], levels = 0:0.01:1, colormap = :matter) #src
diff --git a/docs/src/literate_tutorials/diffusion_equation_on_a_square_plate.jl b/docs/src/literate_tutorials/diffusion_equation_on_a_square_plate.jl
index 2916ce2e..3bf9ed4d 100644
--- a/docs/src/literate_tutorials/diffusion_equation_on_a_square_plate.jl
+++ b/docs/src/literate_tutorials/diffusion_equation_on_a_square_plate.jl
@@ -94,10 +94,10 @@ function exact_solution(x, y, t) #src
 end #src
 function compare_solutions(sol, tri) #src
     n = DelaunayTriangulation.num_points(tri) #src
-    x = zeros(n, length(sol)) #src
-    y = zeros(n, length(sol)) #src
-    u = zeros(n, length(sol)) #src
-    for i in eachindex(sol) #src
+    x = zeros(n, length(sol.u)) #src
+    y = zeros(n, length(sol.u)) #src
+    u = zeros(n, length(sol.u)) #src
+    for i in eachindex(sol.u) #src
         for j in each_solid_vertex(tri) #src
             x[j, i], y[j, i] = get_point(tri, j) #src
             u[j, i] = exact_solution(x[j, i], y[j, i], sol.t[i]) #src
@@ -107,7 +107,7 @@ function compare_solutions(sol, tri) #src
 end #src
 x, y, u = compare_solutions(sol, tri) #src
 fig = Figure(fontsize = 64) #src
-for i in eachindex(sol) #src
+for i in eachindex(sol.u) #src
     local ax #src
     ax = Axis(fig[1, i], width = 600, height = 600) #src
     tricontourf!(ax, tri, sol.u[i], levels = 0:5:50, colormap = :matter) #src
diff --git a/docs/src/literate_tutorials/diffusion_equation_on_an_annulus.jl b/docs/src/literate_tutorials/diffusion_equation_on_an_annulus.jl
index b23d6184..236519aa 100644
--- a/docs/src/literate_tutorials/diffusion_equation_on_an_annulus.jl
+++ b/docs/src/literate_tutorials/diffusion_equation_on_an_annulus.jl
@@ -79,7 +79,7 @@ prob = FVMProblem(
 )
 
 #-
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.2)
 sol |> tc #hide
 
diff --git a/docs/src/literate_tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.jl b/docs/src/literate_tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.jl
index 74b20aaf..9e198af5 100644
--- a/docs/src/literate_tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.jl
+++ b/docs/src/literate_tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.jl
@@ -63,7 +63,7 @@ v_prob = FVMProblem(
 prob = FVMSystem(u_prob, v_prob)
 
 # Now that we have our system, we can solve.
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 10.0, parallel = Val(false))
 sol |> tc #hide
 
@@ -79,7 +79,7 @@ y = LinRange(-1, 1, 200)
 heatmap!(ax, x, y, u, colorrange = (0.0, 0.4))
 hidedecorations!(ax)
 record(
-    fig, joinpath(@__DIR__, "../figures", "gray_scott_patterns.mp4"), eachindex(sol);
+    fig, joinpath(@__DIR__, "../figures", "gray_scott_patterns.mp4"), eachindex(sol.u);
     framerate = 60
 ) do _i
     i[] = _i
diff --git a/docs/src/literate_tutorials/helmholtz_equation_with_inhomogeneous_boundary_conditions.jl b/docs/src/literate_tutorials/helmholtz_equation_with_inhomogeneous_boundary_conditions.jl
index df2600d7..fef54023 100644
--- a/docs/src/literate_tutorials/helmholtz_equation_with_inhomogeneous_boundary_conditions.jl
+++ b/docs/src/literate_tutorials/helmholtz_equation_with_inhomogeneous_boundary_conditions.jl
@@ -64,7 +64,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())))
 sol |> tc #hide
 
diff --git a/docs/src/literate_tutorials/laplaces_equation_with_internal_dirichlet_conditions.jl b/docs/src/literate_tutorials/laplaces_equation_with_internal_dirichlet_conditions.jl
index 9c161b38..82e7eecd 100644
--- a/docs/src/literate_tutorials/laplaces_equation_with_internal_dirichlet_conditions.jl
+++ b/docs/src/literate_tutorials/laplaces_equation_with_internal_dirichlet_conditions.jl
@@ -116,7 +116,7 @@ prob = FVMProblem(
 steady_prob = SteadyFVMProblem(prob)
 
 # Now let's solve the problem.
-using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
+using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization())))
 sol |> tc #hide
 
diff --git a/docs/src/literate_tutorials/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.jl b/docs/src/literate_tutorials/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.jl
index 10516745..59862294 100644
--- a/docs/src/literate_tutorials/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.jl
+++ b/docs/src/literate_tutorials/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.jl
@@ -25,7 +25,7 @@ tc = DisplayAs.withcontext(:displaysize => (15, 80), :limit => true); #hide
 # with \eqref{eq:advdiffeq}. For the mesh, we could use
 # `triangulate_rectangle`, but we want to put most of the triangles
 # near the origin, so we need to use `refine!` on an initial mesh.
-using DelaunayTriangulation, FiniteVolumeMethod, LinearAlgebra, CairoMakie
+using DelaunayTriangulation, FiniteVolumeMethod, LinearAlgebra, CairoMakie, StableRNGs
 L = 30
 tri = triangulate_rectangle(-L, L, -L, L, 2, 2, single_boundary = true)
 tot_area = get_area(tri)
@@ -39,7 +39,7 @@ area_constraint = (_tri, T) -> begin
     flag = A ≥ max_area_function(A, dist_to_origin)
     return flag
 end
-refine!(tri; min_angle = 33.0, custom_constraint = area_constraint)
+refine!(tri; min_angle = 33.0, custom_constraint = area_constraint, rng = StableRNG(123))
 triplot(tri)
 
 #-
@@ -96,7 +96,7 @@ prob = FVMProblem(
 )
 
 # Now we can solve and visualise the solution.
-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
@@ -105,7 +105,7 @@ sol |> tc #hide
 using CairoMakie
 using ReferenceTests #src
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax = Axis(
         fig[1, i], width = 400, height = 400,
         xlabel = "x", ylabel = "y",
@@ -133,9 +133,9 @@ function exact_solution(x, y, t, D, ν) #src
     return 1 / (4 * D * π * t) * exp(1 / (4 * D * t) * (-(x - ν * t)^2 - y^2)) #src
 end #src
 function get_errs(_sol, tri, flux_parameters) #src
-    _errs = zeros(length(_sol)) #src
+    _errs = zeros(length(_sol.u)) #src
     _err = zeros(DelaunayTriangulation.num_points(tri)) #src
-    for i in eachindex(_sol) #src
+    for i in eachindex(_sol.u) #src
         !DelaunayTriangulation.has_vertex(tri, i) && continue #src
         i == 1 && continue #src
         m = maximum(_sol.u[i]) #src
@@ -205,8 +205,8 @@ for j in eachindex(y)
         triangles[i, j] = jump_and_march(tri, (x[i], y[j]))
     end
 end
-interpolated_vals = zeros(length(x), length(y), length(sol))
-for k in eachindex(sol)
+interpolated_vals = zeros(length(x), length(y), length(sol.u))
+for k in eachindex(sol.u)
     for j in eachindex(y)
         for i in eachindex(x)
             interpolated_vals[i, j, k] = pl_interpolate(
@@ -218,9 +218,9 @@ end
 
 # Let's visualise these results to check their accuracy. We compute the triangulation of
 # our grid to make the `tricontourf` call faster.
-_tri = triangulate([[x for x in x, _ in y] |> vec [y for _ in x, y in y] |> vec]')
+_tri = triangulate([[x for x in x, _ in y] |> vec [y for _ in x, y in y] |> vec]', rng = StableRNG(123))
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax = Axis(
         fig[1, i], width = 400, height = 400,
         xlabel = "x", ylabel = "y",
diff --git a/docs/src/literate_tutorials/porous_fisher_equation_and_travelling_waves.jl b/docs/src/literate_tutorials/porous_fisher_equation_and_travelling_waves.jl
index 86960058..21cf2f37 100644
--- a/docs/src/literate_tutorials/porous_fisher_equation_and_travelling_waves.jl
+++ b/docs/src/literate_tutorials/porous_fisher_equation_and_travelling_waves.jl
@@ -47,7 +47,7 @@ tc = DisplayAs.withcontext(:displaysize => (15, 80), :limit => true); #hide
 # 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.
-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)
@@ -81,10 +81,10 @@ c = sqrt(λ / (2D))
 cₘᵢₙ = sqrt(λ * D / 2)
 zᶜ = 0.0
 exact_solution(z) = ifelse(z ≤ zᶜ, 1 - exp(cₘᵢₙ * z), zero(z))
-travelling_wave_values = zeros(ny, length(sol) - large_time_idx + 1)
+travelling_wave_values = zeros(ny, length(sol.u) - large_time_idx + 1)
 z_vals = zero(travelling_wave_values)
 u_mat = [reshape(u, (nx, ny)) for u in sol.u]
-for (i, t_idx) in pairs(large_time_idx:lastindex(sol))
+for (i, t_idx) in pairs(large_time_idx:lastindex(sol.u))
     u = u_mat[t_idx]
     τ = sol.t[t_idx]
     for k in 1:ny
@@ -109,10 +109,10 @@ for (i, j) in zip(1:3, (1, 51, 101))
     tightlimits!(ax)
 end
 ax = Axis(fig[1, 4], width = 900, height = 600)
-colors = cgrad(:matter, length(sol) - large_time_idx + 1; categorical = false)
+colors = cgrad(:matter, length(sol.u) - large_time_idx + 1; categorical = false)
 [
     lines!(ax, z_vals[:, i], travelling_wave_values[:, i], color = colors[i], linewidth = 2)
-        for i in 1:(length(sol) - large_time_idx + 1)
+        for i in 1:(length(sol.u) - large_time_idx + 1)
 ]
 exact_z_vals = collect(LinRange(extrema(z_vals)..., 500))
 exact_travelling_wave_values = exact_solution.(exact_z_vals)
diff --git a/docs/src/literate_tutorials/porous_medium_equation.jl b/docs/src/literate_tutorials/porous_medium_equation.jl
index 54d6708c..e94966eb 100644
--- a/docs/src/literate_tutorials/porous_medium_equation.jl
+++ b/docs/src/literate_tutorials/porous_medium_equation.jl
@@ -60,7 +60,7 @@ prob = FVMProblem(
 
 #-
 ## Step 4: Solve
-using LinearSolve, OrdinaryDiffEq
+using LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()); saveat = 3.0)
 sol |> tc #hide
 
@@ -99,10 +99,10 @@ function exact_solution(x, y, t, m, M, D) #src
 end #src
 function compare_solutions(sol, tri, m, M, D) #src
     n = DelaunayTriangulation.num_points(tri) #src
-    x = zeros(n, length(sol)) #src
-    y = zeros(n, length(sol)) #src
-    u = zeros(n, length(sol)) #src
-    for i in eachindex(sol) #src
+    x = zeros(n, length(sol.u)) #src
+    y = zeros(n, length(sol.u)) #src
+    u = zeros(n, length(sol.u)) #src
+    for i in eachindex(sol.u) #src
         !DelaunayTriangulation.has_vertex(tri, i) && continue #src
         for j in each_solid_vertex(tri) #src
             x[j, i], y[j, i] = get_point(tri, j) #src
@@ -113,7 +113,7 @@ function compare_solutions(sol, tri, m, M, D) #src
 end #src
 x, y, u = compare_solutions(sol, tri, m, M, D) #src
 fig = Figure(fontsize = 64) #src
-for i in eachindex(sol) #src
+for i in eachindex(sol.u) #src
     ax = Axis(fig[1, i], width = 600, height = 600) #src
     tricontourf!(ax, tri, sol.u[i], levels = 0:0.005:0.1, colormap = :matter) #src
     ax = Axis(fig[2, i], width = 600, height = 600) #src
@@ -207,10 +207,10 @@ function exact_solution(x, y, t, m, M, D, λ) #src
 end #src
 function compare_solutions(sol, tri, m, M, D, λ) #src
     n = DelaunayTriangulation.num_solid_vertices(tri) #src
-    x = zeros(n, length(sol)) #src
-    y = zeros(n, length(sol)) #src
-    u = zeros(n, length(sol)) #src
-    for i in eachindex(sol) #src
+    x = zeros(n, length(sol.u)) #src
+    y = zeros(n, length(sol.u)) #src
+    u = zeros(n, length(sol.u)) #src
+    for i in eachindex(sol.u) #src
         for j in each_solid_vertex(tri) #src
             x[j, i], y[j, i] = get_point(tri, j) #src
             u[j, i] = exact_solution(x[j, i], y[j, i], sol.t[i], m, M, D, λ) #src
@@ -220,7 +220,7 @@ function compare_solutions(sol, tri, m, M, D, λ) #src
 end #src
 x, y, u = compare_solutions(sol, tri, m, M, D, λ) #src
 fig = Figure(fontsize = 64) #src
-for i in eachindex(sol) #src
+for i in eachindex(sol.u) #src
     ax = Axis(fig[1, i], width = 600, height = 600) #src
     tricontourf!(
         ax, tri, sol.u[i], levels = 0:0.05:1, extendlow = :auto,
diff --git a/docs/src/literate_tutorials/reaction_diffusion_brusselator_system_of_pdes.jl b/docs/src/literate_tutorials/reaction_diffusion_brusselator_system_of_pdes.jl
index cd8c1b9f..082854aa 100644
--- a/docs/src/literate_tutorials/reaction_diffusion_brusselator_system_of_pdes.jl
+++ b/docs/src/literate_tutorials/reaction_diffusion_brusselator_system_of_pdes.jl
@@ -143,7 +143,7 @@ mesh = FVMGeometry(tri)
 system = FVMSystem(Φ_prob, Ψ_prob)
 
 # We can now solve the problem just as we've done previously.
-using OrdinaryDiffEq, LinearSolve
+using OrdinaryDiffEq, OrdinaryDiffEqSDIRK, LinearSolve
 sol = solve(system, TRBDF2(linsolve = KLUFactorization()), saveat = 1.0)
 sol |> tc #hide
 
@@ -159,7 +159,7 @@ sol.u[3][1, :] |> tc #hide
 # are the values of $\Psi$ at the third time. We can visualise the solutions as follows:
 using CairoMakie
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax1 = Axis(
         fig[1, i], xlabel = L"x", ylabel = L"y",
         width = 400, height = 400,
@@ -170,8 +170,8 @@ for i in eachindex(sol)
         width = 400, height = 400,
         title = L"\Psi: t = %$(sol.t[i])", titlealign = :left
     )
-    tricontourf!(ax1, tri, sol[i][1, :], levels = 0:0.1:1, colormap = :matter)
-    tricontourf!(ax2, tri, sol[i][2, :], levels = 1:10:100, colormap = :matter)
+    tricontourf!(ax1, tri, sol.u[i][1, :], levels = 0:0.1:1, colormap = :matter)
+    tricontourf!(ax2, tri, sol.u[i][2, :], levels = 1:10:100, colormap = :matter)
 end
 resize_to_layout!(fig)
 fig
@@ -180,7 +180,7 @@ using ReferenceTests #src
 
 x = getx.(get_points(tri)) #src
 y = gety.(get_points(tri)) #src
-for i in eachindex(sol) #src
+for i in eachindex(sol.u) #src
     ax3 = Axis(
         fig[3, i], xlabel = L"x", ylabel = L"y", #src
         width = 400, height = 400, #src
diff --git a/docs/src/literate_tutorials/reaction_diffusion_equation_with_a_time_dependent_dirichlet_boundary_condition_on_a_disk.jl b/docs/src/literate_tutorials/reaction_diffusion_equation_with_a_time_dependent_dirichlet_boundary_condition_on_a_disk.jl
index 14b9d65f..6c28a6e6 100644
--- a/docs/src/literate_tutorials/reaction_diffusion_equation_with_a_time_dependent_dirichlet_boundary_condition_on_a_disk.jl
+++ b/docs/src/literate_tutorials/reaction_diffusion_equation_with_a_time_dependent_dirichlet_boundary_condition_on_a_disk.jl
@@ -53,7 +53,7 @@ prob = FVMProblem(
 
 # We can now solve.
 using OrdinaryDiffEq, LinearSolve
-alg = FBDF(linsolve = UMFPACKFactorization(), autodiff = false)
+alg = FBDF(linsolve = UMFPACKFactorization(), autodiff = AutoFiniteDiff())
 sol = solve(prob, alg, saveat = 0.01)
 sol |> tc #hide
 
@@ -85,10 +85,10 @@ function exact_solution(x, y, t) #src
 end #src
 function compare_solutions(sol, tri) #src
     n = DelaunayTriangulation.num_points(tri) #src
-    x = zeros(n, length(sol)) #src
-    y = zeros(n, length(sol)) #src
-    u = zeros(n, length(sol)) #src
-    for i in eachindex(sol) #src
+    x = zeros(n, length(sol.u)) #src
+    y = zeros(n, length(sol.u)) #src
+    u = zeros(n, length(sol.u)) #src
+    for i in eachindex(sol.u) #src
         !DelaunayTriangulation.has_vertex(tri, i) && continue #src
         for j in each_solid_vertex(tri) #src
             x[j, i], y[j, i] = get_point(tri, j) #src
@@ -99,7 +99,7 @@ function compare_solutions(sol, tri) #src
 end #src
 x, y, u = compare_solutions(sol, tri) #src
 fig = Figure(fontsize = 64) #src
-for i in eachindex(sol) #src
+for i in eachindex(sol.u) #src
     ax = Axis(fig[1, i], width = 600, height = 600) #src
     tricontourf!(ax, tri, sol.u[i], levels = 1:0.01:1.4, colormap = :matter) #src
     ax = Axis(fig[2, i], width = 600, height = 600) #src
diff --git a/docs/src/literate_tutorials/solving_mazes_with_laplaces_equation.jl b/docs/src/literate_tutorials/solving_mazes_with_laplaces_equation.jl
index 535987d8..8bfb791d 100644
--- a/docs/src/literate_tutorials/solving_mazes_with_laplaces_equation.jl
+++ b/docs/src/literate_tutorials/solving_mazes_with_laplaces_equation.jl
@@ -77,8 +77,8 @@ prob = FVMProblem(
 steady_prob = SteadyFVMProblem(prob)
 
 #-
-using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
-sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization(), autodiff = false)))
+using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
+sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization(), autodiff = AutoFiniteDiff())))
 sol |> tc #hide
 
 # We now have our solution.
@@ -122,7 +122,7 @@ fig, ax,
 hidedecorations!(ax)
 tightlimits!(ax)
 record(
-    fig, joinpath(@__DIR__, "../figures", "maze_solution_1.mp4"), eachindex(sol);
+    fig, joinpath(@__DIR__, "../figures", "maze_solution_1.mp4"), eachindex(sol.u);
     framerate = 24
 ) do _i
     i[] = _i
diff --git a/docs/src/literate_wyos/diffusion_equations.jl b/docs/src/literate_wyos/diffusion_equations.jl
index e607e3e3..50830562 100644
--- a/docs/src/literate_wyos/diffusion_equations.jl
+++ b/docs/src/literate_wyos/diffusion_equations.jl
@@ -246,7 +246,7 @@ function diffusion_equation(
 end
 
 # Let's now test the function. We use the same problem as in [this tutorial](../tutorials/diffusion_equation_on_a_square_plate.md).
-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)
@@ -383,7 +383,7 @@ sol |> tc #hide
 
 #-
 fig = Figure(fontsize = 38)
-for j in eachindex(sol)
+for j in eachindex(sol.u)
     ax = Axis(
         fig[1, j], width = 600, height = 600,
         xlabel = "x", ylabel = "y",
@@ -416,7 +416,7 @@ fvm_prob = FVMProblem(
 fvm_sol = solve(fvm_prob, TRBDF2(linsolve = KLUFactorization()); saveat = 100.0)
 fvm_sol |> tc #hide
 
-for j in eachindex(fvm_sol)
+for j in eachindex(fvm_sol.u)
     ax = Axis(
         fig[2, j], width = 600, height = 600,
         xlabel = "x", ylabel = "y",
diff --git a/docs/src/literate_wyos/laplaces_equation.jl b/docs/src/literate_wyos/laplaces_equation.jl
index 6ed1c00f..00728da2 100644
--- a/docs/src/literate_wyos/laplaces_equation.jl
+++ b/docs/src/literate_wyos/laplaces_equation.jl
@@ -125,7 +125,7 @@ fvm_prob = SteadyFVMProblem(
 )
 
 #-
-using SteadyStateDiffEq, OrdinaryDiffEq
+using SteadyStateDiffEq, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 fvm_sol = solve(fvm_prob, DynamicSS(TRBDF2()))
 fvm_sol |> tc #hide
 
diff --git a/docs/src/literate_wyos/linear_reaction_diffusion_equations.jl b/docs/src/literate_wyos/linear_reaction_diffusion_equations.jl
index 92e0f5fa..a375972f 100644
--- a/docs/src/literate_wyos/linear_reaction_diffusion_equations.jl
+++ b/docs/src/literate_wyos/linear_reaction_diffusion_equations.jl
@@ -22,7 +22,7 @@ tc = DisplayAs.withcontext(:displaysize => (15, 80), :limit => true); #hide
 # as in the previous sections, but we will need an extra function to add $f(\vb x)$ to the appropriate diagonals.
 # We can also reuse `apply_dirichlet_conditions!`, `apply_dudt_conditions`, and
 # `boundary_edge_contributions!` from the diffusion equation example. Here is our implementation.
-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
@@ -101,7 +101,7 @@ sol |> tc #hide
 #-
 using CairoMakie
 fig = Figure(fontsize = 38)
-for j in eachindex(sol)
+for j in eachindex(sol.u)
     ax = Axis(
         fig[1, j], width = 600, height = 600,
         xlabel = "x", ylabel = "y",
@@ -140,7 +140,7 @@ fvm_prob = FVMProblem(
 fvm_sol = solve(fvm_prob, Tsit5(), saveat = 2.0)
 fvm_sol |> tc #hide
 
-for j in eachindex(fvm_sol) #src
+for j in eachindex(fvm_sol.u) #src
     ax = Axis(
         fig[2, j], width = 600, height = 600, #src
         xlabel = "x", ylabel = "y", #src
diff --git a/docs/src/literate_wyos/mean_exit_time.jl b/docs/src/literate_wyos/mean_exit_time.jl
index a95ad68b..6c599b83 100644
--- a/docs/src/literate_wyos/mean_exit_time.jl
+++ b/docs/src/literate_wyos/mean_exit_time.jl
@@ -177,7 +177,7 @@ fvm_prob = SteadyFVMProblem(
 )
 
 # Let's compare the two solutions.
-using SteadyStateDiffEq, OrdinaryDiffEq
+using SteadyStateDiffEq, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 fvm_sol = solve(fvm_prob, DynamicSS(TRBDF2()))
 fvm_sol |> tc #hide
 
diff --git a/docs/src/literate_wyos/poissons_equation.jl b/docs/src/literate_wyos/poissons_equation.jl
index ed3c5dcd..ecb9df15 100644
--- a/docs/src/literate_wyos/poissons_equation.jl
+++ b/docs/src/literate_wyos/poissons_equation.jl
@@ -133,7 +133,7 @@ fvm_prob = SteadyFVMProblem(
 )
 
 #-
-using SteadyStateDiffEq, OrdinaryDiffEq
+using SteadyStateDiffEq, OrdinaryDiffEq, OrdinaryDiffEqSDIRK
 fvm_sol = solve(fvm_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization())))
 fvm_sol |> tc #hide
 using ReferenceTests #src
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
index a00b355f..d02d9103 100644
--- 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
@@ -115,7 +115,7 @@ ind = findall(DelaunayTriangulation.each_point_index(tri)) do i #hide
     !DelaunayTriangulation.has_vertex(tri, i) #hide
 end #hide
 using Test #hide
-@test sol[ind, :] ≈ reshape(repeat(initial_condition, length(sol)), :, length(sol))[ind, :] # make sure that missing vertices don't change #hide
+@test sol[ind, :] ≈ reshape(repeat(initial_condition, length(sol.u)), :, length(sol.u))[ind, :] # make sure that missing vertices don't change #hide
 sol |> tc #hide
 ````
 
diff --git a/docs/src/tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.md b/docs/src/tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.md
index ec9e9dcb..3a3c2e63 100644
--- a/docs/src/tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.md
+++ b/docs/src/tutorials/gray_scott_model_turing_patterns_from_a_coupled_reaction_diffusion_system.md
@@ -96,7 +96,7 @@ x = LinRange(-1, 1, 200)
 y = LinRange(-1, 1, 200)
 heatmap!(ax, x, y, u, colorrange = (0.0, 0.4))
 hidedecorations!(ax)
-record(fig, joinpath(@__DIR__, "../figures", "gray_scott_patterns.mp4"), eachindex(sol);
+record(fig, joinpath(@__DIR__, "../figures", "gray_scott_patterns.mp4"), eachindex(sol.u);
     framerate = 60) do _i
     i[] = _i
 end
@@ -157,7 +157,7 @@ x = LinRange(-1, 1, 200)
 y = LinRange(-1, 1, 200)
 heatmap!(ax, x, y, u, colorrange = (0.0, 0.4))
 hidedecorations!(ax)
-record(fig, joinpath(@__DIR__, "../figures", "gray_scott_patterns.mp4"), eachindex(sol);
+record(fig, joinpath(@__DIR__, "../figures", "gray_scott_patterns.mp4"), eachindex(sol.u);
     framerate = 60) do _i
     i[] = _i
 end
diff --git a/docs/src/tutorials/keller_segel_chemotaxis.md b/docs/src/tutorials/keller_segel_chemotaxis.md
index 0aebbec7..37f4a0aa 100644
--- a/docs/src/tutorials/keller_segel_chemotaxis.md
+++ b/docs/src/tutorials/keller_segel_chemotaxis.md
@@ -78,7 +78,7 @@ heatmap!(axu, x, y, u, colorrange = (0.0, 2.5), colormap = :turbo)
 heatmap!(axv, x, y, v, colorrange = (0.0, 10.0), colormap = :turbo)
 resize_to_layout!(fig)
 record(
-    fig, joinpath(@__DIR__, "../figures", "keller_segel_chemotaxis.mp4"), eachindex(sol);
+    fig, joinpath(@__DIR__, "../figures", "keller_segel_chemotaxis.mp4"), eachindex(sol.u);
     framerate = 60) do _i
     i[] = _i
 end;
@@ -151,7 +151,7 @@ heatmap!(axu, x, y, u, colorrange = (0.0, 2.5), colormap = :turbo)
 heatmap!(axv, x, y, v, colorrange = (0.0, 10.0), colormap = :turbo)
 resize_to_layout!(fig)
 record(
-    fig, joinpath(@__DIR__, "../figures", "keller_segel_chemotaxis.mp4"), eachindex(sol);
+    fig, joinpath(@__DIR__, "../figures", "keller_segel_chemotaxis.mp4"), eachindex(sol.u);
     framerate = 60) do _i
     i[] = _i
 end;
diff --git a/docs/src/tutorials/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.md b/docs/src/tutorials/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.md
index ae39a594..d7dac821 100644
--- a/docs/src/tutorials/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.md
+++ b/docs/src/tutorials/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.md
@@ -39,7 +39,7 @@ with \eqref{eq:advdiffeq}. For the mesh, we could use
 near the origin, so we need to use `refine!` on an initial mesh.
 
 ````@example piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation
-using DelaunayTriangulation, FiniteVolumeMethod, LinearAlgebra, CairoMakie
+using DelaunayTriangulation, FiniteVolumeMethod, LinearAlgebra, CairoMakie, StableRNGs
 L = 30
 tri = triangulate_rectangle(-L, L, -L, L, 2, 2, single_boundary = true)
 tot_area = get_area(tri)
@@ -53,7 +53,7 @@ area_constraint = (_tri, T) -> begin
     flag = A ≥ max_area_function(A, dist_to_origin)
     return flag
 end
-refine!(tri; min_angle = 33.0, custom_constraint = area_constraint)
+refine!(tri; min_angle = 33.0, custom_constraint = area_constraint, rng = StableRNG(123))
 triplot(tri)
 ````
 
@@ -135,7 +135,7 @@ sol |> tc #hide
 ````@example piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation
 using CairoMakie
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax = Axis(fig[1, i], width = 400, height = 400,
         xlabel = "x", ylabel = "y",
         title = "t = $(sol.t[i])",
@@ -177,8 +177,8 @@ for j in eachindex(y)
         triangles[i, j] = jump_and_march(tri, (x[i], y[j]))
     end
 end
-interpolated_vals = zeros(length(x), length(y), length(sol))
-for k in eachindex(sol)
+interpolated_vals = zeros(length(x), length(y), length(sol.u))
+for k in eachindex(sol.u)
     for j in eachindex(y)
         for i in eachindex(x)
             interpolated_vals[i, j, k] = pl_interpolate(
@@ -192,9 +192,9 @@ Let's visualise these results to check their accuracy. We compute the triangulat
 our grid to make the `tricontourf` call faster.
 
 ````@example piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation
-_tri = triangulate([[x for x in x, _ in y] |> vec [y for _ in x, y in y] |> vec]')
+_tri = triangulate([[x for x in x, _ in y] |> vec [y for _ in x, y in y] |> vec]', rng = StableRNG(123))
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax = Axis(fig[1, i], width = 400, height = 400,
         xlabel = "x", ylabel = "y",
         title = "t = $(sol.t[i])",
@@ -311,7 +311,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/piecewise_linear_and_natural_neighbour_interpolation_for_an_advection_diffusion_equation.jl).
 
 ```julia
-using DelaunayTriangulation, FiniteVolumeMethod, LinearAlgebra, CairoMakie
+using DelaunayTriangulation, FiniteVolumeMethod, LinearAlgebra, CairoMakie, StableRNGs
 L = 30
 tri = triangulate_rectangle(-L, L, -L, L, 2, 2, single_boundary = true)
 tot_area = get_area(tri)
@@ -325,7 +325,7 @@ area_constraint = (_tri, T) -> begin
     flag = A ≥ max_area_function(A, dist_to_origin)
     return flag
 end
-refine!(tri; min_angle = 33.0, custom_constraint = area_constraint)
+refine!(tri; min_angle = 33.0, custom_constraint = area_constraint, rng = StableRNG(123))
 triplot(tri)
 
 mesh = FVMGeometry(tri)
@@ -357,7 +357,7 @@ sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = times)
 
 using CairoMakie
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax = Axis(fig[1, i], width = 400, height = 400,
         xlabel = "x", ylabel = "y",
         title = "t = $(sol.t[i])",
@@ -378,8 +378,8 @@ for j in eachindex(y)
         triangles[i, j] = jump_and_march(tri, (x[i], y[j]))
     end
 end
-interpolated_vals = zeros(length(x), length(y), length(sol))
-for k in eachindex(sol)
+interpolated_vals = zeros(length(x), length(y), length(sol.u))
+for k in eachindex(sol.u)
     for j in eachindex(y)
         for i in eachindex(x)
             interpolated_vals[i, j, k] = pl_interpolate(
@@ -388,9 +388,9 @@ for k in eachindex(sol)
     end
 end
 
-_tri = triangulate([[x for x in x, _ in y] |> vec [y for _ in x, y in y] |> vec]')
+_tri = triangulate([[x for x in x, _ in y] |> vec [y for _ in x, y in y] |> vec]', rng = StableRNG(123))
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax = Axis(fig[1, i], width = 400, height = 400,
         xlabel = "x", ylabel = "y",
         title = "t = $(sol.t[i])",
diff --git a/docs/src/tutorials/porous_fisher_equation_and_travelling_waves.md b/docs/src/tutorials/porous_fisher_equation_and_travelling_waves.md
index 9079dd88..f1492c6c 100644
--- a/docs/src/tutorials/porous_fisher_equation_and_travelling_waves.md
+++ b/docs/src/tutorials/porous_fisher_equation_and_travelling_waves.md
@@ -99,10 +99,10 @@ c = sqrt(λ / (2D))
 cₘᵢₙ = sqrt(λ * D / 2)
 zᶜ = 0.0
 exact_solution(z) = ifelse(z ≤ zᶜ, 1 - exp(cₘᵢₙ * z), zero(z))
-travelling_wave_values = zeros(ny, length(sol) - large_time_idx + 1)
+travelling_wave_values = zeros(ny, length(sol.u) - large_time_idx + 1)
 z_vals = zero(travelling_wave_values)
 u_mat = [reshape(u, (nx, ny)) for u in sol.u]
-for (i, t_idx) in pairs(large_time_idx:lastindex(sol))
+for (i, t_idx) in pairs(large_time_idx:lastindex(sol.u))
     u = u_mat[t_idx]
     τ = sol.t[t_idx]
     for k in 1:ny
@@ -128,9 +128,9 @@ for (i, j) in zip(1:3, (1, 51, 101))
     tightlimits!(ax)
 end
 ax = Axis(fig[1, 4], width = 900, height = 600)
-colors = cgrad(:matter, length(sol) - large_time_idx + 1; categorical = false)
+colors = cgrad(:matter, length(sol.u) - large_time_idx + 1; categorical = false)
 [lines!(ax, z_vals[:, i], travelling_wave_values[:, i], color = colors[i], linewidth = 2)
- for i in 1:(length(sol) - large_time_idx + 1)]
+ for i in 1:(length(sol.u) - large_time_idx + 1)]
 exact_z_vals = collect(LinRange(extrema(z_vals)..., 500))
 exact_travelling_wave_values = exact_solution.(exact_z_vals)
 lines!(ax, exact_z_vals, exact_travelling_wave_values,
@@ -172,10 +172,10 @@ c = sqrt(λ / (2D))
 cₘᵢₙ = sqrt(λ * D / 2)
 zᶜ = 0.0
 exact_solution(z) = ifelse(z ≤ zᶜ, 1 - exp(cₘᵢₙ * z), zero(z))
-travelling_wave_values = zeros(ny, length(sol) - large_time_idx + 1)
+travelling_wave_values = zeros(ny, length(sol.u) - large_time_idx + 1)
 z_vals = zero(travelling_wave_values)
 u_mat = [reshape(u, (nx, ny)) for u in sol.u]
-for (i, t_idx) in pairs(large_time_idx:lastindex(sol))
+for (i, t_idx) in pairs(large_time_idx:lastindex(sol.u))
     u = u_mat[t_idx]
     τ = sol.t[t_idx]
     for k in 1:ny
@@ -197,9 +197,9 @@ for (i, j) in zip(1:3, (1, 51, 101))
     tightlimits!(ax)
 end
 ax = Axis(fig[1, 4], width = 900, height = 600)
-colors = cgrad(:matter, length(sol) - large_time_idx + 1; categorical = false)
+colors = cgrad(:matter, length(sol.u) - large_time_idx + 1; categorical = false)
 [lines!(ax, z_vals[:, i], travelling_wave_values[:, i], color = colors[i], linewidth = 2)
- for i in 1:(length(sol) - large_time_idx + 1)]
+ for i in 1:(length(sol.u) - large_time_idx + 1)]
 exact_z_vals = collect(LinRange(extrema(z_vals)..., 500))
 exact_travelling_wave_values = exact_solution.(exact_z_vals)
 lines!(ax, exact_z_vals, exact_travelling_wave_values,
diff --git a/docs/src/tutorials/reaction_diffusion_brusselator_system_of_pdes.md b/docs/src/tutorials/reaction_diffusion_brusselator_system_of_pdes.md
index 045a8d35..7e7bb658 100644
--- a/docs/src/tutorials/reaction_diffusion_brusselator_system_of_pdes.md
+++ b/docs/src/tutorials/reaction_diffusion_brusselator_system_of_pdes.md
@@ -208,15 +208,15 @@ are the values of $\Psi$ at the third time. We can visualise the solutions as fo
 ````@example reaction_diffusion_brusselator_system_of_pdes
 using CairoMakie
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax1 = Axis(fig[1, i], xlabel = L"x", ylabel = L"y",
         width = 400, height = 400,
         title = L"\Phi: t = %$(sol.t[i])", titlealign = :left)
     ax2 = Axis(fig[2, i], xlabel = L"x", ylabel = L"y",
         width = 400, height = 400,
         title = L"\Psi: t = %$(sol.t[i])", titlealign = :left)
-    tricontourf!(ax1, tri, sol[i][1, :], levels = 0:0.1:1, colormap = :matter)
-    tricontourf!(ax2, tri, sol[i][2, :], levels = 1:10:100, colormap = :matter)
+    tricontourf!(ax1, tri, sol.u[i][1, :], levels = 0:0.1:1, colormap = :matter)
+    tricontourf!(ax2, tri, sol.u[i][2, :], levels = 1:10:100, colormap = :matter)
 end
 resize_to_layout!(fig)
 fig
@@ -275,15 +275,15 @@ sol.u[3][1, :]
 
 using CairoMakie
 fig = Figure(fontsize = 38)
-for i in eachindex(sol)
+for i in eachindex(sol.u)
     ax1 = Axis(fig[1, i], xlabel = L"x", ylabel = L"y",
         width = 400, height = 400,
         title = L"\Phi: t = %$(sol.t[i])", titlealign = :left)
     ax2 = Axis(fig[2, i], xlabel = L"x", ylabel = L"y",
         width = 400, height = 400,
         title = L"\Psi: t = %$(sol.t[i])", titlealign = :left)
-    tricontourf!(ax1, tri, sol[i][1, :], levels = 0:0.1:1, colormap = :matter)
-    tricontourf!(ax2, tri, sol[i][2, :], levels = 1:10:100, colormap = :matter)
+    tricontourf!(ax1, tri, sol.u[i][1, :], levels = 0:0.1:1, colormap = :matter)
+    tricontourf!(ax2, tri, sol.u[i][2, :], levels = 1:10:100, colormap = :matter)
 end
 resize_to_layout!(fig)
 fig
diff --git a/docs/src/tutorials/reaction_diffusion_equation_with_a_time_dependent_dirichlet_boundary_condition_on_a_disk.md b/docs/src/tutorials/reaction_diffusion_equation_with_a_time_dependent_dirichlet_boundary_condition_on_a_disk.md
index 61d4193d..6d99df08 100644
--- a/docs/src/tutorials/reaction_diffusion_equation_with_a_time_dependent_dirichlet_boundary_condition_on_a_disk.md
+++ b/docs/src/tutorials/reaction_diffusion_equation_with_a_time_dependent_dirichlet_boundary_condition_on_a_disk.md
@@ -74,7 +74,7 @@ We can now solve.
 
 ````@example reaction_diffusion_equation_with_a_time_dependent_dirichlet_boundary_condition_on_a_disk
 using OrdinaryDiffEq, LinearSolve
-alg = FBDF(linsolve = UMFPACKFactorization(), autodiff = false)
+alg = FBDF(linsolve = UMFPACKFactorization(), autodiff = AutoFiniteDiff())
 sol = solve(prob, alg, saveat = 0.01)
 sol |> tc #hide
 ````
@@ -127,7 +127,7 @@ prob = FVMProblem(mesh, BCs;
     initial_condition)
 
 using OrdinaryDiffEq, LinearSolve
-alg = FBDF(linsolve = UMFPACKFactorization(), autodiff = false)
+alg = FBDF(linsolve = UMFPACKFactorization(), autodiff = AutoFiniteDiff())
 sol = solve(prob, alg, saveat = 0.01)
 
 fig = Figure(fontsize = 38)
diff --git a/docs/src/tutorials/solving_mazes_with_laplaces_equation.md b/docs/src/tutorials/solving_mazes_with_laplaces_equation.md
index 5640e40d..ab5383be 100644
--- a/docs/src/tutorials/solving_mazes_with_laplaces_equation.md
+++ b/docs/src/tutorials/solving_mazes_with_laplaces_equation.md
@@ -91,7 +91,7 @@ steady_prob = SteadyFVMProblem(prob)
 
 ````@example solving_mazes_with_laplaces_equation
 using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
-sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization(), autodiff = false)))
+sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization(), autodiff = AutoFiniteDiff())))
 sol |> tc #hide
 ````
 
@@ -137,7 +137,7 @@ sc = tricontourf(tri, ∇norms, colormap = :matter, levels = LinRange(0, 0.0035,
     extendlow = :auto, extendhigh = :auto)
 hidedecorations!(ax)
 tightlimits!(ax)
-record(fig, joinpath(@__DIR__, "../figures", "maze_solution_1.mp4"), eachindex(sol);
+record(fig, joinpath(@__DIR__, "../figures", "maze_solution_1.mp4"), eachindex(sol.u);
     framerate = 24) do _i
     i[] = _i
 end;
@@ -201,7 +201,7 @@ prob = FVMProblem(mesh, BCs;
 steady_prob = SteadyFVMProblem(prob)
 
 using SteadyStateDiffEq, LinearSolve, OrdinaryDiffEq
-sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization(), autodiff = false)))
+sol = solve(steady_prob, DynamicSS(TRBDF2(linsolve = KLUFactorization(), autodiff = AutoFiniteDiff())))
 
 tricontourf(tri, sol.u, colormap = :matter)
 
@@ -228,7 +228,7 @@ sc = tricontourf(tri, ∇norms, colormap = :matter, levels = LinRange(0, 0.0035,
     extendlow = :auto, extendhigh = :auto)
 hidedecorations!(ax)
 tightlimits!(ax)
-record(fig, joinpath(@__DIR__, "../figures", "maze_solution_1.mp4"), eachindex(sol);
+record(fig, joinpath(@__DIR__, "../figures", "maze_solution_1.mp4"), eachindex(sol.u);
     framerate = 24) do _i
     i[] = _i
 end;
diff --git a/docs/src/wyos/diffusion_equations.md b/docs/src/wyos/diffusion_equations.md
index 6308c679..f3761401 100644
--- a/docs/src/wyos/diffusion_equations.md
+++ b/docs/src/wyos/diffusion_equations.md
@@ -414,7 +414,7 @@ sol |> tc #hide
 
 ````@example diffusion_equations
 fig = Figure(fontsize = 38)
-for j in eachindex(sol)
+for j in eachindex(sol.u)
     ax = Axis(fig[1, j], width = 600, height = 600,
         xlabel = "x", ylabel = "y",
         title = "t = $(sol.t[j])",
@@ -444,7 +444,7 @@ fvm_prob = FVMProblem(mesh, BCs_prob;
 fvm_sol = solve(fvm_prob, TRBDF2(linsolve = KLUFactorization()); saveat = 100.0)
 fvm_sol |> tc #hide
 
-for j in eachindex(fvm_sol)
+for j in eachindex(fvm_sol.u)
     ax = Axis(fig[2, j], width = 600, height = 600,
         xlabel = "x", ylabel = "y",
         title = "t = $(fvm_sol.t[j])",
@@ -657,7 +657,7 @@ prob = DiffusionEquation(mesh, BCs;
 sol = solve(prob, Tsit5(); saveat = 100.0)
 
 fig = Figure(fontsize = 38)
-for j in eachindex(sol)
+for j in eachindex(sol.u)
     ax = Axis(fig[1, j], width = 600, height = 600,
         xlabel = "x", ylabel = "y",
         title = "t = $(sol.t[j])",
@@ -679,7 +679,7 @@ fvm_prob = FVMProblem(mesh, BCs_prob;
     final_time)
 fvm_sol = solve(fvm_prob, TRBDF2(linsolve = KLUFactorization()); saveat = 100.0)
 
-for j in eachindex(fvm_sol)
+for j in eachindex(fvm_sol.u)
     ax = Axis(fig[2, j], width = 600, height = 600,
         xlabel = "x", ylabel = "y",
         title = "t = $(fvm_sol.t[j])",
diff --git a/docs/src/wyos/linear_reaction_diffusion_equations.md b/docs/src/wyos/linear_reaction_diffusion_equations.md
index 38ec02e1..833e8cb3 100644
--- a/docs/src/wyos/linear_reaction_diffusion_equations.md
+++ b/docs/src/wyos/linear_reaction_diffusion_equations.md
@@ -116,7 +116,7 @@ sol |> tc #hide
 ````@example linear_reaction_diffusion_equations
 using CairoMakie
 fig = Figure(fontsize = 38)
-for j in eachindex(sol)
+for j in eachindex(sol.u)
     ax = Axis(fig[1, j], width = 600, height = 600,
         xlabel = "x", ylabel = "y",
         title = "t = $(sol.t[j])")
@@ -248,7 +248,7 @@ sol = solve(prob, Tsit5(); saveat = 2)
 
 using CairoMakie
 fig = Figure(fontsize = 38)
-for j in eachindex(sol)
+for j in eachindex(sol.u)
     ax = Axis(fig[1, j], width = 600, height = 600,
         xlabel = "x", ylabel = "y",
         title = "t = $(sol.t[j])")
diff --git a/src/FiniteVolumeMethod.jl b/src/FiniteVolumeMethod.jl
index 034cf758..ef7a3594 100644
--- a/src/FiniteVolumeMethod.jl
+++ b/src/FiniteVolumeMethod.jl
@@ -1,7 +1,7 @@
 module FiniteVolumeMethod
 
-using ChunkSplitters: ChunkSplitters, chunks
-using CommonSolve: CommonSolve
+using ChunkSplitters: ChunkSplitters, index_chunks
+using CommonSolve: CommonSolve, solve
 using DelaunayTriangulation: DelaunayTriangulation, Triangulation,
     add_ghost_triangles!,
     convert_boundary_points_to_indices,
diff --git a/src/solve.jl b/src/solve.jl
index 8b28a813..4440e536 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -10,9 +10,9 @@ function get_multithreading_parameters(prob::Union{FVMProblem, FVMSystem})
     end
     solid_triangles = collect(each_solid_triangle(prob.mesh.triangulation))
     solid_vertices = collect(DelaunayTriangulation.each_point_index(prob.mesh.triangulation)) # we check for points in the vertex inside the source contribution codes
-    chunked_solid_triangles = chunks(solid_triangles, nt)
+    chunked_solid_triangles = [(range, i) for (i, range) in enumerate(index_chunks(solid_triangles; n = nt))]
     boundary_edges = collect(keys(get_boundary_edge_map(prob.mesh.triangulation)))
-    chunked_boundary_edges = chunks(boundary_edges, nt)
+    chunked_boundary_edges = [(range, i) for (i, range) in enumerate(index_chunks(boundary_edges; n = nt))]
     return (
         duplicated_du = duplicated_du,
         dirichlet_nodes = dirichlet_nodes,
diff --git a/test/Project.toml b/test/Project.toml
index 33c5f4cc..60751dca 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -14,17 +14,20 @@ DisplayAs = "0b91fe84-8a4c-11e9-3e1d-67c38462b6d6"
 ElasticArrays = "fdbdab4c-e67f-52f5-8c3f-e7b388dad3d4"
 ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
+FiniteVolumeMethod = "d4f04ab7-4f65-4d72-8a28-7087bc7f46f4"
 Krylov = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 NaturalNeighbours = "f16ad982-4edb-46b1-8125-78e5a8b5a9e6"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
 PolygonOps = "647866c9-e3ac-4575-94e7-e3d426903924"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
+SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SciMLTesting = "09d9d899-5365-40a9-917a-5f67fddea283"
 SimpleGraphs = "55797a34-41de-5266-9ec1-32ac4eb504d3"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
@@ -55,10 +58,12 @@ LinearSolve = "3.84.0"
 NaturalNeighbours = "1.3.6"
 NonlinearSolve = "4.19.1"
 OrdinaryDiffEq = "7.0.0"
+OrdinaryDiffEqSDIRK = "2.7.0"
 PolygonOps = "0.1.2"
 ReferenceTests = "0.10.6"
 Roots = "3.0.0"
 SafeTestsets = "0.1, 1"
+SciMLOperators = "1.22.0"
 SciMLTesting = "1"
 SimpleGraphs = "0.8.10"
 SpecialFunctions = "2.8.0"
diff --git a/test/equations.jl b/test/equations.jl
index 5b861405..e8372d3c 100644
--- a/test/equations.jl
+++ b/test/equations.jl
@@ -4,6 +4,7 @@ using PolygonOps
 using LinearAlgebra
 using DelaunayTriangulation
 using OrdinaryDiffEq
+using OrdinaryDiffEqSDIRK: TRBDF2
 using LinearSolve
 const FVM = FiniteVolumeMethod
 const DT = DelaunayTriangulation
@@ -112,9 +113,9 @@ end
     solprobser = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.05, parallel = Val(false))
     solprobu = reduce(hcat, solprob.u)
     solprobseru = reduce(hcat, solprobser.u)
-    solsysu1 = reduce(hcat, [solsys.u[i][1, :] for i in eachindex(solsys)])
-    solsysu2 = reduce(hcat, [solsys.u[i][2, :] for i in eachindex(solsys)])
-    solsysser1 = reduce(hcat, [solsysser.u[i][1, :] for i in eachindex(solsysser)])
-    solsysser2 = reduce(hcat, [solsysser.u[i][2, :] for i in eachindex(solsysser)])
+    solsysu1 = reduce(hcat, [solsys.u[i][1, :] for i in eachindex(solsys.u)])
+    solsysu2 = reduce(hcat, [solsys.u[i][2, :] for i in eachindex(solsys.u)])
+    solsysser1 = reduce(hcat, [solsysser.u[i][1, :] for i in eachindex(solsysser.u)])
+    solsysser2 = reduce(hcat, [solsysser.u[i][2, :] for i in eachindex(solsysser.u)])
     @test all(≈(solprobu), (solprobseru, solsysu1, solsysu2, solsysser1, solsysser2))
 end
diff --git a/test/test_figures/conditions.png b/test/test_figures/conditions.png
index fc42eda2..b0fcd97d 100644
Binary files a/test/test_figures/conditions.png and b/test/test_figures/conditions.png differ
diff --git a/test/test_figures/internal_conditions.png b/test/test_figures/internal_conditions.png
index a58f9f3d..44e05bf8 100644
Binary files a/test/test_figures/internal_conditions.png and b/test/test_figures/internal_conditions.png differ