diff --git a/docs/src/examples.md b/docs/src/examples.md
index 65d01b3..eaeda33 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -128,7 +128,8 @@ Now we solve and plot the solution.
 ```julia
 using OrdinaryDiffEq
 using LinearSolve
-sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.001)
+sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.001,
+    abstol = 1.0e-10, reltol = 1.0e-8)
 
 using CairoMakie
 let t_range = LinRange(0.0, final_time, 250)
diff --git a/docs/src/figures/dirichlet_source_surface.png b/docs/src/figures/dirichlet_source_surface.png
index fa81269..93e1720 100644
Binary files a/docs/src/figures/dirichlet_source_surface.png and b/docs/src/figures/dirichlet_source_surface.png differ
diff --git a/docs/src/figures/porous_surface.png b/docs/src/figures/porous_surface.png
index a89a105..29eb7d8 100644
Binary files a/docs/src/figures/porous_surface.png and b/docs/src/figures/porous_surface.png differ
diff --git a/docs/src/figures/robin_diffusion_surface.png b/docs/src/figures/robin_diffusion_surface.png
index d4d06cd..b18453a 100644
Binary files a/docs/src/figures/robin_diffusion_surface.png and b/docs/src/figures/robin_diffusion_surface.png differ
diff --git a/src/precompilation.jl b/src/precompilation.jl
index 1179e27..fb556c5 100644
--- a/src/precompilation.jl
+++ b/src/precompilation.jl
@@ -1,4 +1,4 @@
-using PrecompileTools
+using PrecompileTools: @setup_workload, @compile_workload
 
 @setup_workload begin
     # Minimal setup - these are lightweight and don't add much to precompilation time
diff --git a/test/Project.toml b/test/Project.toml
index 60c1e99..da217d6 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -8,6 +8,7 @@ MethodOfLines = "94925ecb-adb7-4558-8ed8-f975c56a0bf4"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
 ReferenceTests = "324d217c-45ce-50fc-942e-d289b448e8cf"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
diff --git a/test/dirichlet_source.jl b/test/dirichlet_source.jl
index cf4606b..05bc6d5 100644
--- a/test/dirichlet_source.jl
+++ b/test/dirichlet_source.jl
@@ -2,6 +2,7 @@ using ..FiniteVolumeMethod1D
 using LinearAlgebra
 using LinearSolve
 using OrdinaryDiffEq
+using OrdinaryDiffEqSDIRK
 using CairoMakie
 using ReferenceTests
 
@@ -32,8 +33,11 @@ prob = FVMProblem(
     initial_condition,
     final_time
 )
-sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.001)
-exact_sol = [exact_u.(mesh_points, sol.t[i]) for i in eachindex(sol)]
+sol = solve(
+    prob, TRBDF2(linsolve = KLUFactorization()),
+    saveat = 0.001, abstol = 1.0e-10, reltol = 1.0e-8
+)
+exact_sol = [exact_u.(mesh_points, sol.t[i]) for i in eachindex(sol.t)]
 @test reduce(hcat, sol.u) ≈ reduce(hcat, exact_sol) rtol = 1.0e-1
 
 let t_range = LinRange(0.0, final_time, 250)
diff --git a/test/fisher.jl b/test/fisher.jl
index b72676a..6ef6eb6 100644
--- a/test/fisher.jl
+++ b/test/fisher.jl
@@ -1,6 +1,7 @@
 using ..FiniteVolumeMethod1D
 using MethodOfLines
 using OrdinaryDiffEq
+using OrdinaryDiffEqSDIRK
 using ReferenceTests
 using CairoMakie
 using ModelingToolkit
diff --git a/test/heat.jl b/test/heat.jl
index a2d73b9..e3d5832 100644
--- a/test/heat.jl
+++ b/test/heat.jl
@@ -1,5 +1,6 @@
 using ..FiniteVolumeMethod1D
 using OrdinaryDiffEq
+using OrdinaryDiffEqSDIRK
 using ReferenceTests
 using CairoMakie
 using LinearSolve
@@ -27,7 +28,7 @@ prob = FVMProblem(
 )
 
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()))
-exact_sol = [exact2.(mesh_points, sol.t[i]) for i in eachindex(sol)]
+exact_sol = [exact2.(mesh_points, sol.t[i]) for i in eachindex(sol.t)]
 @test reduce(hcat, sol.u) ≈ reduce(hcat, exact_sol) rtol = 1.0e-2
 
 let t_range = LinRange(0.0, final_time, 250)
diff --git a/test/porous_medium.jl b/test/porous_medium.jl
index 64ac28d..c51230a 100644
--- a/test/porous_medium.jl
+++ b/test/porous_medium.jl
@@ -1,5 +1,6 @@
 using ..FiniteVolumeMethod1D
 using OrdinaryDiffEq
+using OrdinaryDiffEqSDIRK
 using CairoMakie
 using ReferenceTests
 using SpecialFunctions
@@ -42,7 +43,7 @@ prob = FVMProblem(
     reaction_function = (u, x, t, p) -> zero(u)
 )
 sol = solve(prob, TRBDF2(); saveat = saveat)
-exact = [exact_solution.(mesh_points, sol.t[i]) for i in eachindex(sol)]
+exact = [exact_solution.(mesh_points, sol.t[i]) for i in eachindex(sol.t)]
 @test exact ≈ sol.u rtol = 1.0e-1
 
 fig = Figure(size = (2150, 460), fontsize = 34)
@@ -51,14 +52,14 @@ ax = Axis(
     titlealign = :left, width = 600, height = 300
 )
 colors = [:red, :black, :blue, :darkgreen, :magenta, :orange]
-[lines!(ax, mesh_points, sol.u[i], color = colors[i]) for i in eachindex(sol)]
+[lines!(ax, mesh_points, sol.u[i], color = colors[i]) for i in eachindex(sol.t)]
 ylims!(ax, -1.0e-6, 5)
 xlims!(ax, -1, 1)
 ax = Axis(
     fig[1, 2], xlabel = L"x", ylabel = L"u(x)", title = L"(b):$ $ Exact solution",
     titlealign = :left, width = 600, height = 300
 )
-[lines!(ax, mesh_points, exact[i], color = colors[i]) for i in eachindex(sol)]
+[lines!(ax, mesh_points, exact[i], color = colors[i]) for i in eachindex(sol.t)]
 ylims!(ax, -1.0e-6, 5)
 xlims!(ax, -1, 1)
 ax = Axis(
@@ -67,7 +68,7 @@ ax = Axis(
 )
 [
     lines!(ax, mesh_points, abs.(exact[i] .- sol.u[i]), color = colors[i])
-        for i in eachindex(sol)
+        for i in eachindex(sol.t)
 ]
 ylims!(ax, -1.0e-6, 0.5)
 xlims!(ax, -1, 1)
diff --git a/test/robin_diffusion.jl b/test/robin_diffusion.jl
index d79c922..8c8cd9e 100644
--- a/test/robin_diffusion.jl
+++ b/test/robin_diffusion.jl
@@ -2,6 +2,7 @@ using ..FiniteVolumeMethod1D
 using LinearAlgebra
 using LinearSolve
 using OrdinaryDiffEq
+using OrdinaryDiffEqSDIRK
 using NonlinearSolve
 using CairoMakie
 using ReferenceTests
@@ -49,7 +50,7 @@ prob = FVMProblem(
 sol = solve(prob, TRBDF2(linsolve = KLUFactorization()), saveat = 0.01)
 
 μ = compute_μₙ.(1:100)
-exact_sol = [exact_u.(mesh_points, sol.t[i], Ref(μ)) for i in eachindex(sol)]
+exact_sol = [exact_u.(mesh_points, sol.t[i], Ref(μ)) for i in eachindex(sol.t)]
 @test reduce(hcat, sol.u) ≈ reduce(hcat, exact_sol) rtol = 1.0e-2
 
 let t_range = LinRange(0.0, final_time, 250)