diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
index e83089e464..eb4cc169b9 100644
--- a/src/OrdinaryDiffEq.jl
+++ b/src/OrdinaryDiffEq.jl
@@ -389,7 +389,7 @@ export MagnusMidpoint, LinearExponential, MagnusLeapfrog, LieEuler, CayleyEuler,
 
 export Rosenbrock23, Rosenbrock32, RosShamp4, Veldd4, Velds4, GRK4T, GRK4A,
        Ros4LStab, ROS3P, Rodas3, Rodas23W, Rodas3P, Rodas4, Rodas42, Rodas4P, Rodas4P2,
-       Rodas5, Rodas5P,
+       Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr,
        RosenbrockW6S4OS, ROS34PW1a, ROS34PW1b, ROS34PW2, ROS34PW3, ROS34PRw, ROS3PRL,
        ROS3PRL2,
        ROS2, ROS2PR, ROS2S, ROS3, ROS3PR, Scholz4_7
diff --git a/src/alg_utils.jl b/src/alg_utils.jl
index 40a31190a8..62ad86e2af 100644
--- a/src/alg_utils.jl
+++ b/src/alg_utils.jl
@@ -34,6 +34,8 @@ isfsal(alg::Rodas3P) = false
 isfsal(alg::Rodas23W) = false
 isfsal(alg::Rodas5) = false
 isfsal(alg::Rodas5P) = false
+isfsal(alg::Rodas5Pr) = false
+isfsal(alg::Rodas5Pe) = false
 isfsal(alg::Rodas4) = false
 isfsal(alg::Rodas42) = false
 isfsal(alg::Rodas4P) = false
@@ -664,6 +666,8 @@ alg_order(alg::Rodas4P) = 4
 alg_order(alg::Rodas4P2) = 4
 alg_order(alg::Rodas5) = 5
 alg_order(alg::Rodas5P) = 5
+alg_order(alg::Rodas5Pr) = 5
+alg_order(alg::Rodas5Pe) = 5
 
 alg_order(alg::AB3) = 3
 alg_order(alg::AB4) = 4
@@ -1037,6 +1041,7 @@ isstandard(alg::VCABM) = true
 isWmethod(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = false
 isWmethod(alg::Rosenbrock23) = true
 isWmethod(alg::Rosenbrock32) = true
+isWmethod(alg::Rodas23W) = true
 isWmethod(alg::ROS2S) = true
 isWmethod(alg::ROS34PW1a) = true
 isWmethod(alg::ROS34PW1b) = true
diff --git a/src/algorithms.jl b/src/algorithms.jl
index a60cbb2f90..b8ff73dee7 100644
--- a/src/algorithms.jl
+++ b/src/algorithms.jl
@@ -2967,6 +2967,11 @@ University of Geneva, Switzerland.
 - Steinebach G.   Construction of Rosenbrock–Wanner method Rodas5P and numerical benchmarks within the Julia Differential Equations package. 
  In: BIT Numerical Mathematics, 63(2), 2023
 
+ #### Rodas23W, Rodas3P, Rodas5Pe, Rodas5Pr
+- Steinebach G. Rosenbrock methods within OrdinaryDiffEq.jl - Overview, recent developments and applications - 
+ Preprint 2024
+ https://github.com/hbrs-cse/RosenbrockMethods/blob/main/paper/JuliaPaper.pdf
+
 =#
 
 for Alg in [
@@ -3000,7 +3005,9 @@ for Alg in [
     :Rodas4P,
     :Rodas4P2,
     :Rodas5,
-    :Rodas5P
+    :Rodas5P,
+    :Rodas5Pe,
+    :Rodas5Pr
 ]
     @eval begin
         struct $Alg{CS, AD, F, P, FDT, ST, CJ} <:
diff --git a/src/caches/rosenbrock_caches.jl b/src/caches/rosenbrock_caches.jl
index 29eef0f9d5..7aeac95ebe 100644
--- a/src/caches/rosenbrock_caches.jl
+++ b/src/caches/rosenbrock_caches.jl
@@ -1045,7 +1045,7 @@ function alg_cache(alg::Rodas5, u, rate_prototype, ::Type{uEltypeNoUnits},
             constvalue(tTypeNoUnits)), J, W, linsolve)
 end
 
-function alg_cache(alg::Rodas5P, u, rate_prototype, ::Type{uEltypeNoUnits},
+function alg_cache(alg::Union{Rodas5P, Rodas5Pe, Rodas5Pr}, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
         dt, reltol, p, calck,
         ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
@@ -1093,7 +1093,7 @@ function alg_cache(alg::Rodas5P, u, rate_prototype, ::Type{uEltypeNoUnits},
         linsolve, jac_config, grad_config, reltol, alg)
 end
 
-function alg_cache(alg::Rodas5P, u, rate_prototype, ::Type{uEltypeNoUnits},
+function alg_cache(alg::Union{Rodas5P, Rodas5Pe, Rodas5Pr}, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
         dt, reltol, p, calck,
         ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
diff --git a/src/dense/stiff_addsteps.jl b/src/dense/stiff_addsteps.jl
index c75a06f81a..b90afc8182 100644
--- a/src/dense/stiff_addsteps.jl
+++ b/src/dense/stiff_addsteps.jl
@@ -350,7 +350,7 @@ function _ode_addsteps!(k, t, uprev, u, dt, f, p,
         always_calc_begin = false, allow_calc_end = true,
         force_calc_end = false)
     if length(k) < 2 || always_calc_begin
-        @unpack du, du1, du2, tmp, k1, k2, k3, k4, k5, k6, dT, J, W, uf, tf, linsolve_tmp, jac_config, fsalfirst = cache
+        @unpack du, du1, du2, tmp, k1, k2, k3, k4, k5, dT, J, W, uf, tf, linsolve_tmp, jac_config, fsalfirst = cache
         @unpack a21, a41, a42, a43, C21, C31, C32, C41, C42, C43, C51, C52, C53, C54, gamma, c2, c3, d1, d2, d3 = cache.tab
 
         # Assignments
diff --git a/src/derivative_utils.jl b/src/derivative_utils.jl
index 55f8ec8e17..272482dfc2 100644
--- a/src/derivative_utils.jl
+++ b/src/derivative_utils.jl
@@ -750,7 +750,7 @@ end
         if J isa StaticArray &&
            integrator.alg isa
            Union{
-            Rosenbrock23, Rodas23W, Rodas3P, Rodas4, Rodas4P, Rodas4P2, Rodas5, Rodas5P}
+            Rosenbrock23, Rodas23W, Rodas3P, Rodas4, Rodas4P, Rodas4P2, Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr}
             W = W_transform ? J - mass_matrix * inv(dtgamma) :
                 dtgamma * J - mass_matrix
         else
@@ -775,7 +775,7 @@ end
                 W_full
             elseif len !== nothing &&
                    integrator.alg isa
-                   Union{Rosenbrock23, Rodas4, Rodas4P, Rodas4P2, Rodas5, Rodas5P}
+                   Union{Rosenbrock23, Rodas23W, Rodas3P, Rodas4, Rodas4P, Rodas4P2, Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr}
                 StaticWOperator(W_full)
             else
                 DiffEqBase.default_factorize(W_full)
@@ -923,7 +923,7 @@ function build_J_W(alg, u, uprev, p, t, dt, f::F, ::Type{uEltypeNoUnits},
             len = StaticArrayInterface.known_length(typeof(J))
             if len !== nothing &&
                alg isa
-               Union{Rosenbrock23, Rodas4, Rodas4P, Rodas4P2, Rodas5, Rodas5P}
+               Union{Rosenbrock23, Rodas23W, Rodas3P, Rodas4, Rodas4P, Rodas4P2, Rodas5, Rodas5P, Rodas5Pe, Rodas5Pr}
                 StaticWOperator(J, false)
             else
                 ArrayInterface.lu_instance(J)
diff --git a/src/integrators/integrator_interface.jl b/src/integrators/integrator_interface.jl
index 92e2d797e1..b5df01e109 100644
--- a/src/integrators/integrator_interface.jl
+++ b/src/integrators/integrator_interface.jl
@@ -97,7 +97,7 @@ end
                   !(integrator.alg isa
                     Union{Rosenbrock23, Rosenbrock32, Rodas23W,
             Rodas3P, Rodas4, Rodas4P, Rodas4P2, Rodas5,
-            Rodas5P})
+            Rodas5P, Rodas5Pe, Rodas5Pr})
             # Special stiff interpolations do not store the right value in fsallast
             out .= integrator.fsallast
         else
diff --git a/src/perform_step/rosenbrock_perform_step.jl b/src/perform_step/rosenbrock_perform_step.jl
index bcbfc88865..4937407918 100644
--- a/src/perform_step/rosenbrock_perform_step.jl
+++ b/src/perform_step/rosenbrock_perform_step.jl
@@ -1068,21 +1068,26 @@ end
         integrator.k[2] = h31 * k1 + h32 * k2 + h33 * k3 + h34 * k4 + h35 * k5
         integrator.k[3] = h2_21 * k1 + h2_22 * k2 + h2_23 * k3 + h2_24 * k4 + h2_25 * k5
         if integrator.opts.adaptive
-            calculate_interpoldiff!(
-                k1, k2, uprev, du, u, integrator.k[1], integrator.k[2], integrator.k[3])
-            atmp = calculate_residuals(k2, uprev, k1, integrator.opts.abstol,
+            if  isa(linsolve_tmp,AbstractFloat)
+                u_int, u_diff = calculate_interpoldiff(uprev, du, u, integrator.k[1], integrator.k[2], integrator.k[3])
+            else
+                u_int = linsolve_tmp 
+                u_diff = linsolve_tmp .+ 0
+                calculate_interpoldiff!(u_int, u_diff, uprev, du, u, integrator.k[1], integrator.k[2], integrator.k[3])
+            end
+            atmp = calculate_residuals(u_diff, uprev, u_int, integrator.opts.abstol,
                 integrator.opts.reltol, integrator.opts.internalnorm, t)
-            EEst = max(EEst, integrator.opts.internalnorm(atmp, t))  #-- role of t unclear
+            EEst = max(EEst,integrator.opts.internalnorm(atmp, t))  #-- role of t unclear
         end
     end
 
     if (integrator.alg isa Rodas23W)
-        k1[:] = u[:]
-        u[:] = du[:]
-        du[:] = k1[:]
+        k1 = u .+ 0
+        u = du .+ 0
+        du = k1 .+ 0
         if integrator.opts.calck
-            integrator.k[1][:] = integrator.k[3][:]
-            integrator.k[2][:] .= 0.0
+            integrator.k[1] = integrator.k[3] .+ 0
+            integrator.k[2] = 0*integrator.k[2]
         end
     end
 
@@ -1432,7 +1437,6 @@ end
             calculate_residuals!(atmp, du2, uprev, du1, integrator.opts.abstol,
                 integrator.opts.reltol, integrator.opts.internalnorm, t)
             EEst = max(EEst, integrator.opts.internalnorm(atmp, t))  #-- role of t unclear
-            #println(t," ",EEst," ",du2)
         end
     end
 
@@ -1451,11 +1455,40 @@ end
         calculate_residuals!(atmp, u - du, uprev, u, integrator.opts.abstol,
             integrator.opts.reltol, integrator.opts.internalnorm, t)
         integrator.EEst = max(EEst, integrator.opts.internalnorm(atmp, t))
-        #println(t," ",EEst," ",integrator.EEst)
     end
     cache.linsolve = linres.cache
 end
 
+function calculate_interpoldiff(uprev, up2, up3, c_koeff, d_koeff, c2_koeff)
+    u_int = 0.0
+    u_diff = 0.0
+    a1 = up3 + c_koeff - up2 - c2_koeff
+    a2 = d_koeff - c_koeff + c2_koeff
+    a3 = -d_koeff
+    dis = a2^2 - 3*a1*a3
+    u_int = up3
+    u_diff = 0.0
+    if dis > 0.0 #-- Min/Max occurs
+        tau1 = (-a2 - sqrt(dis))/(3*a3)
+        tau2 = (-a2 + sqrt(dis))/(3*a3)
+        if tau1 > tau2 
+            tau1,tau2 = tau2,tau1 
+        end
+        for tau in (tau1,tau2)
+            if (tau > 0.0) && (tau < 1.0)
+                y_tau = (1 - tau)*uprev + 
+                        tau*(up3 + (1 - tau)*(c_koeff + tau*d_koeff))
+                dy_tau = ((a3*tau + a2)*tau + a1)*tau
+                if abs(dy_tau) > abs(u_diff)
+                    u_diff = dy_tau
+                    u_int = y_tau
+                 end
+            end
+        end
+    end
+    return u_int, u_diff
+end
+
 function calculate_interpoldiff!(u_int, u_diff, uprev, up2, up3, c_koeff, d_koeff, c2_koeff)
     for i in eachindex(up2)
         a1 = up3[i] + c_koeff[i] - up2[i] - c2_koeff[i]
@@ -2179,9 +2212,14 @@ end
     k8 = _reshape(W \ -_vec(linsolve_tmp), axes(uprev))
     integrator.stats.nsolve += 1
     u = u + k8
+    linsolve_tmp = k8
 
     if integrator.opts.adaptive
-        atmp = calculate_residuals(k8, uprev, u, integrator.opts.abstol,
+	    if (integrator.alg isa Rodas5Pe)
+            linsolve_tmp = 0.2606326497975715*k1 - 0.005158627295444251*k2 + 1.3038988631109731*k3 + 1.235000722062074*k4 +
+               - 0.7931985603795049*k5 - 1.005448461135913*k6 - 0.18044626132120234*k7 + 0.17051519239113755*k8
+	    end
+        atmp = calculate_residuals(linsolve_tmp, uprev, u, integrator.opts.abstol,
             integrator.opts.reltol, integrator.opts.internalnorm, t)
         integrator.EEst = integrator.opts.internalnorm(atmp, t)
     end
@@ -2194,6 +2232,19 @@ end
                           h37 * k7 + h38 * k8
         integrator.k[3] = h41 * k1 + h42 * k2 + h43 * k3 + h44 * k4 + h45 * k5 + h46 * k6 +
                           h47 * k7 + h48 * k8
+        if (integrator.alg isa Rodas5Pr) && integrator.opts.adaptive && (integrator.EEst < 1.0)
+            k2 = 0.5*(uprev + u + 0.5 * (integrator.k[1] + 0.5 * (integrator.k[2] + 0.5 * integrator.k[3])))
+            du1 = ( 0.25*(integrator.k[2] + integrator.k[3]) - uprev + u) / dt
+            du = f(k2, p, t + dt/2)
+            integrator.stats.nf += 1
+            if mass_matrix === I
+                du2 = du1 - du
+            else
+                du2 = mass_matrix*du1 - du
+            end
+	        EEst = norm(du2) / (integrator.opts.abstol + integrator.opts.reltol*norm(k2)) 
+            integrator.EEst = max(EEst,integrator.EEst)
+        end
     end
 
     integrator.u = u
@@ -2409,10 +2460,15 @@ end
     @.. broadcast=false veck8=-vecu
     integrator.stats.nsolve += 1
 
+    du .= k8
     u .+= k8
 
     if integrator.opts.adaptive
-        calculate_residuals!(atmp, k8, uprev, u, integrator.opts.abstol,
+	    if (integrator.alg isa Rodas5Pe)
+            du = 0.2606326497975715*k1 - 0.005158627295444251*k2 + 1.3038988631109731*k3 + 1.235000722062074*k4 +
+               - 0.7931985603795049*k5 - 1.005448461135913*k6 - 0.18044626132120234*k7 + 0.17051519239113755*k8
+	    end
+        calculate_residuals!(atmp, du, uprev, u, integrator.opts.abstol,
             integrator.opts.reltol, integrator.opts.internalnorm, t)
         integrator.EEst = integrator.opts.internalnorm(atmp, t)
     end
@@ -2425,6 +2481,20 @@ end
                                             h35 * k5 + h36 * k6 + h37 * k7 + h38 * k8
         @.. broadcast=false integrator.k[3]=h41 * k1 + h42 * k2 + h43 * k3 + h44 * k4 +
                                             h45 * k5 + h46 * k6 + h47 * k7 + h48 * k8
+        if (integrator.alg isa Rodas5Pr) && integrator.opts.adaptive && (integrator.EEst < 1.0)
+            k2 = 0.5*(uprev + u + 0.5 * (integrator.k[1] + 0.5 * (integrator.k[2] + 0.5 * integrator.k[3])))
+            du1 = ( 0.25*(integrator.k[2] + integrator.k[3]) - uprev + u) / dt
+            f(du, k2, p, t + dt/2)
+            integrator.stats.nf += 1
+            if mass_matrix === I
+                du2 = du1 - du
+            else
+                mul!(_vec(du2), mass_matrix, _vec(du1))
+                du2 = du2 - du
+            end
+	        EEst = norm(du2) / (integrator.opts.abstol + integrator.opts.reltol*norm(k2)) 
+            integrator.EEst = max(EEst,integrator.EEst)
+        end
     end
     cache.linsolve = linres.cache
 end
@@ -2709,10 +2779,17 @@ end
 
     @inbounds @simd ivdep for i in eachindex(u)
         u[i] += k8[i]
+        du[i] = k8[i]
     end
 
     if integrator.opts.adaptive
-        calculate_residuals!(atmp, k8, uprev, u, integrator.opts.abstol,
+	    if (integrator.alg isa Rodas5Pe)
+	    	@inbounds @simd ivdep for i in eachindex(u)
+               	du[i] = 0.2606326497975715*k1[i] - 0.005158627295444251*k2[i] + 1.3038988631109731*k3[i] + 1.235000722062074*k4[i] +
+                   	- 0.7931985603795049*k5[i] - 1.005448461135913*k6[i] - 0.18044626132120234*k7[i] + 0.17051519239113755*k8[i]
+	        end					 
+        end
+        calculate_residuals!(atmp, du, uprev, u, integrator.opts.abstol,
             integrator.opts.reltol, integrator.opts.internalnorm, t)
         integrator.EEst = integrator.opts.internalnorm(atmp, t)
     end
@@ -2726,6 +2803,26 @@ end
                                  h35 * k5[i] + h36 * k6[i] + h37 * k7[i] + h38 * k8[i]
             integrator.k[3][i] = h41 * k1[i] + h42 * k2[i] + h43 * k3[i] + h44 * k4[i] +
                                  h45 * k5[i] + h46 * k6[i] + h47 * k7[i] + h48 * k8[i]
+	        if (integrator.alg isa Rodas5Pr)
+                k2[i] = 0.5*(uprev[i] + u[i] + 0.5 * (integrator.k[1][i] + 0.5 * (integrator.k[2][i] + 0.5 * integrator.k[3][i])))
+                du1[i] = ( 0.25*(integrator.k[2][i] + integrator.k[3][i]) - uprev[i] + u[i]) / dt
+	        end
+        end
+        if integrator.opts.adaptive && (integrator.EEst < 1.0) && (integrator.alg isa Rodas5Pr) 
+            f(du, k2, p, t + dt/2)
+            integrator.stats.nf += 1
+            if mass_matrix === I
+                @inbounds @simd ivdep for i in eachindex(u)
+                    du2[i] = du1[i] - du[i]
+                end
+            else
+                mul!(_vec(du2), mass_matrix, _vec(du1))
+                @inbounds @simd ivdep for i in eachindex(u)
+                    du2[i] = du2[i] - du[i]
+                end
+            end
+	        EEst = norm(du2) / (integrator.opts.abstol + integrator.opts.reltol*norm(k2)) 
+            integrator.EEst = max(EEst,integrator.EEst)
         end
     end
     cache.linsolve = linres.cache
diff --git a/test/algconvergence/ode_rosenbrock_tests.jl b/test/algconvergence/ode_rosenbrock_tests.jl
index 067bde3dfc..257d424af9 100644
--- a/test/algconvergence/ode_rosenbrock_tests.jl
+++ b/test/algconvergence/ode_rosenbrock_tests.jl
@@ -435,6 +435,49 @@ import LinearSolve
     sim = test_convergence(dts, prob, RosenbrockW6S4OS())#test non-inplace
     @test sim.𝒪est[:final]≈4 atol=testTol
 
+    ### Rodas23W, Rodas3P
+
+    println("Rodas23W")
+
+    prob = prob_ode_linear
+
+    dts = (1 / 2) .^ (6:-1:3)
+    sim = test_convergence(dts, prob, Rodas23W(), dense_errors = true)
+    @test sim.𝒪est[:final]≈2 atol=testTol
+    @test sim.𝒪est[:L2]≈2 atol=testTol
+
+    sol = solve(prob, Rodas23W())
+    @test length(sol) < 20
+
+    prob = prob_ode_2Dlinear
+
+    sim = test_convergence(dts, prob, Rodas23W(), dense_errors = true)
+    @test sim.𝒪est[:final]≈2 atol=testTol
+    @test sim.𝒪est[:L2]≈2 atol=testTol
+
+    sol = solve(prob, Rodas23W())
+    @test length(sol) < 20
+
+    println("Rodas3P")
+
+    prob = prob_ode_linear
+
+    sim = test_convergence(dts, prob, Rodas3P(), dense_errors = true)
+    @test sim.𝒪est[:final]≈3 atol=testTol
+    @test sim.𝒪est[:L2]≈3 atol=testTol
+
+    sol = solve(prob, Rodas3P())
+    @test length(sol) < 20
+
+    prob = prob_ode_2Dlinear
+
+    sim = test_convergence(dts, prob, Rodas3P(), dense_errors = true)
+    @test sim.𝒪est[:final]≈3 atol=testTol
+    @test sim.𝒪est[:L2]≈3 atol=testTol
+
+    sol = solve(prob, Rodas3P())
+    @test length(sol) < 20
+
     ### Rodas4 Algorithms
 
     println("RODAS")
@@ -570,21 +613,61 @@ import LinearSolve
 
     prob = prob_ode_linear
 
-    dts = (1 / 2) .^ (6:-1:3)
-    sim = test_convergence(dts, prob, Rodas5(), dense_errors = true)
-    @test sim.𝒪est[:final]≈5 atol=testTol
+    dts = (1 / 2) .^ (5:-1:2)
+    sim = test_convergence(dts, prob, Rodas5P(), dense_errors = true)
+    #@test sim.𝒪est[:final]≈5 atol=testTol #-- observed order > 6
     @test sim.𝒪est[:L2]≈5 atol=testTol
 
-    sol = solve(prob, Rodas5())
+    sol = solve(prob, Rodas5P())
     @test length(sol) < 20
 
     prob = prob_ode_2Dlinear
 
-    sim = test_convergence(dts, prob, Rodas5(), dense_errors = true)
-    @test sim.𝒪est[:final]≈5 atol=testTol
+    sim = test_convergence(dts, prob, Rodas5P(), dense_errors = true)
+    #@test sim.𝒪est[:final]≈5 atol=testTol #-- observed order > 6
     @test sim.𝒪est[:L2]≈5 atol=testTol
 
-    sol = solve(prob, Rodas5())
+    sol = solve(prob, Rodas5P())
+    @test length(sol) < 20
+    
+    println("Rodas5Pe")
+
+    prob = prob_ode_linear
+
+    sim = test_convergence(dts, prob, Rodas5Pe(), dense_errors = true)
+    #@test sim.𝒪est[:final]≈5 atol=testTol #-- observed order > 6
+    @test sim.𝒪est[:L2]≈5 atol=testTol
+
+    sol = solve(prob, Rodas5Pe())
+    @test length(sol) < 20
+
+    prob = prob_ode_2Dlinear
+
+    sim = test_convergence(dts, prob, Rodas5Pe(), dense_errors = true)
+    #@test sim.𝒪est[:final]≈5 atol=testTol #-- observed order > 6
+    @test sim.𝒪est[:L2]≈5 atol=testTol
+
+    sol = solve(prob, Rodas5Pe())
+    @test length(sol) < 20
+
+    println("Rodas5Pr")
+
+    prob = prob_ode_linear
+
+    sim = test_convergence(dts, prob, Rodas5Pr(), dense_errors = true)
+    #@test sim.𝒪est[:final]≈5 atol=testTol #-- observed order > 6
+    @test sim.𝒪est[:L2]≈5 atol=testTol
+
+    sol = solve(prob, Rodas5Pr())
+    @test length(sol) < 20
+
+    prob = prob_ode_2Dlinear
+
+    sim = test_convergence(dts, prob, Rodas5Pr(), dense_errors = true)
+    #@test sim.𝒪est[:final]≈5 atol=testTol #-- observed order > 6
+    @test sim.𝒪est[:L2]≈5 atol=testTol
+
+    sol = solve(prob, Rodas5Pr())
     @test length(sol) < 20
 
     prob = ODEProblem((u, p, t) -> 0.9u, 0.1, (0.0, 1.0))
diff --git a/test/interface/mass_matrix_tests.jl b/test/interface/mass_matrix_tests.jl
index 955ebb604e..ff177d5f96 100644
--- a/test/interface/mass_matrix_tests.jl
+++ b/test/interface/mass_matrix_tests.jl
@@ -131,6 +131,11 @@ dependent_M2 = MatrixOperator(ones(3, 3), update_func = update_func2,
         @test _norm_dsol(Rodas42(), prob, prob2)≈0 atol=1e-9
         @test _norm_dsol(Rodas4P(), prob, prob2)≈0 atol=1e-9
         @test _norm_dsol(Rodas5(), prob, prob2)≈0 atol=1e-7
+        @test _norm_dsol(Rodas23W(), prob, prob2)≈0 atol=1e-11
+        @test _norm_dsol(Rodas3P(), prob, prob2)≈0 atol=1e-11
+        @test _norm_dsol(Rodas5P(), prob, prob2)≈0 atol=1e-7
+        @test _norm_dsol(Rodas5Pe(), prob, prob2)≈0 atol=1e-7
+        @test _norm_dsol(Rodas5Pr(), prob, prob2)≈0 atol=1e-7
     end
 end
 
diff --git a/test/regression/iipvsoop_tests.jl b/test/regression/iipvsoop_tests.jl
index fdf3e40a7f..df0fb77696 100644
--- a/test/regression/iipvsoop_tests.jl
+++ b/test/regression/iipvsoop_tests.jl
@@ -116,6 +116,7 @@ end
 working_rosenbrock_algs = [Rosenbrock23(), ROS3P(), Rodas3(),
     RosShamp4(), Veldd4(), Velds4(), GRK4T(), GRK4A(),
     Ros4LStab(), Rodas4(), Rodas42(), Rodas4P(), Rodas5(),
+    Rodas23W(), Rodas3P(), Rodas5Pe(), Rodas5P(),
     ROS2(), ROS2PR(), ROS2S(), ROS3(), ROS3PR(), Scholz4_7(),
     ROS34PW1a(), ROS34PW1b(), ROS34PW2(), ROS34PW3(), 
     ROS34PRw(), ROS3PRL(), ROS3PRL2()] 
diff --git a/test/regression/ode_dense_tests.jl b/test/regression/ode_dense_tests.jl
index 67ac5a738b..b633716254 100644
--- a/test/regression/ode_dense_tests.jl
+++ b/test/regression/ode_dense_tests.jl
@@ -525,6 +525,12 @@ regression_test(Rosenbrock23(), 3e-3, 6e-3; test_diff1 = true, nth_der = 1, dert
 # Rosenbrock32
 regression_test(Rosenbrock32(), 6e-4, 9e-4; test_diff1 = true, nth_der = 1, dertol = 1e-14)
 
+# Rodas23W
+regression_test(Rodas23W(), 2e-3, 4e-3, test_diff1 = true, nth_der = 1, dertol = 1e-14)
+
+# Rodas3P
+regression_test(Rodas3P(), 2e-4, 4e-4, test_diff1 = true, nth_der = 1, dertol = 1e-14)
+
 # Rodas4
 regression_test(Rodas4(), 8.5e-6, 2e-5, test_diff1 = true, nth_der = 1, dertol = 1e-14)
 
@@ -543,6 +549,12 @@ regression_test(Rodas5(), 2e-6, 3e-6, test_diff1 = true, nth_der = 3, dertol = 5
 # Rodas5P
 regression_test(Rodas5P(), 2e-5, 3e-5, test_diff1 = true, nth_der = 3, dertol = 5e-1)
 
+# Rodas5Pe
+regression_test(Rodas5Pe(), 2e-5, 3e-5, test_diff1 = true, nth_der = 3, dertol = 5e-1)
+
+# Rodas5Pr
+regression_test(Rodas5Pr(), 2e-5, 3e-5, test_diff1 = true, nth_der = 3, dertol = 5e-1)
+
 # ExplicitRK
 regression_test(ExplicitRK(), 7e-5, 2e-4)