diff --git a/Project.toml b/Project.toml
index 1de121ee..c4b2d792 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,7 +10,6 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 SimpleUnPack = "ce78b400-467f-4804-87d8-8f486da07d0a"
@@ -29,4 +28,6 @@ Reexport = "1"
 SciMLBase = "2"
 SimpleUnPack = "1"
 SparseArrays = "1.7"
+StaticArrays = "1.5"
+SymbolicIndexingInterface = "0.2, 0.3"
 julia = "1.7"
diff --git a/examples/04_example_problemlibrary.jl b/examples/04_example_problemlibrary.jl
index b2c0e1ae..42a01e88 100644
--- a/examples/04_example_problemlibrary.jl
+++ b/examples/04_example_problemlibrary.jl
@@ -49,7 +49,7 @@ plot(sol_nonlinmod, legend = :right)
 myplot!(sol_nonlinmod_MPE, "MPE")
 plot!(sol_nonlinmod_MPE, idxs = (f3, 0, 1, 2, 3))
 
-# convergence order 
+# convergence order
 test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
 convergence_tab_plot(prob_pds_nonlinmod, [MPE(), Euler()], test_setup; dts = 0.5 .^ (3:17),
                      order_plot = true)
@@ -74,7 +74,7 @@ myplot!(sol_brusselator_MPE, "MPE")
 plot!(sol_brusselator_MPE, idxs = (f_brusselator, 0, 1, 2, 3, 4, 5, 6),
       label = "f_brusselator")
 
-# convergence order 
+# convergence order
 test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
 convergence_tab_plot(prob_pds_brusselator, [MPE()], test_setup; dts = 0.5 .^ (3:17),
                      order_plot = true)
@@ -93,7 +93,7 @@ myplot!(sol_sir_Euler, "Euler")
 plot!(sol_sir_Euler, idxs = (f3, 0, 1, 2, 3), label = "f3")
 plot(p1, p2)
 
-# convergence order 
+# convergence order
 test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
 convergence_tab_plot(prob_pds_sir, [MPE(), Euler()], test_setup; dts = 0.5 .^ (1:15),
                      order_plot = true)
@@ -102,13 +102,13 @@ convergence_tab_plot(prob_pds_sir, [MPE(), Euler()], test_setup; dts = 0.5 .^ (1
 sol_bertolazzi = solve(prob_pds_bertolazzi, TRBDF2());
 sol_bertolazzi_MPE = solve(prob_pds_bertolazzi, MPE(), dt = 0.01);
 
-# plot 
+# plot
 plot(sol_bertolazzi, legend = :right)
 myplot!(sol_bertolazzi_MPE, "MPE")
 ylims!((-0.5, 3.5))
 plot!(sol_bertolazzi_MPE, idxs = (f3, 0, 1, 2, 3))
 
-# convergence order 
+# convergence order
 test_setup = Dict(:alg => Rosenbrock23(), :reltol => 1e-8, :abstol => 1e-8)
 convergence_tab_plot(prob_pds_bertolazzi, [MPE(), ImplicitEuler()], test_setup;
                      dts = 0.5 .^ (10:15), order_plot = true)
@@ -123,13 +123,13 @@ myplot!(sol_npzd_MPE, "MPE")
 plot!(sol_npzd_MPE, idxs = (f_npzd, 0, 1, 2, 3, 4), label = "f_npzd")
 plot!(legend = :bottomright)
 
-# convergence order 
-# error should take all time steps into account, not only the final time! 
+# convergence order
+# error should take all time steps into account, not only the final time!
 test_setup = Dict(:alg => Rosenbrock23(), :reltol => 1e-14, :abstol => 1e-14)
 convergence_tab_plot(prob_pds_npzd, [MPE(), ImplicitEuler()], test_setup;
                      dts = 0.5 .^ (5:17), order_plot = true)
 
-### strat reac problem ####################################################
+### stratospheric reaction problem ####################################################
 sol_stratreac = solve(prob_pds_stratreac, TRBDF2(autodiff = false));
 #sol_stratreac_MPE = solve(prob_pds_stratreac, MPE(), dt = 25920.0);
 #tmp = reduce(hcat, sol_stratreac_MPE.u)
diff --git a/examples/Project.toml b/examples/Project.toml
index 51556992..c173b243 100644
--- a/examples/Project.toml
+++ b/examples/Project.toml
@@ -8,6 +8,7 @@ Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 PositiveIntegrators = "d1b20bf0-b083-4985-a874-dc5121669aa5"
 PrettyTables = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 
 [compat]
diff --git a/src/PDSProblemLibrary.jl b/src/PDSProblemLibrary.jl
index cd396311..d9120c11 100644
--- a/src/PDSProblemLibrary.jl
+++ b/src/PDSProblemLibrary.jl
@@ -19,7 +19,7 @@ end
 u0_nonlinmod = @SVector [9.98; 0.01; 0.01]
 prob_pds_nonlinmod = ConservativePDSProblem(P_nonlinmod, u0_nonlinmod, (0.0, 30.0))
 
-# robertson problem 
+# robertson problem
 function P_robertson(u, p, t)
     @SMatrix [0.0 1e4*u[2]*u[3] 0.0; 4e-2*u[1] 0.0 0.0; 0.0 3e7*u[2]^2 0.0]
 end
@@ -40,7 +40,7 @@ end;
 u0_brusselator = @SVector [10.0, 10.0, 0.0, 0.0, 0.1, 0.1]
 prob_pds_brusselator = ConservativePDSProblem(P_brusselator, u0_brusselator, (0.0, 10.0))
 
-# SIR problem 
+# SIR problem
 P_sir(u, p, t) = @SMatrix [0.0 0.0 0.0; 2*u[1]*u[2] 0.0 0.0; 0.0 u[2] 0.0]
 u0_sir = @SVector [0.99, 0.005, 0.005]
 prob_pds_sir = ConservativePDSProblem(P_sir, u0_sir, (0.0, 20.0))
@@ -72,7 +72,7 @@ end
 u0_npzd = @SVector [8.0, 2.0, 1.0, 4.0]
 prob_pds_npzd = ConservativePDSProblem(P_npzd, u0_npzd, (0.0, 10.0))
 
-#strat reac problem
+# stratospheric reaction problem
 function P_stratreac(u, p, t)
     O1D, O, O3, O2, NO, NO2 = u