diff --git a/Project.toml b/Project.toml
index 48f7eec0..da14ab5e 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,11 +9,13 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 
 [compat]
@@ -22,10 +24,12 @@ LinearAlgebra = "1.7"
 LinearSolve = "2.21"
 MuladdMacro = "0.2.1"
 OrdinaryDiffEq = "6.59"
+RecipesBase = "1.3"
 Reexport = "1"
 SciMLBase = "2"
 SimpleUnPack = "1"
 SparseArrays = "1.7"
 StaticArrays = "1.5"
+Statistics = "1"
 SymbolicIndexingInterface = "0.2, 0.3"
 julia = "1.9"
diff --git a/docs/Project.toml b/docs/Project.toml
index c7256d09..b3ace33f 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
@@ -7,11 +8,13 @@ LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [compat]
 BenchmarkTools = "1"
+DiffEqDevTools = "2.44"
 Documenter = "1"
 InteractiveUtils = "1"
 LinearAlgebra = "1.7"
@@ -19,5 +22,6 @@ LinearSolve = "2.21"
 OrdinaryDiffEq = "6.59"
 Pkg = "1"
 Plots = "1"
+PrettyTables = "2.3"
 SparseArrays = "1.7"
 StaticArrays = "1.5"
diff --git a/docs/make.jl b/docs/make.jl
index d59aa7f8..6e6c0988 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -81,6 +81,12 @@ makedocs(modules = [PositiveIntegrators],
                  "Heat Equation, Neumann BCs" => "heat_equation_neumann.md",
                  "Heat Equation, Dirichlet BCs" => "heat_equation_dirichlet.md",
              ],
+             "Benchmarks" => [
+                 "Experimental order of convergence" => "convergence.md",
+                 "NPZD model" => "npzd_model_benchmark.md",
+                 "Robertson problem" => "robertson_benchmark.md",
+                 "Stratospheric reaction problem" => "stratospheric_reaction_benchmark.md",
+             ],
              "Troubleshooting, FAQ" => "faq.md",
              "API reference" => "api_reference.md",
              "Contributing" => "contributing.md",
diff --git a/docs/src/api_reference.md b/docs/src/api_reference.md
index 75672ae5..ac1a6f9f 100644
--- a/docs/src/api_reference.md
+++ b/docs/src/api_reference.md
@@ -42,4 +42,12 @@ SSPMPRK43
 ```@docs
 isnegative
 isnonnegative
+rel_max_error_tend
+rel_max_error_overall
+rel_l1_error_tend
+rel_l2_error_tend
+work_precision_adaptive
+work_precision_adaptive!
+work_precision_fixed
+work_precision_fixed!
 ```
diff --git a/docs/src/convergence.md b/docs/src/convergence.md
new file mode 100644
index 00000000..9f9ed239
--- /dev/null
+++ b/docs/src/convergence.md
@@ -0,0 +1,180 @@
+# [Experimental convergence order of MPRK schemes](@id convergence_mprk)
+
+In this tutorial we check that the implemented MPRK schemes have the expected order of convergence. 
+
+## Conservative production-destruction systems
+
+First, we consider conservative production-destruction systems (PDS). To investigate the convergence order we define the non-autonomous test problem 
+
+```math
+\begin{aligned}
+u_1' &= \cos(\pi t)^2 u_2 - \sin(2\pi t)^2 u_1,\\
+u_2' & = \sin(2\pi t)^2 u_1 - \cos(\pi t)^2 u_2.
+\end{aligned}
+```
+
+The PDS is conservative, since the sum of the right hand side terms equals zero. 
+An implementation of the problem is given next.
+
+
+```@example eoc
+using PositiveIntegrators
+
+# define problem
+P(u, p, t) = [0.0 cos.(π * t) .^ 2 * u[2]; sin.(2 * π * t) .^ 2 * u[1] 0.0]
+prob = ConservativePDSProblem(P, [0.9; 0.1], (0.0, 1.0))
+
+nothing # hide output
+```
+
+To use `analyticless_test_convergence` from [`DiffEqDevTools`](https://github.com/SciML/DiffEqDevTools.jl) we need to pick a solver to compute the reference solution and specify tolerances. Since the problem is not stiff we use the high order explicit solver `Vern9()` from [`OrdinaryDiffEq`](https://docs.sciml.ai/OrdinaryDiffEq/stable/). Moreover, we need to choose the different time step sizes which are used to investigate the convergence behavior. 
+
+```@example eoc
+using OrdinaryDiffEq
+using DiffEqDevTools # load analyticless_test_convergence
+
+# solver and tolerances to compute reference solution
+test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
+
+# choose step sizes
+dts = 0.5 .^ (5:10)
+
+nothing # hide
+```
+
+### Second order MPRK schemes
+
+First, we test several second order MPRK schemes.
+
+```@example eoc
+# select schemes
+algs2 = [MPRK22(0.5); MPRK22(2.0 / 3.0); MPRK22(1.0); SSPMPRK22(0.5, 1.0)]
+labels2 = ["MPRK22(0.5)"; "MPRK22(2.0/3.0)"; "MPRK22(1.0)"; "SSPMPRK22(0.5, 1.0)"]
+
+#compute errors and experimental order of convergence
+err_eoc = []
+for i in eachindex(algs2)
+     sim = analyticless_test_convergence(dts, prob, algs2[i], test_setup)
+
+     err = sim.errors[:l∞]
+     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+
+     push!(err_eoc, tuple.(err, eoc))
+end
+```
+
+Next, we print a table with the computed data. The table lists the errors obtained with the respective time step size ``Δ t`` as well as the estimated order of convergence in parenthesis.
+
+```@example eoc
+using Printf # load @sprintf
+using PrettyTables # load pretty_table
+
+# gather data for table
+data = hcat(dts, reduce(hcat,err_eoc))
+
+# print table
+formatter = (v, i, j) ->  (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v)
+pretty_table(data, formatters = formatter, header = ["Δt"; labels2])                  
+```
+
+The table shows that all schemes converge as expected.
+
+### Third order MPRK schemes
+
+In this section, we proceed as above, but consider third order MPRK schemes instead.
+
+```@example eoc
+# select 3rd order schemes
+algs3 = [MPRK43I(1.0, 0.5); MPRK43I(0.5, 0.75); MPRK43II(0.5); MPRK43II(2.0 / 3.0); 
+         SSPMPRK43()]
+labels3 = ["MPRK43I(1.0,0.5)"; "MPRK43I(0.5, 0.75)"; "MPRK43II(0.5)"; "MPRK43II(2.0/3.0)";
+          "SSPMPRK43()"]
+
+#compute errors and experimental order of convergence
+err_eoc = []
+for i in eachindex(algs3)
+     sim = analyticless_test_convergence(dts, prob, algs3[i], test_setup)
+
+     err = sim.errors[:l∞]
+     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+
+     push!(err_eoc, tuple.(err, eoc))
+end
+
+# gather data for table
+data = hcat(dts, reduce(hcat,err_eoc))
+
+# print table
+formatter = (v, i, j) ->  (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v)
+pretty_table(data, formatters = formatter, header = ["Δt"; labels3])  
+```
+
+As above, the table shows that all schemes converge as expected.
+
+## Non-conservative PDS
+
+In this section we consider the non-autonomous but non-conservative test problem 
+
+```math
+\begin{aligned}
+u_1' &= \cos(\pi t)^2 u_2 - \sin(2\pi t)^2 u_1 - \cos(2\pi t)^2 u_1,\\
+u_2' & = \sin(2\pi t)^2 u_1 - \cos(\pi t)^2 u_2 - \sin(\pi t)^2 u_2.
+\end{aligned}
+```
+
+Since the sum of the right hand side terms don't cancel, the PDS is indeed non-conservative. Hence, we need to use `PDSProblem` for its implementation.
+
+```@example eoc
+# choose problem
+P(u, p, t) = [0.0 cos.(π * t) .^ 2 * u[2]; sin.(2 * π * t) .^ 2 * u[1] 0.0]
+D(u, p, t) = [cos.(2 * π * t) .^ 2 * u[1]; sin.(π * t) .^ 2 * u[2]]
+prob = PDSProblem(P, D, [0.9; 0.1], (0.0, 1.0))
+
+nothing # hide
+```
+
+The following sections will show that the selected MPRK schemes show the expected convergence order also for this non-conservative PDS.
+
+### Second order MPRK schemes
+
+```@example eoc
+#compute errors and experimental order of convergence
+err_eoc = []
+for i in eachindex(algs2)
+     sim = analyticless_test_convergence(dts, prob, algs2[i], test_setup)
+
+     err = sim.errors[:l∞]
+     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+
+     push!(err_eoc, tuple.(err, eoc))
+end
+
+# gather data for table
+data = hcat(dts, reduce(hcat,err_eoc))
+
+# print table
+formatter = (v, i, j) ->  (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v)
+pretty_table(data, formatters = formatter, header = ["Δt"; labels2])                  
+```
+
+### Third order MPRK schemes
+
+```@example eoc
+#compute errors and experimental order of convergence
+err_eoc = []
+for i in eachindex(algs3)
+     sim = analyticless_test_convergence(dts, prob, algs3[i], test_setup)
+
+     err = sim.errors[:l∞]
+     eoc = [NaN; -log2.(err[2:end] ./ err[1:(end - 1)])]
+
+     push!(err_eoc, tuple.(err, eoc))
+end
+
+# gather data for table
+data = hcat(dts, reduce(hcat,err_eoc))
+
+# print table
+formatter = (v, i, j) ->  (j>1) ? (@sprintf "%5.2e (%4.2f) " v[1] v[2]) : (@sprintf "%5.2e " v)
+pretty_table(data, formatters = formatter, header = ["Δt"; labels3])  
+```
\ No newline at end of file
diff --git a/docs/src/npzd_model_benchmark.md b/docs/src/npzd_model_benchmark.md
new file mode 100644
index 00000000..230d408b
--- /dev/null
+++ b/docs/src/npzd_model_benchmark.md
@@ -0,0 +1,392 @@
+# [Benchmark: Solution of an NPZD model](@id benchmark-npzd)
+
+We use the NPZD model [`prob_pds_npzd`](@ref) to assess the efficiency of different solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) and [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl).
+
+```@example NPZD
+using OrdinaryDiffEq, PositiveIntegrators
+
+# select NPZD problem
+prob = prob_pds_npzd
+nothing # hide
+```
+
+To keep the following code as clear as possible, we define a helper function `npzd_plot` that we use for plotting.
+
+```@example NPZD
+using Plots
+
+npzd_plot = function(sol, sol_ref = nothing, title = "")
+     colors = palette(:default)[1:4]'
+     if !isnothing(sol_ref)
+          p = plot(sol_ref, linestyle = :dash, label = "", color = colors, 
+          linewidth = 2)
+          plot!(p, sol; denseplot = false, markers = :circle, ylims = (-1.0, 10.0),
+          color = colors, title, label = ["N" "P" "Z" "D"], legend = :right,
+          linewidth = 2);
+     else
+          p = plot(sol; denseplot = false, markers = :circle, ylims = (-1.0, 10.0),
+          color = colors, title, label = ["N" "P" "Z" "D"], legend = :right, 
+          linewidths = 2);
+     end
+     return p
+end
+nothing  # hide
+```
+
+Standard methods have difficulties to solve the NPZD problem accurately for loose tolerances or large time step sizes. 
+This is because once negative values in the ``N``-component occur, this will lead to a further decrease in ``N`` and thus completely inaccurate solutions.  
+
+```@example NPZD
+# compute reference solution for plotting
+ref_sol = solve(prob, Vern7(); abstol = 1e-14, reltol = 1e-13);
+
+# compute solutions with loose tolerances
+abstol = 1e-2
+reltol = 1e-1
+sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol);
+sol_MPRK = solve(prob, MPRK22(1.0); abstol, reltol);
+
+# plot solutions
+p1 = npzd_plot(sol_Ros23, ref_sol, "Rosenbrock23"); # helper function defined above
+p2 = npzd_plot(sol_MPRK, ref_sol, "MPRK22(1.0)"); 
+plot(p1, p2)
+```
+
+Nevertheless, [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) provides the solver option `isoutofdomain`, which can be used in combination with [`isnegative`](@ref) to guarantee nonnegative solutions. 
+
+```@example NPZD
+sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol, 
+                  isoutofdomain = isnegative); #reject negative solutions
+
+npzd_plot(sol_Ros23, ref_sol) #auxiliary function defined above
+```
+
+## Work-Precision diagrams
+
+In the following sections we show several work-precision diagrams, which compare different methods with respect to computing time and error. 
+First we focus on adaptive methods, afterwards we also show results obtained with fixed time step sizes.
+
+Since the NPZD problem is not stiff, we can use an explicit high-order scheme to compute a reference solution.
+
+```@example NPZD
+# select solver to compute reference solution
+alg_ref = Vern7()
+nothing  # hide
+```
+
+### Adaptive schemes
+
+We use the functions [`work_precision_adaptive`](@ref) and [`work_precision_adaptive!`](@ref) to compute the data for the diagrams.
+Furthermore, the following absolute and relative tolerances are used.
+
+```@example NPZD
+# set absolute and relative tolerances
+abstols = 1.0 ./ 10.0 .^ (2:1:8)
+reltols = abstols .* 10.0
+nothing  # hide
+```
+
+#### Relative maximum error at the final time
+
+In this section the chosen error is the relative maximum error at the final time ``t = 10.0``.
+
+```@example NPZD
+# select relative maximum error at the end of the problem's time span.
+compute_error = rel_max_error_tend
+nothing # hide
+```
+
+We start with a comparison of different adaptive MPRK schemes.
+
+```@example NPZD
+# choose methods to compare
+algs = [MPRK22(0.5); MPRK22(2.0 / 3.0); MPRK22(1.0); SSPMPRK22(0.5, 1.0); 
+        MPRK43I(1.0, 0.5); MPRK43I(0.5, 0.75); MPRK43II(0.5); MPRK43II(2.0 / 3.0)]
+labels = ["MPRK22(0.5)"; "MPPRK22(2/3)"; "MPRK22(1.0)"; "SSPMPRK22(0.5,1.0)"; 
+          "MPRK43I(1.0, 0.5)"; "MPRK43I(0.5, 0.75)"; "MPRK43II(0.5)"; "MPRK43II(2.0/3.0)"]
+
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
+                               compute_error)
+
+# plot work-precision diagram
+plot(wp, labels; title = "NPZD benchmark", legend = :topright,     
+     color = permutedims([repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...]),
+     xlims = (10^-7, 2*10^-1), xticks = 10.0 .^ (-8:1:0),
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+The second- and third-order methods behave very similarly. For comparisons with other schemes from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) we choose `MPRK22(1.0)` and `MPRK43I(1.0, 0.5)`.
+
+```@example NPZD
+sol_MPRK22 = solve(prob, MPRK22(1.0); abstol, reltol)
+sol_MPRK43 = solve(prob, MPRK43I(1.0, 0.5); abstol, reltol)
+
+p1 = npzd_plot(sol_MPRK22, ref_sol, "MPRK22(1.0)");
+p2 = npzd_plot(sol_MPRK43, ref_sol, "MPRK43I(1.0, 0.5)");
+plot(p1, p2)
+```
+
+Next we compare `MPRK22(1.0)` and `MPRK43I(1.0, 0.5)` to explicit and implicit methods of second and third order from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). 
+To guarantee nonnegative solutions, we select the solver option `isoutofdomain = isnegative`.
+
+```@example NPZD
+# select MPRK methods for reference
+algs1 = [MPRK22(1.0); MPRK43I(1.0, 0.5)]
+labels1 = ["MPRK22(1.0)"; "MPRK43I(1.0,0.5)"]
+
+# select methods from OrdinaryDiffEq
+algs2 = [Midpoint(); Heun(); Ralston(); TRBDF2(); SDIRK2(); Kvaerno3(); KenCarp3(); Rodas3();
+         ROS2(); ROS3(); Rosenbrock23()]
+labels2 = ["Midpoint"; "Heun"; "Ralston"; "TRBDF2"; "SDIRK2"; "Kvearno3"; "KenCarp3"; "Rodas3";
+          "ROS2"; "ROS3"; "Rosenbrock23"]
+
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
+                               compute_error)
+# add work-precision data with isoutofdomain=isnegative
+work_precision_adaptive!(wp, prob, algs2, labels2, abstols, reltols, alg_ref;
+                               compute_error, isoutofdomain=isnegative)
+
+plot(wp, [labels1; labels2]; title = "NPZD benchmark", legend = :topright,
+     color = permutedims([1, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 4)...]),
+     xlims = (5*10^-8, 2*10^-1), xticks = 10.0 .^ (-8:1:0),
+     ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+We see that for the NPZD problem the use of adaptive MPRK schemes is only beneficial when using the loosest tolerances.
+
+Now we compare `MPRK22(1.0)` and `MPRK43I(1.0, 0.5)` to [recommended solvers](https://docs.sciml.ai/DiffEqDocs/dev/solvers/ode_solve/) from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). Again, to guarantee positive solutions we select the solver option `isoutofdomain = isnegative`.
+
+```@example NPZD
+algs3 = [Tsit5(); BS3(); Vern6(); Vern7(); Vern8(); TRBDF2(); Rosenbrock23(); 
+         Rodas5P(); Rodas4P()]
+labels3 = ["Tsit5"; "BS3"; "Vern6"; "Vern7"; "Vern8"; "TRBDF2"; "Rosenbrock23";
+          "Rodas5P"; "Rodas4P"]
+
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
+                               compute_error) 
+# add work-precision data with isoutofdomain = isnegative
+work_precision_adaptive!(wp, prob, algs3, labels3, abstols, reltols, alg_ref;
+                               compute_error, isoutofdomain = isnegative)
+
+# plot work-precision diagram
+plot(wp, [labels1; labels3]; title = "NPZD benchmark", legend = :topright,
+     color = permutedims([1, 3, repeat([4], 5)...,5, repeat([6], 1)...,repeat([7],2)...]),
+     xlims = (10^-11, 10^1), xticks = 10.0 .^ (-11:1:1),
+     ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+We see that it is advisable to use a high order explicit method like `Vern7()` with `isoutofdomain = isnegative` to obtain nonnegative solutions of such a non-stiff problem.
+
+#### Relative maximum error over all time steps
+
+In this section we do not compare the relative maximum errors at the final time ``t = 10.0``, but the relative maximum errors over all time steps. 
+
+```@example NPZD
+# select relative maximum error over all time steps
+compute_error = rel_max_error_overall
+nothing  # hide
+```
+
+The results are very similar to those from above. 
+We therefore only show the work-precision diagrams without further comments. 
+The main difference are significantly increased errors which mainly occur around time ``t = 2.0`` where there is a sharp kink in the solution.
+
+```@example NPZD
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
+                               compute_error)
+
+# plot work-precision diagram
+plot(wp, labels; title = "NPZD benchmark", legend = :topright,     
+          color = permutedims([repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...]),
+          xlims = (10^-5, 10^4), xticks = 10.0 .^ (-5:1:4),
+          ylims = (10^-6, 10^-1), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+```@example NPZD
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
+                               compute_error)
+# add work-precision data with isoutofdomain = isnegative
+work_precision_adaptive!(wp, prob, algs2, labels2, abstols, reltols, alg_ref;
+                               compute_error, isoutofdomain=isnegative)
+
+# plot work-precision diagram
+plot(wp, [labels1; labels2]; title = "NPZD benchmark", legend = :topright,
+     color = permutedims([1, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 4)...]),
+     xlims = (10^-5, 10^4), xticks = 10.0 .^ (-5:1:4),
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-6:1:0), minorticks = 10)
+```
+
+```@example NPZD
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
+                               compute_error)
+# add work-precision data with isoutofdomain = isnegative                             
+work_precision_adaptive!(wp, prob, algs3, labels3, abstols, reltols, alg_ref;
+                               compute_error, isoutofdomain=isnegative)
+
+# plot work-precision diagram
+plot(wp, [labels1; labels3]; title = "NPZD benchmark", legend = :topright,
+     color = permutedims([1, 3, repeat([4], 5)...,5, repeat([6], 1)...,repeat([7],2)...]),
+     xlims = (10^-7, 10^5), xticks = 10.0 .^ (-7:1:5),
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-6:1:0), minorticks = 10)
+```
+
+### Fixed time steps sizes
+
+Here we use fixed time step sizes instead of adaptive time stepping. 
+Similar to the adaptive situation above, standard schemes are likely to compute negative solutions for the NPZD problem. 
+
+```@example NPZD
+sol_Ros23 = solve(prob, Rosenbrock23(), dt = 1.0, adaptive = false);
+sol_MPRK = solve(prob, MPRK22(1.0), dt = 1.0, adaptive = false);
+
+p1 = npzd_plot(sol_Ros23, ref_sol, "Rosenbrock23");
+p2 = npzd_plot(sol_MPRK, ref_sol, "MPRK22(1.0)");
+plot(p1, p2)
+```
+
+We use the functions [`work_precision_fixed`](@ref) and [`work_precision_fixed!`](@ref) to compute the data for the diagrams.
+Please note that these functions set error and computing time to `Inf`, whenever a solution contains negative elements. 
+Consequently, such cases are not visible in the work-precision diagrams.
+
+Within the work-precision diagrams we use the following time step sizes.
+
+```@example NPZD
+# set time step sizes
+dts = 1.0 ./ 2.0 .^ (0:1:12)
+nothing # hide
+```
+
+#### Relative maximum error at the end of the problem's time span
+
+Again, we start with the relative maximum error at the final time ``t = 10.0``.
+
+```@example NPZD
+# select relative maximum error at the end of the problem's time span.
+compute_error = rel_max_error_tend
+nothing  # hide
+```
+
+First, we compare different MPRK methods. 
+For fixed time step sizes we can also consider `MPE()` and `SSPMPRK43()`.
+
+```@example NPZD
+# choose MPRK methods to compare
+algs = [MPE(); algs; SSPMPRK43()]
+labels = ["MPE()"; labels; "SSPMPRK43"]
+
+# compute work-precision data
+wp = work_precision_fixed(prob, algs, labels, dts, alg_ref;
+                               compute_error)
+
+# plot work-precision diagram
+plot(wp, labels; title = "NPZD benchmark", legend = :bottomleft,     
+     color = permutedims([5,repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...,6]),
+     xlims = (10^-10, 1*10^0), xticks = 10.0 .^ (-10:1:0),
+     ylims = (1*10^-6, 10^-1), yticks = 10.0 .^ (-6:1:0), minorticks = 10)
+```
+
+Apart from `MPE()` the schemes behave very similar and a difference in order can only be observed for the smaller step sizes. 
+We choose `MPRK22(1.0)` and `MPRK43I(1.0, 0.5)` for comparisons with other second and third order schemes from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
+
+```@example NPZD
+# compute work-precision data
+wp = work_precision_fixed(prob, [algs1; algs2], [labels1; labels2], dts, alg_ref;
+                               compute_error)
+
+# plot work-precision diagram
+plot(wp, [labels1; labels2]; title = "NPZD benchmark", legend = :topright,
+     color = permutedims([1, 3, repeat([4], 3)..., repeat([5],4)...,repeat([6],4)...]),
+     xlims = (10^-13, 10^2), xticks = 10.0 .^ (-12:2:6),
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)         
+```
+
+We see that the MPRK schemes are to be preferred for the rather large step sizes ``\Delta t \in\lbrace 1.0, 0.5, 0.25, 0.125\rbrace``, for which the other schemes cannot provide nonnegative solutions.
+
+```@example NPZD
+# solution computed with MPRK43I(1.0, 0.5) and dt = 0.125
+sol_MPRK = solve(prob, MPRK43I(1.0, 0.5); dt = dts[4], adaptive = false);
+
+# plot solution
+npzd_plot(sol_MPRK, ref_sol)
+```
+
+Finally, we show a comparison of `MPRK22(1.0)`, `MPRK43I(1.0, 0.5)` and [recommended solvers](https://docs.sciml.ai/DiffEqDocs/dev/solvers/ode_solve/) from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/).
+
+```@example NPZD
+# compute work-precision data
+wp = work_precision_fixed(prob, [algs1; algs3], [labels1; labels3], dts, alg_ref;
+                               compute_error)
+
+# plot work-precision diagram
+plot(wp, [labels1; labels3]; title = "NPZD benchmark", legend = :topright,
+     color = permutedims([1, 3, repeat([4], 3)..., repeat([5],4)...,repeat([6],4)...]),
+     xlims = (10^-14, 10^0), xticks = 10.0 .^ (-14:2:10),
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)         
+```
+
+#### Relative maximum error over all time steps
+
+As for the adaptive schemes, we also show work-precisions diagrams where the error is the relative maximum error over all time steps.
+
+```@example NPZD
+# select relative maximum error over all time steps
+compute_error = rel_max_error_overall
+nothing  # hide
+```
+
+```@example NPZD
+
+# compute work-precision
+wp = work_precision_fixed(prob, algs, labels, dts, alg_ref;
+                               compute_error)
+
+#plot work-precision diagram
+plot(wp, labels; title = "NPZD benchmark", legend = :bottomleft,     
+     color = permutedims([5,repeat([1], 3)..., 2, repeat([3], 2)..., repeat([4], 2)...,6]),
+     xlims = (10^-4, 10^5), xticks = 10.0 .^ (-4:1:5),
+     ylims = (10^-6, 10^-1), yticks = 10.0 .^ (-6:1:0), minorticks = 10)
+```
+
+```@example NPZD
+wp = work_precision_fixed(prob, algs1, labels1, dts, alg_ref;
+                               compute_error)
+work_precision_fixed!(wp, prob, algs2, labels2, dts, alg_ref;
+                     compute_error)                               
+
+plot(wp, [labels1; labels2]; title = "NPZD benchmark", legend = :topright,
+     color = permutedims([1, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 4)...]),
+     xlims = (10^-4, 10^6), xticks = 10.0 .^ (-12:1:6),
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)         
+```
+
+```@example NPZD
+wp = work_precision_fixed(prob, algs1, labels1, dts, alg_ref;
+                               compute_error)
+work_precision_fixed!(wp, prob, algs3, labels3, dts, alg_ref;
+                     compute_error)                               
+
+plot(wp, [labels1; labels3]; title = "NPZD benchmark", legend = :bottomleft,
+     color = permutedims([1, 3, repeat([4], 5)..., 5, repeat([7], 3)...]),
+     xlims = (10^-12, 10^6), xticks = 10.0 .^ (-12:2:6),
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)         
+```
+
+## Package versions
+
+These results were obtained using the following versions.
+```@example NPZD
+using InteractiveUtils
+versioninfo()
+println()
+
+using Pkg
+Pkg.status(["PositiveIntegrators", "StaticArrays", "LinearSolve", "OrdinaryDiffEq"],
+           mode=PKGMODE_MANIFEST)
+nothing # hide
+```
diff --git a/docs/src/robertson.md b/docs/src/robertson.md
index 886fe3a9..8bc2b8a3 100644
--- a/docs/src/robertson.md
+++ b/docs/src/robertson.md
@@ -44,7 +44,7 @@ tspan = (0.0, 1.0e11)
 prob = ConservativePDSProblem(prod, u0, tspan)
 
 sol = solve(prob, MPRK43I(1.0, 0.5))
-nothing #hide
+nothing  # hide
 ```
 ```@example robertson
 using Plots
@@ -72,7 +72,7 @@ stepsize_callback = DiscreteCallback(
     save_positions = (false, false),
     initialize = (c, u, t, integrator) -> set_proposed_dt!(integrator, 1.0e-5))
 sol_cb = solve(prob, SSPMPRK43(); dt = Inf, callback =  stepsize_callback);
-nothing #hide
+nothing  # hide
 ```
 ```@example robertson
 plot(sol_cb, tspan = (1e-6, 1e11),  xaxis = :log,
diff --git a/docs/src/robertson_benchmark.md b/docs/src/robertson_benchmark.md
new file mode 100644
index 00000000..764d71a2
--- /dev/null
+++ b/docs/src/robertson_benchmark.md
@@ -0,0 +1,268 @@
+# [Benchmark: Solution of the Robertson problem](@id benchmark-robertson)
+
+Here we use the stiff Robertson problem [`prob_pds_robertson`](@ref) to assess the efficiency of different solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) and [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl).
+
+```@example ROBER
+using OrdinaryDiffEq, PositiveIntegrators
+
+# select Robertson problem
+prob = prob_pds_robertson
+nothing # hide
+```
+
+To keep the following code as clear as possible, we define a helper function `robertson_plot` that we use for plotting.
+
+```@example ROBER
+using Plots
+
+robertson_plot = function (sol, sol_ref = nothing, title = "")
+    colors = palette(:default)[1:3]'
+    if !isnothing(sol_ref)
+        p = plot(sol_ref, tspan = (1e-6, 1e11), xaxis = :log,
+                 idxs = [(0, 1), ((x, y) -> (x, 1e4 .* y), 0, 2), (0, 3)],
+                 linestyle = :dash, label = "", color = colors, linewidth = 2)
+        plot!(p, sol; tspan = (1e-6, 1e11), xaxis = :log, denseplot = false,
+              markers = :circle, ylims = (-0.2, 1.2),
+              idxs = [(0, 1), ((x, y) -> (x, 1e4 .* y), 0, 2), (0, 3)],
+              title, xticks = 10.0 .^ (-6:4:10), color = colors,
+              linewidht = 2, legend = :right, label = ["u₁" "u₂" "u₃"])
+    else
+        p = plot(sol; tspan = (1e-6, 1e11), xaxis = :log, denseplot = false,
+                 markers = :circle, ylims = (-0.2, 1.2),
+                 idxs = [(0, 1), ((x, y) -> (x, 1e4 .* y), 0, 2), (0, 3)],
+                 title, xticks = 10.0 .^ (-6:4:10), color = colors,
+                 linewidht = 2, legend = :right, label = ["u₁" "u₂" "u₃"])
+    end
+    return p
+end
+nothing # hide
+```
+
+For this stiff problem the computation of negative approximations may lead to inaccurate solutions. 
+This typically occurs when adaptive time stepping uses loose tolerances.
+
+```@example ROBER
+# compute reference solution for plotting
+ref_sol = solve(prob, Rodas4P(); abstol = 1e-14, reltol = 1e-13);
+
+# compute solutions with loose tolerances
+abstol = 1e-2
+reltol = 1e-1
+sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol);
+sol_MPRK = solve(prob, MPRK22(1.0); abstol, reltol);
+
+# plot solutions
+p1 = robertson_plot(sol_Ros23, ref_sol, "Rosenbrock23");
+p2 = robertson_plot(sol_MPRK, ref_sol, "MPRK22(1.0)");
+plot(p1, p2)
+```
+
+Nevertheless, [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) provides the solver option `isoutofdomain`, which can be used in combination with [`isnegative`](@ref) to guarantee nonnegative solutions. 
+
+```@example ROBER
+sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol, 
+                  isoutofdomain = isnegative) #reject negative solutions
+
+robertson_plot(sol_Ros23, ref_sol, "Rosenbrock23")
+```
+
+## Work-Precision diagrams
+
+In the following we show several work-precision diagrams, which compare different methods with respect to computing time and the respective error. 
+We focus solely on adaptive methods, since the time interval ``(0, 10^{11})`` is too large to generate accurate solutions with fixed step sizes.
+
+Since the Robertson problem is stiff, we need to use a suited implicit scheme to compute a reference solution, see the [solver guide](https://docs.sciml.ai/DiffEqDocs/dev/solvers/ode_solve/#Stiff-Problems). Note that we cannot use the recommended method `radau()`, since [`prob_pds_robertson`](@ref) uses [StaticArrays](https://juliaarrays.github.io/StaticArrays.jl/stable/) instead of arrays of type `Float64`.
+
+```@example ROBER
+# select solver to compute reference solution
+alg_ref = Rodas4P()
+nothing # hide
+```
+
+We use the functions [`work_precision_adaptive`](@ref) and [`work_precision_adaptive!`](@ref) to compute the data for the diagrams.
+Furthermore, the following absolute and relative tolerances are used.
+
+```@example ROBER
+# set absolute and relative tolerances
+abstols = 1.0 ./ 10.0 .^ (2:1:10)
+reltols = abstols .* 10.0
+nothing # hide
+```
+
+### Relative maximum error at the final time
+
+In this section the chosen error is the relative maximum error at the final time ``t = 10^{11}``.
+
+```@example ROBER
+# select relative maximum error at the end of the problem's time span.
+compute_error = rel_max_error_tend
+nothing # hide
+```
+
+We start with a comparison of different adaptive MPRK schemes.
+```@example ROBER
+# choose methods to compare
+algs = [MPRK22(0.5); MPRK22(2.0 / 3.0); MPRK22(1.0); MPRK43I(1.0, 0.5); MPRK43I(0.5, 0.75);
+        MPRK43II(0.5); MPRK43II(2.0 / 3.0)]
+labels = ["MPRK22(0.5)"; "MPPRK22(2/3)"; "MPRK22(1.0)"; "MPRK43I(1.0,0.5)"; "MPRK43I(0.5,0.75)";
+         "MPRK43II(0.5)"; "MPRK43II(2.0/3.0)"]
+
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
+                            adaptive_ref = true, compute_error)
+
+# plot work-precision diagram
+plot(wp, labels; title = "Robertson benchmark", legend = :topright,     
+     color = permutedims([repeat([1], 3)..., repeat([3], 2)..., repeat([4], 2)...]),
+     xlims = (10^-10, 10^0), xticks = 10.0 .^ (-10:1:0),
+     ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+We see that the second and third order schemes perform very similar, with the exception of `MPRK22(0.5)`.
+This superior performance of `MPRK22(0.5)` cannot be seen in other benchmarks is therefore an exception here.
+
+The scheme `SSPMPRK22(0.5, 1.0)` has not been considered above, since it generates oscillatory solutions.
+
+```@example ROBER
+sol1 = solve(prob, SSPMPRK22(0.5, 1.0), abstol=1e-5, reltol = 1e-4);
+
+# plot solutions
+robertson_plot(sol1, ref_sol, "SSPMPRK22(0.5, 1.0)")
+```
+
+For comparisons with schemes from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) we choose the second order schemes `MPRK22(0.5)` and `MPRK22(1.0)` as well as the third order scheme `MPRK43I(0.5, 0.75)`.
+
+```@example ROBER
+sol_MPRK22_½ = solve(prob, MPRK22(0.5); abstol, reltol)
+sol_MPRK22_1 = solve(prob, MPRK22(1.0); abstol, reltol)
+sol_MPRK43 = solve(prob, MPRK43I(0.5, 0.75); abstol, reltol)
+
+p1 = robertson_plot(sol_MPRK22_½, ref_sol, "MPRK22(0.5)");
+p2 = robertson_plot(sol_MPRK22_1, ref_sol, "MPRK22(1.0)");
+p3 = robertson_plot(sol_MPRK43, ref_sol, "MPRK43I(0.5, 0.75)");
+plot(p1, p2, p3)
+```
+
+Now we compare these three schemes with a selection of second and third order stiff solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). To guarantee nonnegative solutions, we use the solver option `isoutofdomain = isnegative`.
+
+```@example ROBER
+# select reference MPRK methods
+algs1 = [MPRK22(0.5); MPRK22(1.0); MPRK43I(0.5, 0.75)]
+labels1 = ["MPRK22(0.5)"; "MPRK22(1.0)"; "MPRK43I(0.5,0.75)"]
+
+# select methods from OrdinaryDiffEq
+algs2 = [TRBDF2(); Kvaerno3(); KenCarp3(); Rodas3(); ROS2(); ROS3(); Rosenbrock23()]
+labels2 = ["TRBDF2"; "Kvearno3"; "KenCarp3"; "Rodas3"; "ROS2"; "ROS3"; "Rosenbrock23"]
+
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
+                               adaptive_ref = true, compute_error)
+# add work-precision data with isoutofdomain = isnegative                               
+work_precision_adaptive!(wp, prob, algs2, labels2, abstols, reltols, alg_ref;
+                               adaptive_ref = true, compute_error, isoutofdomain=isnegative)
+
+# plot work-precision diagram
+plot(wp, [labels1; labels2]; title = "Robertson benchmark", legend = :topright,
+     color = permutedims([repeat([1], 2)..., 3, repeat([5], 3)..., repeat([6], 4)...]),
+     xlims = (10^-10, 10^3), xticks = 10.0 .^ (-14:1:3),
+     ylims = (10^-6, 10^1), yticks = 10.0 .^ (-6:1:0), minorticks = 10)
+```
+
+We see that `MPRK22(1.0)` and `MPRK43I(0.5, 0.75)` perform similar to `Ros3()` or `Rosenbrock23()` and are a good choice as long as low accuracy is acceptable. For high accuracy we should employ a scheme like `KenCarp3()`. As for `MPRK22(0.5)` the  superior performance of `Rodas3()` seems to be an exception here.
+
+In addition,  we compare `MPRK22(1.0)` and `MPRK43I(0.5, 0.75)` to some [recommended solvers](https://docs.sciml.ai/DiffEqDocs/dev/solvers/ode_solve/) of higher order from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). Again, to guarantee positive solutions we select the solver option `isoutofdomain = isnegative`.
+
+```@example ROBER
+algs3 = [Rodas5P(); Rodas4P(); RadauIIA5()]
+labels3 = ["Rodas5P"; "Rodas4P"; "RadauIIA5"]
+
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
+                               adaptive_ref = true, compute_error)
+# add work-precision data with isoutofdomain = isnegative                                 
+work_precision_adaptive!(wp, prob, algs3, labels3, abstols, reltols, alg_ref;
+                               adaptive_ref = true, compute_error, isoutofdomain=isnegative)
+
+# plot work-precision diagram
+plot(wp, [labels1; labels3]; title = "Robertson benchmark", legend = :topright,
+     color = permutedims([repeat([1],2)..., 3, repeat([4], 2)..., 5]),
+     xlims = (10^-10, 2*10^0), xticks = 10.0 .^ (-10:1:0),
+     ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+Again, we see that the MPRK schemes are in general only beneficial if low accuracy is acceptable.
+
+### Relative maximum error over all time steps
+
+In this section we do not compare the relative maximum errors at the final time ``t = 10^{11}``, but the relative maximum errors over all time steps. 
+
+```@example ROBER
+# select relative maximum error at the end of the problem's time span.
+compute_error = rel_max_error_overall
+nothing # hide
+```
+
+First, we compare different MPRK schemes. As above, we omit `MPRK22(0.5)` and `SSPMPRK22(0.5, 1.0)`.
+
+```@example ROBER
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
+                            adaptive_ref = true, compute_error)
+
+# plot work-precision diagram
+plot(wp, labels; title = "Robertson benchmark", legend = :top,     
+     color = permutedims([repeat([1], 2)..., repeat([3], 2)..., repeat([4], 2)...]),
+     xlims = (10^-4, 5*10^1), xticks = 10.0 .^ (-5:1:2),
+     ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+Notably, the error of the second-order methods does not decrease when stricter tolerances are used. We choose the second order scheme `MPRK22(1.0)` and the third order scheme `MPRK43I(0.5, 0.75)` for comparison with solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). To guarantee nonnegative solutions of these methods, we select the solver option `isoutofdomain = isnegative`.
+
+```@example ROBER
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
+                               adaptive_ref = true, compute_error)
+# add work-precision data with isoutofdomain = isnegative                               
+work_precision_adaptive!(wp, prob, algs2, labels2, abstols, reltols, alg_ref;
+                               adaptive_ref = true, compute_error, isoutofdomain=isnegative)
+
+# plot work-precision diagram
+plot(wp, [labels1; labels2]; title = "Robertson benchmark", legend = :bottomleft,
+     color = permutedims([1, 3, repeat([5], 3)..., repeat([6], 4)...]),
+     xlims = (10^-5, 10^2), xticks = 10.0 .^ (-14:1:3),
+     ylims = (10^-6, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+Here too, some methods show that the error does not decrease even though stricter tolerances are used.
+Interestingly, the `MPRK43I(0.5, 0.75)` is superior to almost all other methods in this comparison. Only `Rodas3()` is preferable when higher accuracy is demanded.
+
+Finally, we compare `MPRK43I(0.5, 0.75)` and `MPRK22(1.0)` to [recommended solvers](https://docs.sciml.ai/DiffEqDocs/dev/solvers/ode_solve/) of higher order from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). Again, to guarantee positive solutions we select the solver option `isoutofdomain = isnegative`.
+
+```@example ROBER
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref;
+                               adaptive_ref = true, compute_error)
+# add work-precision data with isoutofdomain = isnegative                                 
+work_precision_adaptive!(wp, prob, algs3, labels3, abstols, reltols, alg_ref;
+                               adaptive_ref = true, compute_error, isoutofdomain=isnegative)
+
+# plot work-precision diagram
+plot(wp, [labels1; labels3]; title = "Robertson benchmark", legend = :topright,
+     color = permutedims([1, 3, repeat([4], 2)..., 5]),
+     xlims = (10^-4, 2*10^0), xticks = 10.0 .^ (-11:1:0),
+     ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+## Package versions
+
+These results were obtained using the following versions.
+```@example ROBER
+using InteractiveUtils
+versioninfo()
+println()
+
+using Pkg
+Pkg.status(["PositiveIntegrators", "StaticArrays", "LinearSolve", "OrdinaryDiffEq"],
+           mode=PKGMODE_MANIFEST)
+nothing # hide
+```
diff --git a/docs/src/stratospheric_reaction_benchmark.md b/docs/src/stratospheric_reaction_benchmark.md
new file mode 100644
index 00000000..c35e41c0
--- /dev/null
+++ b/docs/src/stratospheric_reaction_benchmark.md
@@ -0,0 +1,305 @@
+# [Benchmark: Solution of a stratospheric reaction problem](@id benchmark-stratos)
+
+We use the stiff stratospheric reaction problem [`prob_pds_stratreac`](@ref) to assess the efficiency of different solvers from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) and [PositiveIntegrators.jl](https://github.com/SKopecz/PositiveIntegrators.jl).
+
+```@example stratreac
+using PositiveIntegrators, OrdinaryDiffEq
+# select problem
+prob = prob_pds_stratreac
+nothing # hide
+```
+
+To keep the following code as clear as possible, we define a helper function `stratreac_plot` that we use for plotting.
+
+```@example stratreac
+using Plots
+
+function stratreac_plot(sols, labels = fill("", length(sols)), sol_ref = nothing)
+    if !(sols isa Vector)
+        sols = [sols]
+    end
+    if !(labels isa Vector)
+        labels = [labels]
+    end
+
+    tspan = prob.tspan
+    layout = (3, 2)
+    linewidth = 2
+    xticks = (range(first(tspan), last(tspan), 4), range(12.0, 84.0, 4))
+    tickfontsize = 7
+    xguide = "t [h]" #fill("t [h]", 1, 6)
+    xguidefontsize = 8
+    yguide = ["O¹ᴰ" "O" "O₃" "O₂" "NO" "NO₂"]
+    ylims = [(-20, 120) (-1e8, 7e8) (2e11, 6e11) (1.69699e16, 1.69705e16) (-2e6, 1.2e7) (1.084e9,
+                                                                                         1.098e9)]
+    legend = :outertop
+    legend_column = -1
+    widen = true
+
+    if !isnothing(sol_ref)
+        p = plot(ref_sol; layout, linestyle = :dash, label = "Ref.", linewidth)
+        for (sol, label) in zip(sols, labels)
+            plot!(p, sol; xguide, xguidefontsize, xticks, tickfontsize, yguide, legend,
+                  legend_column, widen, ylims, linewidth, label, denseplot = false)
+        end
+    else
+        p = plot(sols[1]; layout, xguide, xguidefontsize, xticks, tickfontsize, yguide,
+                 legend, legend_column, widen, ylims, linewidth,
+                 label = labels[1])
+        if length(sols) > 1
+            for (sol, label) in zip(sols[2:end], labels[2:end])
+                plot!(p, sol; layout, xguide, xguidefontsize, xticks, tickfontsize, yguide,
+                      legend, legend_column, widen, label, denseplot = false, linewidth,
+                      ylims)
+            end
+        end
+    end
+    return p
+end
+nothing # hide
+```
+First, we show approximations of `Rosenbrock23()` using loose tolerances. 
+
+```@example stratreac
+# compute reference solution for plotting
+ref_sol = solve(prob, Rodas4P(); abstol = 1e-12, reltol = 1e-11);
+
+# compute solution with low tolerances
+abstol = 1e-2
+reltol = 1e-1
+sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol);
+
+# plot solution
+stratreac_plot(sol_Ros23,  "Ros23", ref_sol)
+```
+
+Although not visible in the plots, the `Rosenbrock23` solution contains negative values.
+
+```@example stratreac
+isnonnegative(sol_Ros23)
+```
+
+Nevertheless, [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) provides the solver option `isoutofdomain`, which can be used in combination with [`isnegative`](@ref) to guarantee nonnegative solutions. 
+
+```@example stratreac
+# compute solution with isoutofdomain = isnegative
+sol_Ros23 = solve(prob, Rosenbrock23(); abstol, reltol, 
+                  isoutofdomain = isnegative); #reject negative solutions
+
+# plot solution
+stratreac_plot(sol_Ros23,  "Ros23", ref_sol)
+```
+
+For this problem the use of adaptive MPRK schemes with loose tolerances will generally result in poor approximations. In particular with respect to the O₂ component.
+
+```@example stratreac
+sol_MPRK = solve(prob, MPRK22(1.0); abstol, reltol);
+
+# plot solutions
+stratreac_plot(sol_MPRK, "MPRK22(1.0)", ref_sol)
+```
+
+To improve the solution of the MPRK scheme we can inrecase the method's `small_constant`. Trial and error has shown that `small_constant = 1e-6` is a good value for this problem and the given tolerances.
+
+```@example stratreac
+# compute MPRK solution with modified small_constant
+sol_MPRK = solve(prob, MPRK22(1.0, small_constant = 1e-6); abstol, reltol);
+
+# plot solution
+stratreac_plot(sol_MPRK, "MPRK22(1.0)", ref_sol)
+```
+
+The remaining poor approximation of the O₂ component could be due to the fact that the MPRK methods do not preserve all linear invariants, as is the case with standard methods like Runge-Kutta or Rosenbrock schemes.
+
+## Work-Precision diagrams
+
+In the following we show several work-precision diagrams, which compare different methods with respect to computing times and errors. 
+First we focus on adaptive methods, afterwards we also show results obtained with fixed time step sizes.
+
+Since the stratospheric reaction problem is stiff, we need to use a suited implicit scheme to compute its reference solution.
+
+```@example stratreac
+# select solver to compute reference solution
+alg_ref = Rodas4P()
+nothing  # hide
+```
+
+The error chosen to compare the performances of different solvers is the relative maximum error at the final time ``t = 84`` hours (``t = 302400`` seconds).
+
+```@example stratreac
+# select relative maximum error at the end of the problem's time span.
+compute_error = rel_max_error_tend
+nothing # hide
+```
+
+### Adaptive time stepping
+
+We use the functions [`work_precision_adaptive`](@ref) and [`work_precision_adaptive!`](@ref) to compute the data for the diagrams.
+Furthermore, the following absolute and relative tolerances are used.
+
+```@example stratreac
+abstols = 1.0 ./ 10.0 .^ (2:1:5)
+reltols = 10.0 .* abstols
+nothing # hide
+```
+
+We also note that using MPRK schemes with stricter tolerances, quickly requires more than a million time steps, which makes these schemes inefficient in such situations.
+
+First we compare different MPRK schemes. In addition to the default version we also use the schemes with `small_constant = 1e-6`. 
+
+```@example stratreac
+# choose methods to compare
+algs = [MPRK22(1.0); MPRK22(1.0, small_constant = 1e-6); SSPMPRK22(0.5, 1.0); SSPMPRK22(0.5, 1.0, small_constant = 1e-6);
+        MPRK43I(1.0, 0.5); MPRK43I(1.0, 0.5, small_constant = 1e-6); MPRK43I(0.5, 0.75); MPRK43I(0.5, 0.75, small_constant = 1e-6)
+        MPRK43II(0.5); MPRK43II(0.5, small_constant = 1e-6); MPRK43II(2.0 / 3.0); MPRK43II(2.0 / 3.0, small_constant = 1e-6)]
+labels = ["MPRK22(1.0)"; "MPRK22(1.0, sc=1e-6)"; "SSPMPRK22(0.5,1.0)"; "SSPMPRK22(0.5,1.0, sc=1e-6)"; 
+          "MPRK43I(1.0,0.5)"; "MPRK43I(1.0,0.5, sc=1e-6)"; "MPRK43I(0.5,0.75)"; "MPRK43I(0.5,0.75, sc=1e-6)"; "MPRK43II(0.5)"; "MPRK43II(0.5, sc=1e-6)"
+          "MPRK43II(2.0/3.0)"; "MPRK43II(2.0/3.0, sc=1e-6)"]
+
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref; compute_error)
+
+# plot work-precision diagram
+plot(wp, labels; title = "Stratospheric reaction benchmark", legend = :bottomleft,     
+     color = permutedims([repeat([1],2)..., repeat([2],2)..., repeat([3],4)..., repeat([4],4)...]),
+     xlims = (10^-7, 10^0), xticks = 10.0 .^ (-8:1:0),
+     ylims = (10^-5, 10^1), yticks = 10.0 .^ (-5:1:1), minorticks = 10)
+```
+
+We see that using `small_constant = 1e-6` clearly improves the performance of most methods. 
+For comparisons with other second and third order schemes from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) we choose the second order scheme `MPRK22(1.0, small_constant = 1e-6)` and the third order scheme `MPRK43I(0.5, 0.75)`. To guarantee positive solutions of the [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/) methods, we select the solver option `isoutofdomain = isnegative`.
+
+```@example stratreac
+# select reference MPRK methods
+algs1 = [MPRK22(1.0, small_constant = 1e-6); MPRK43I(0.5, 0.75)]
+labels1 = ["MPRK22(1.0, sc=1e-6)"; "MPRK43I(0.5,0.75)"]
+
+# select OrdinaryDiffEq methods
+algs2 = [TRBDF2(); Kvaerno3(); KenCarp3(); Rodas3(); ROS2(); ROS3(); Rosenbrock23()]
+labels2 = ["TRBDF2"; "Kvearno3"; "KenCarp3"; "Rodas3"; "ROS2"; "ROS3"; "Rosenbrock23"]
+
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref; compute_error)
+work_precision_adaptive!(wp, prob, algs2, labels2, abstols, reltols, alg_ref; compute_error, 
+                        isoutofdomain = isnegative)
+
+# plot work-precision diagram
+plot(wp, [labels1; labels2]; title = "Stratospheric reaction benchmark", legend = :bottomleft,     
+     color = permutedims([1, 3, repeat([4], 3)..., repeat([5], 4)...]),
+     xlims = (10^-7, 10^0), xticks = 10.0 .^ (-8:1:0),
+     ylims = (2*10^-4, 5*10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+Wee see that MPRK methods are advantageous if low accuracy is acceptable.
+
+In addition,  we compare `MPRK22(1.0, small_constant = 1e-6)` and  `MPRK43I(0.5, 0.75)` to some [recommended solvers](https://docs.sciml.ai/DiffEqDocs/dev/solvers/ode_solve/) of higher order from [OrdinaryDiffEq.jl](https://docs.sciml.ai/OrdinaryDiffEq/stable/). Again, to guarantee positive solutions we select the solver option `isoutofdomain = isnegative`.
+
+```@example stratreac
+# select OrdinaryDiffEq methods
+algs3 = [Rodas5P(); Rodas4P(); RadauIIA5()]
+labels3 = ["Rodas5P"; "Rodas4P"; "RadauIIA5"]
+
+# compute work-precision data
+wp = work_precision_adaptive(prob, algs1, labels1, abstols, reltols, alg_ref; compute_error)
+work_precision_adaptive!(wp, prob, algs3, labels3, abstols, reltols, alg_ref; compute_error, 
+                        isoutofdomain = isnegative)
+
+# plot work-precision diagram
+plot(wp, [labels1; labels3]; title = "Stratospheric reaction benchmark", legend = :topright,     
+     color = permutedims([1, 3, repeat([4], 3)...]),
+     xlims = (10^-7, 10^0), xticks = 10.0 .^ (-8:1:0),
+     ylims = (2*10^-4, 5*10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+Again, it can be seen that MPRK methods are only advantageous if low accuracy is acceptable.
+
+### Fixed time steps sizes
+
+Here we use fixed time step sizes instead of adaptive time stepping.
+We use the functions [`work_precision_fixed`](@ref) and [`work_precision_fixed!`](@ref) to compute the data for the diagrams.
+Please note that these functions set error and computing time to `Inf`, whenever a solution contains negative elements. 
+Consequently, such cases are not visible in the work-precision diagrams.
+
+Within the work-precision diagrams we use the following time step sizes.
+
+```@example stratreac
+# set time step sizes
+dt0 = 48 * 60 # 48 minutes
+dts = dt0 ./ 2.0 .^ (0:1:10)
+nothing # hide
+```
+
+Other than for adaptive schemes increasing `small_constant` has no positive effect, as the following examples show.
+
+```@example stratreac
+# solve prob with large step size
+sol1 = solve(prob, MPRK22(1.0); dt = dt0, adaptive = false)
+
+# plot solution
+stratreac_plot(sol1, "MPRK22(1.0)", ref_sol)
+```
+
+Choosing `small_constant = 1e-100` gives poorer results.
+
+```@example stratreac
+# solve prob with large step size and slightly increased small_constant
+sol2 = solve(prob, MPRK22(1.0, small_constant = 1e-100); dt = dt0, adaptive = false)
+
+# plot solution
+stratreac_plot(sol2, "MPRK22(1.0)", ref_sol)
+```
+
+Increasing `small_constant` only a little bit, already results in worse solutions.
+For this reason we will only consider schemes with default values of `small_constant`.
+
+```@example stratreac
+# select schemes
+algs = [MPRK22(1.0); SSPMPRK22(0.5, 1.0); MPRK43I(1.0, 0.5); MPRK43I(0.5, 0.75); MPRK43II(0.5); MPRK43II(2.0 / 3.0);
+        SSPMPRK43()]
+labels = ["MPRK22(1.0)"; "SSPMPRK22(0.5,1.0)"; "MPRK43I(1.0,0.5)"; "MPRK43I(0.5,0.75)";  "MPRK43II(0.5)"; "MPRK43II(2.0/3.0)"; 
+          "SSPMPRK43()"]
+
+# compute work-precision data
+wp = work_precision_fixed(prob, algs, labels, dts, alg_ref; compute_error)
+
+# plot work-precision diagram
+plot(wp, labels; title = "Stratospheric reaction benchmark", legend = :bottomleft,     
+     color = permutedims([1, 2, repeat([3],2)..., repeat([4],2)..., 5]),
+     xlims = (10^-6, 10^2), xticks = 10.0 .^ (-8:1:2),
+     ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:1), minorticks = 10)
+```
+
+Apart from `SSPMPRK22(0.5, 1.0)` all schemes perform quite similar. We choose `MPRK22(1.0)` and `MPRK43II(0.5)` for comparisons with other schemes.
+
+For the chosen time step sizes none of the above used standard schemes provides nonnegative solutions.
+
+```@example stratreac
+# select reference MPRK methods
+algs = [MPRK22(1.0); MPRK43II(0.5); TRBDF2(); Kvaerno3(); KenCarp3(); Rodas3(); ROS2(); ROS3(); Rosenbrock23(); 
+         Rodas5P(); Rodas4P()]
+labels = ["MPRK22(1.0)"; "MPRK43II(0.5)"; "TRBDF2"; "Kvearno3"; "KenCarp3"; "Rodas3"; "ROS2"; "ROS3"; "Rosenbrock23"; 
+          "Rodas5P"; "Rodas4P"]
+
+# compute work-precision data
+wp = work_precision_fixed(prob, algs, labels, dts, alg_ref; compute_error)
+
+# plot work-precision diagram
+plot(wp, labels; title = "Stratospheric reaction benchmark", legend = :bottomleft,     
+     color = permutedims([1, 3, repeat([4], 3)..., repeat([5], 4)..., repeat([6], 3)...]),
+     xlims = (10^-6, 10^1), xticks = 10.0 .^ (-12:2:4),
+     ylims = (10^-5, 10^0), yticks = 10.0 .^ (-5:1:0), minorticks = 10)
+```
+
+## Package versions
+
+These results were obtained using the following versions.
+```@example STRATREAC
+using InteractiveUtils
+versioninfo()
+println()
+
+using Pkg
+Pkg.status(["PositiveIntegrators", "StaticArrays", "LinearSolve", "OrdinaryDiffEq"],
+           mode=PKGMODE_MANIFEST)
+nothing # hide
+```
diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl
index c2c11ef2..81e488eb 100644
--- a/src/PositiveIntegrators.jl
+++ b/src/PositiveIntegrators.jl
@@ -2,6 +2,8 @@ module PositiveIntegrators
 
 # 1. Load dependencies
 using LinearAlgebra: LinearAlgebra, Tridiagonal, I, diag, mul!
+using Statistics: median
+
 using SparseArrays: SparseArrays, AbstractSparseMatrix,
                     issparse, nonzeros, nzrange, rowvals, spdiagm
 using StaticArrays: SVector, SMatrix, StaticArray, @SVector, @SMatrix
@@ -38,6 +40,8 @@ import OrdinaryDiffEq: alg_order, isfsal,
                        _ode_interpolant, _ode_interpolant!,
                        get_fsalfirstlast
 
+using RecipesBase: @recipe
+
 # 2. Export functionality defining the public API
 export PDSFunction, PDSProblem
 export ConservativePDSFunction, ConservativePDSProblem
@@ -50,6 +54,9 @@ export prob_pds_linmod, prob_pds_linmod_inplace, prob_pds_nonlinmod,
        prob_pds_bertolazzi, prob_pds_npzd, prob_pds_stratreac, prob_pds_minmapk
 
 export isnegative, isnonnegative
+export work_precision_adaptive, work_precision_adaptive!, work_precision_fixed,
+       work_precision_fixed!
+export rel_max_error_overall, rel_max_error_tend, rel_l1_error_tend, rel_l2_error_tend
 
 # 3. Load source code
 
diff --git a/src/utilities.jl b/src/utilities.jl
index 965af5ed..0a0fa1a2 100644
--- a/src/utilities.jl
+++ b/src/utilities.jl
@@ -30,7 +30,9 @@ end
 """
     isnegative(sol::ODESolution)
 
-Returns `true` if `sol` contains negative elements.
+Returns `true` if `sol.u` contains negative elements.
+    
+Please note that negative values may occur when plotting the solution, depending on the interpolation used.
 
 See also [`isnonnegative`](@ref).
 """
@@ -44,3 +46,259 @@ end
 Negation of [`isnegative`](@ref).    
 """
 isnonnegative(args...) = !isnegative(args...)
+
+### Errors #########################################################################
+"""
+    rel_max_error_tend(sol, ref_sol)
+
+Returns the relative maximum error between sol and ref_sol at time `sol.t[end]`.   
+"""
+function rel_max_error_tend(sol, ref_sol)
+    return maximum(abs.((sol[end] .- ref_sol[end]) ./ ref_sol[end]))
+end
+
+"""
+    rel_max_error_overall(sol, ref_sol)
+
+Returns the maximum of the relative maximum errors between sol and ref_sol over all time steps.    
+"""
+function rel_max_error_overall(sol, ref_sol)
+    err = zero(eltype(eltype(sol)))
+    for i in eachindex(sol)
+        max_err_i = maximum(abs.((abs.(sol[i]) .- abs.(ref_sol[i])) ./ ref_sol[i]))
+        if max_err_i > err
+            err = max_err_i
+        end
+    end
+    return err
+end
+
+"""
+    rel_l1_error_tend(sol, ref_sol)
+
+Returns the relative l1 error between sol and ref_sol at time `sol.t[end]`.   
+"""
+function rel_l1_error_tend(sol, ref_sol)
+    return sum(abs.((sol[end] .- ref_sol[end]) ./ ref_sol[end])) / length(ref_sol[end])
+end
+
+"""
+    rel_l2_error_tend(sol, ref_sol)
+
+Returns the relative l2 error between sol and ref_sol at time `sol.t[end]`.   
+"""
+function rel_l2_error_tend(sol, ref_sol)
+    return sqrt(sum(((sol[end] .- ref_sol[end]) ./ ref_sol[end]) .^ 2) /
+                length(ref_sol[end]))
+end
+
+### Functions to compute work-precision diagrams ########################################## 
+function _compute_time(benchmark_f, seconds, numruns)
+    benchmark_f() # pre-compile
+
+    time = benchmark_f()
+
+    if time ≤ seconds
+        time = median([time; [benchmark_f() for i in 2:numruns]])
+    end
+
+    return time
+end
+
+function compute_time_fixed(dt, prob, alg, seconds, numruns)
+    benchmark_f = let dt = dt, prob = prob, alg = alg
+        () -> @elapsed solve(prob, alg; dt, adaptive = false,
+                             save_everystep = false)
+    end
+
+    return _compute_time(benchmark_f, seconds, numruns)
+end
+
+function compute_time_adaptive(abstol, reltol, prob, alg, seconds, numruns, kwargs...)
+    benchmark_f = let abstol = abstol, reltol = reltol, prob = prob, alg = alg,
+        kwargs = kwargs
+
+        () -> @elapsed solve(prob, alg; abstol, reltol, adaptive = true,
+                             save_everystep = false, kwargs...)
+    end
+
+    return _compute_time(benchmark_f, seconds, numruns)
+end
+
+"""
+    work_precision_fixed!(dict, prob, algs, labels, dts, alg_ref;
+                          compute_error = rel_max_error_tend, 
+                          seconds = 2,
+                          numruns = 10000)
+    )
+
+Adds work-precision data to the dictionary `dict`, which was created with `work_precion_fixed`.
+See [`work_precision_fixed`](@ref) for the meaning of the inputs.
+"""
+function work_precision_fixed!(dict, prob, algs, labels, dts, alg_ref;
+                               compute_error = rel_max_error_tend, seconds = 2,
+                               numruns = 10000)
+    tspan = prob.tspan
+    dt_ref = (last(tspan) - first(tspan)) ./ 1e5
+    ref_sol = solve(prob, alg_ref; dt = dt_ref, adaptive = false, save_everystep = true)
+
+    let ref_sol = ref_sol
+        for (alg, label) in zip(algs, labels)
+            println(label)
+            error_time = Vector{Tuple{Float64, Float64}}(undef, length(dts))
+
+            for (i, dt) in enumerate(dts)
+                error_time[i] = (Inf, Inf)
+                try
+                    sol = solve(prob, alg; dt, adaptive = false, save_everystep = true)
+                    if Int(sol.retcode) == 1 && isnonnegative(sol)
+                        error = compute_error(sol, ref_sol(sol.t))
+                        time = compute_time_fixed(dt, prob, alg, seconds, numruns)
+
+                        error_time[i] = (error, time)
+                    else
+                        error_time[i] = (Inf, Inf)
+                    end
+                catch e
+                end
+            end
+            dict[label] = error_time
+        end
+    end
+end
+
+"""
+    work_precision_fixed(prob, algs, labels, dts, alg_ref;
+                         compute_error = rel_max_error_tend,
+                         seconds = 2, 
+                         numruns = 10000)
+
+Returns a dictionary to create work-precision diagrams. 
+The problem `prob` is solved by each algorithm in `algs` for all the step sizes defined in `dts`. 
+For each step size the error and computing time are stored in the dictionary. 
+If the solve is not successful for a given step size, then `(Inf, Inf)` is stored in the dictionary.
+The strings in the array `labels` are used as keys of the dictionary. 
+The reference solution used for error computations is computed with the algorithm `alg_ref`.
+
+### Keyword arguments: ###
+
+- compute_error(sol::ODESolution, ref_sol::ODESolution): Function to compute the error between `sol` and `ref_sol`.
+- seconds: If the measured computing time of a single solve is larger than `seconds`, then this computing time is stored in the dictionary.
+- numruns: If the measured computing time of a single solve is less or equal to `seconds`, then `numruns` solves are performed and the median of the respective computing times is stored in the dictionary.
+"""
+function work_precision_fixed(prob, algs, labels, dts, alg_ref;
+                              compute_error = rel_max_error_tend,
+                              seconds = 2, numruns = 10000)
+    dict = Dict(label => [] for label in labels)
+    work_precision_fixed!(dict, prob, algs, labels, dts, alg_ref; compute_error, seconds,
+                          numruns)
+    return dict
+end
+
+"""
+    work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
+                            adaptive_ref = false,
+                            abstol_ref = 1e-14, 
+                            reltol_ref = 1e-13,
+                            compute_error = rel_max_error_tend,
+                            seconds = 2, 
+                            numruns = 10000,
+                            kwargs...)
+
+Adds work-precision data to the dictionary `dict`, which was created with `work_precion_fixed_adaptive`.
+See [`work_precision_adaptive`](@ref) for the meaning of the inputs.
+"""
+function work_precision_adaptive!(dict, prob, algs, labels, abstols, reltols, alg_ref;
+                                  adaptive_ref = false,
+                                  abstol_ref = 1e-14, reltol_ref = 1e-13,
+                                  compute_error = rel_max_error_tend,
+                                  seconds = 2, numruns = 10000, kwargs...)
+    if adaptive_ref
+        ref_sol = solve(prob, alg_ref; adaptive = true, save_everystep = true,
+                        abstol = abstol_ref, reltol = reltol_ref)
+    else
+        tspan = prob.tspan
+        dt_ref = (last(tspan) - first(tspan)) ./ 1e5
+        ref_sol = solve(prob, alg_ref; dt = dt_ref, adaptive = false, save_everystep = true)
+    end
+
+    for (alg, label) in zip(algs, labels)
+        println(label)
+        error_time = Vector{Tuple{Float64, Float64}}(undef, length(abstols))
+
+        for (i, dt) in enumerate(abstols)
+            abstol = abstols[i]
+            reltol = reltols[i]
+            sol = solve(prob, alg; abstol, reltol, save_everystep = true,
+                        kwargs...)
+
+            if Int(sol.retcode) == 1 && isnonnegative(sol)
+                error = compute_error(sol, ref_sol(sol.t))
+                time = compute_time_adaptive(abstol, reltol, prob, alg, seconds, numruns,
+                                             kwargs...)
+
+                error_time[i] = (error, time)
+            else
+                error_time[i] = (Inf, Inf)
+            end
+        end
+        dict[label] = error_time
+    end
+
+    return nothing
+end
+
+"""
+    work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
+                            adaptive_ref = false,
+                            abstol_ref = 1e-14, 
+                            reltol_ref = 1e-13,
+                            compute_error = rel_max_error_tend,
+                            seconds = 2, 
+                            numruns = 10000, 
+                            kwargs...)
+
+Returns a dictionary to create work-precision diagrams. 
+The problem `prob` is solved by each algorithm in `algs` for all tolerances defined in `abstols` and `reltols`. 
+For the respective tolerances the error and computing time are stored in the dictionary. 
+If the solve is not successful for the given tolerances, then `(Inf, Inf)` is stored in the dictionary.
+The strings in the array `labels` are used as keys of the dictionary. 
+The reference solution used for error computations is computed with the algorithm `alg_ref`.
+Additional keyword arguments are passed on to `solve`.
+
+### Keyword arguments: ###
+
+- `adaptive_ref`: If `true` the refenerce solution is computed adaptively with tolerances `abstol_ref` and `reltol_ref`. Otherwise ``10^5`` steps are used.
+- `abstol_ref`: See `adaptive_ref`.
+- `reltol_ref`: See `adaptive_ref`.
+- `compute_error(sol::ODESolution, ref_sol::ODESolution)`: A function to compute the error between `sol` and `ref_sol`.
+- `seconds`: If the measured computing time of a single solve is larger than `seconds`, then this computing time is stored in the dictionary.
+- `numruns`: If the measured computing time of a single solve is less or equal to `seconds`, then `numruns` solves are performed and the median of the respective computing times is stored in the dictionary.
+"""
+function work_precision_adaptive(prob, algs, labels, abstols, reltols, alg_ref;
+                                 adaptive_ref = false,
+                                 abstol_ref = 1e-14, reltol_ref = 1e-13,
+                                 compute_error = rel_max_error_tend, seconds = 2,
+                                 numruns = 10000, kwargs...)
+    dict = Dict(label => [] for label in labels)
+    work_precision_adaptive!(dict, prob, algs, labels, abstols, reltols, alg_ref;
+                             adaptive_ref, abstol_ref, reltol_ref,
+                             compute_error, seconds,
+                             numruns, kwargs...)
+    return dict
+end
+
+# plot recipe to plot work-precision dictionaries
+@recipe function f(wp::Dict, labels; color = nothing)
+    seriestype --> :path
+    linewidth --> 3
+    xscale --> :log10
+    yscale --> :log10
+    markershape --> :auto
+    xs = [first.(wp[label]) for label in labels]
+    ys = [last.(wp[label]) for label in labels]
+    xguide --> "error"
+    yguide --> "times (s)"
+    label --> reshape(labels, 1, length(wp))
+    return xs, ys
+end
diff --git a/test/Project.toml b/test/Project.toml
index 5182e22f..968f25b2 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -4,6 +4,7 @@ ExplicitImports = "7d51a73a-1435-4ff3-83d9-f097790105c7"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -17,6 +18,7 @@ ExplicitImports = "1.0.1"
 LinearSolve = "2.21"
 OrdinaryDiffEq = "6.59"
 Plots = "1.11.1"
+RecipesBase = "1.3"
 StaticArrays = "1.5"
 Statistics = "1"
 Unitful = "1.21"
diff --git a/test/runtests.jl b/test/runtests.jl
index 3c0b1829..5b20d345 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,7 +1,7 @@
 using Test
 using LinearAlgebra
 using SparseArrays
-using Statistics: mean
+using Statistics: mean, median
 
 using StaticArrays: MVector, @SVector, SA
 
@@ -13,6 +13,7 @@ using PositiveIntegrators
 using LinearSolve: RFLUFactorization, LUFactorization, KrylovJL_GMRES
 
 using Aqua: Aqua
+using RecipesBase: RecipesBase # only for Aqua tests
 using ExplicitImports: check_no_implicit_imports, check_no_stale_explicit_imports
 
 """
@@ -210,7 +211,8 @@ end
         # We do not test ambiguities since we get a lot of
         # false positives from dependencies
         Aqua.test_all(PositiveIntegrators;
-                      ambiguities = false,)
+                      ambiguities = false,
+                      piracies = (; treat_as_own = [RecipesBase.apply_recipe],),)
     end
 
     @testset "ExplicitImports.jl" begin
@@ -2215,8 +2217,19 @@ end
 
     @testset "plot" begin
         using Plots
+
+        # plot ode solution
         sol = solve(prob_pds_linmod, MPRK22(1.0))
         @test_nowarn plot(sol)
+
+        # plot work-precision diagram
+        prob = prob_pds_nonlinmod
+        algs = [MPE(); MPRK22(1.0)]
+        labels = ["MPE"; "MPRK"]
+        dts = [0.1; 0.01; 0.001]
+        alg_ref = Vern7()
+        wp = work_precision_fixed(prob, algs, labels, dts, alg_ref)
+        @test_nowarn plot(wp, labels; minorticks = 10)
     end
 
     @testset "utilities" begin
@@ -2235,5 +2248,82 @@ end
             @test isnegative(sol_bruss.u)
             @test isnegative(last(sol_bruss.u))
         end
+
+        @testset "errors" begin
+            v = [[1.0; 2.0; 3.0], [4.0; 5.0; 6.0]]
+            w = [[8.0; 10.0; 12.0], [2.0; 4.0; 6.0]]
+
+            @test rel_max_error_tend(v, w) == 1.0
+            @test rel_l1_error_tend(v, w) == 1.25 / 3
+            @test rel_l2_error_tend(v, w) == sqrt(1.0625 / 3)
+            @test rel_max_error_overall(w, v) == 7.0
+        end
+
+        # Here we run a single scheme multiple times and check
+        # that error agree and computing times are close enough
+        @testset "work-precision fixed" begin
+            prob = prob_pds_nonlinmod
+            alg = MPRK22(1.0)
+            algs = [alg; alg; alg; alg; alg]
+            labelsA = ["A_1"; "A_2"; "A_3"; "A_4"; "A_5"]
+            labelsB = ["B_1"; "B_2"; "B_3"; "B_4"; "B_5"]
+            dts = (last(prob.tspan) - first(prob.tspan)) / 10.0 * 0.5 .^ (4:13)
+            alg_ref = Vern7()
+            wp = work_precision_fixed(prob, [alg; alg; alg; alg; alg], labelsA, dts,
+                                      alg_ref)
+            work_precision_fixed!(wp, prob, [alg; alg; alg; alg; alg], labelsB, dts,
+                                  alg_ref)
+
+            # check that errors agree
+            for (i, _) in enumerate(dts)
+                v = [value[i][1] for (key, value) in wp]
+                @test all(y -> y == v[1], v)
+            end
+
+            # check that computing times are close enough 
+            for (i, _) in enumerate(dts)
+                @show i
+                v = [value[i][2] for (key, value) in wp]
+                m1 = mean(v)
+                m2 = median(v)
+
+                @test maximum((v .- m1) ./ m1) < 0.3
+                @test maximum((v .- m2) ./ m2) < 0.3
+            end
+        end
+
+        # Here we run a single scheme multiple times and check
+        # that error agree and computing times are close enough
+        @testset "work-precision adaptive" begin
+            prob = prob_pds_nonlinmod
+            alg = MPRK22(1.0)
+            algs = [alg; alg; alg; alg; alg]
+            labelsA = ["A_1"; "A_2"; "A_3"; "A_4"; "A_5"]
+            labelsB = ["B_1"; "B_2"; "B_3"; "B_4"; "B_5"]
+            abstols = 1 ./ 10 .^ (4:8)
+            reltols = 1 ./ 10 .^ (3:7)
+            alg_ref = Vern7()
+            wp = work_precision_adaptive(prob, [alg; alg; alg; alg; alg], labelsA, abstols,
+                                         reltols, alg_ref)
+            work_precision_adaptive!(wp, prob, [alg; alg; alg; alg; alg], labelsB, abstols,
+                                     reltols, alg_ref)
+
+            # check that errors agree
+            for (i, _) in enumerate(abstols)
+                v = [value[i][1] for (key, value) in wp]
+                @test all(y -> y == v[1], v)
+            end
+
+            # check that computing times are close enough 
+            for (i, _) in enumerate(abstols)
+                @show i
+                v = [value[i][2] for (key, value) in wp]
+                m1 = mean(v)
+                m2 = median(v)
+
+                @test maximum((v .- m1) ./ m1) < 0.3
+                @test maximum((v .- m2) ./ m2) < 0.3
+            end
+        end
     end
 end;