diff --git a/Project.toml b/Project.toml
index fdc71584..32939ff8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "PositiveIntegrators"
 uuid = "d1b20bf0-b083-4985-a874-dc5121669aa5"
 authors = ["Stefan Kopecz, Hendrik Ranocha, and contributors"]
-version = "0.2.0"
+version = "0.2.1"
 
 [deps]
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
diff --git a/docs/src/index.md b/docs/src/index.md
index 22d89feb..5e0a8aed 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -174,7 +174,7 @@ Finally, we can use [Plots.jl](https://docs.juliaplots.org/stable/) to visualize
 using Plots
 
 plot(sol, label = ["S" "I" "R"], legend=:right)
-plot!(sol, idxs = ((t, S, I, R) -> (t, S + I + R), 0, 1, 2, 3), label = "S+I+R") #Plot S+I+R over time.
+plot!(sol.t, sum.(sol.u), label = "S+I+R") # Plot S+I+R over time.
 ```
 We see that there is always a nonnegative number of people in each compartment, while the population ``S+I+R`` remains constant over time.
 
diff --git a/src/PositiveIntegrators.jl b/src/PositiveIntegrators.jl
index dc87fb75..87c3f589 100644
--- a/src/PositiveIntegrators.jl
+++ b/src/PositiveIntegrators.jl
@@ -40,7 +40,8 @@ using OrdinaryDiffEq: @cache,
                       recursivefill!, _vec, wrapprecs, dolinsolve
 import OrdinaryDiffEq: alg_order, isfsal,
                        calculate_residuals, calculate_residuals!,
-                       alg_cache, initialize!, perform_step!,
+                       alg_cache, get_tmp_cache,
+                       initialize!, perform_step!,
                        _ode_interpolant, _ode_interpolant!
 
 # 2. Export functionality defining the public API
diff --git a/src/mprk.jl b/src/mprk.jl
index 6e54e4b2..83f11491 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -337,6 +337,8 @@ struct MPEConservativeCache{PType, uType, tabType, F} <: OrdinaryDiffEqMutableCa
     linsolve::F
 end
 
+get_tmp_cache(integrator, ::MPE, cache::OrdinaryDiffEqMutableCache) = (cache.σ,)
+
 # In-place
 function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
@@ -638,6 +640,8 @@ struct MPRK22ConservativeCache{uType, PType, tabType, F} <:
     linsolve::F
 end
 
+get_tmp_cache(integrator, ::MPRK22, cache::OrdinaryDiffEqMutableCache) = (cache.σ,)
+
 # In-place
 function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
@@ -1218,6 +1222,11 @@ struct MPRK43ConservativeCache{uType, PType, tabType, F} <: OrdinaryDiffEqMutabl
     linsolve::F
 end
 
+function get_tmp_cache(integrator, ::Union{MPRK43I, MPRK43II},
+                       cache::OrdinaryDiffEqMutableCache)
+    (cache.σ,)
+end
+
 # In-place
 function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index c25c79a3..5473033b 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -223,6 +223,8 @@ struct SSPMPRK22ConservativeCache{uType, PType, tabType, F} <:
     linsolve::F
 end
 
+get_tmp_cache(integrator, ::SSPMPRK22, cache::OrdinaryDiffEqMutableCache) = (cache.σ,)
+
 # In-place
 function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
@@ -734,6 +736,8 @@ struct SSPMPRK43ConservativeCache{uType, PType, tabType, F} <: OrdinaryDiffEqMut
     linsolve::F
 end
 
+get_tmp_cache(integrator, ::SSPMPRK43, cache::OrdinaryDiffEqMutableCache) = (cache.σ,)
+
 # In-place
 function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
diff --git a/test/runtests.jl b/test/runtests.jl
index 0148c6d6..1907244e 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1198,7 +1198,7 @@ end
             end
         end
 
-        # Here we check that the type of p_prototype actually 
+        # Here we check that the type of p_prototype actually
         # defines the types of the Ps inside the algorithm caches.
         # We test sparse, tridiagonal, and dense matrices.
         @testset "Prototype type check" begin
@@ -1246,7 +1246,7 @@ end
                                                 p_prototype = P_dense)
             prob_sparse = ConservativePDSProblem(prod_sparse!, u0, tspan;
                                                  p_prototype = P_sparse)
-            ## nonconservative PDS                                         
+            ## nonconservative PDS
             prob_default2 = PDSProblem(prod_dense!, dest!, u0, tspan)
             prob_tridiagonal2 = PDSProblem(prod_tridiagonal!, dest!, u0, tspan;
                                            p_prototype = P_tridiagonal)
@@ -1262,7 +1262,19 @@ end
                 for prob in (prob_default, prob_tridiagonal, prob_dense, prob_sparse,
                              prob_default2,
                              prob_tridiagonal2, prob_dense2, prob_sparse2)
-                    solve(prob, alg; dt, adaptive = false)
+                    sol1 = solve(prob, alg; dt, adaptive = false)
+
+                    # test get_tmp_cache and integrator interface - modifying
+                    # values from the cache should not change the final results
+                    integrator = init(prob, alg; dt, adaptive = false)
+                    step!(integrator)
+                    cache = @inferred get_tmp_cache(integrator)
+                    @test !isempty(cache)
+                    tmp = first(cache)
+                    fill!(tmp, NaN)
+                    sol2 = solve!(integrator)
+                    @test sol1.t ≈ sol2.t
+                    @test sol1.u ≈ sol2.u
                 end
             end
         end