From 52f791c8312187dc49bd205ad51e8c60663b968e Mon Sep 17 00:00:00 2001
From: jClugstor <jadonclugston@gmail.com>
Date: Mon, 10 Nov 2025 15:16:40 -0500
Subject: [PATCH 01/18] add verbose argument to alg_cache to be compatible with
 OrdinaryDiffEqCore __init

---
 src/mpdec.jl   |  4 ++--
 src/mprk.jl    | 12 ++++++------
 src/sspmprk.jl |  8 ++++----
 3 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/mpdec.jl b/src/mpdec.jl
index 5fb83dbbd..449dd28de 100644
--- a/src/mpdec.jl
+++ b/src/mpdec.jl
@@ -242,7 +242,7 @@ end
 function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("MPDeC can only be applied to production-destruction systems"))
     end
@@ -641,7 +641,7 @@ get_tmp_cache(integrator, ::MPDeC, cache::OrdinaryDiffEqMutableCache) = (cache.
 function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     nodes, theta = get_constant_parameters(alg, uEltypeNoUnits)
     tab = MPDeCConstantCache(alg.K, alg.M, nodes, theta,
                              alg.small_constant_function(uEltypeNoUnits), Val(alg.M))
diff --git a/src/mprk.jl b/src/mprk.jl
index 1328d86d1..04032c9c8 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -282,7 +282,7 @@ end
 function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("MPE can only be applied to production-destruction systems"))
     end
@@ -347,7 +347,7 @@ get_tmp_cache(integrator, ::MPE, cache::OrdinaryDiffEqMutableCache) = (cache.σ,
 function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     P = p_prototype(u, f)
     σ = zero(u)
     tab = MPEConstantCache(alg.small_constant_function(uEltypeNoUnits))
@@ -540,7 +540,7 @@ end
 function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("MPRK22 can only be applied to production-destruction systems"))
     end
@@ -652,7 +652,7 @@ get_tmp_cache(integrator, ::MPRK22, cache::OrdinaryDiffEqMutableCache) = (cache.
 function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     a21, b1, b2 = get_constant_parameters(alg)
     tab = MPRK22ConstantCache(a21, b1, b2, alg.small_constant_function(uEltypeNoUnits))
     tmp = zero(u)
@@ -1084,7 +1084,7 @@ end
 function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("MPRK43 can only be applied to production-destruction systems"))
     end
@@ -1250,7 +1250,7 @@ end
 function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2 = get_constant_parameters(alg)
     tab = MPRK43ConstantCache(a21, a31, a32, b1, b2, b3, c2, c3,
                               beta1, beta2, q1, q2,
diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index f3a325a53..a8f2f362d 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -106,7 +106,7 @@ end
 function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("SSPMPRK22 can only be applied to production-destruction systems"))
     end
@@ -227,7 +227,7 @@ get_tmp_cache(integrator, ::SSPMPRK22, cache::OrdinaryDiffEqMutableCache) = (cac
 function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     a21, a10, a20, b10, b20, b21, s, τ = get_constant_parameters(alg)
     tab = SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, τ,
                                  alg.small_constant_function(uEltypeNoUnits))
@@ -589,7 +589,7 @@ end
 function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("SSPMPRK43 can only be applied to production-destruction systems"))
     end
@@ -761,7 +761,7 @@ get_tmp_cache(integrator, ::SSPMPRK43, cache::OrdinaryDiffEqMutableCache) = (cac
 function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = get_constant_parameters(alg)
     tab = SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30,
                                  α31, α32,

From b2d4ab6e92125842955eb12b45b21b52f7eb1476 Mon Sep 17 00:00:00 2001
From: jClugstor <jadonclugston@gmail.com>
Date: Thu, 15 Jan 2026 17:55:51 -0500
Subject: [PATCH 02/18] bump OrdinaryDiffEqCore compat for alg_cache

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 982e9a49b..fba65a98f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -30,7 +30,7 @@ JuMP = "1.28"
 LinearAlgebra = "1"
 LinearSolve = "3.7.1"
 MuladdMacro = "0.2.4"
-OrdinaryDiffEqCore = "1.21"
+OrdinaryDiffEqCore = "1.21, 3"
 RecipesBase = "1.3.4"
 Reexport = "1.2.2"
 SciMLBase = "2.78"

From cb08702c16075a87320ace8854a6ed5c8d846b1a Mon Sep 17 00:00:00 2001
From: jClugstor <jadonclugston@gmail.com>
Date: Thu, 15 Jan 2026 18:16:59 -0500
Subject: [PATCH 03/18] use just version 3 for OrdinaryDiffEqCore

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index fba65a98f..16092a56a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -30,7 +30,7 @@ JuMP = "1.28"
 LinearAlgebra = "1"
 LinearSolve = "3.7.1"
 MuladdMacro = "0.2.4"
-OrdinaryDiffEqCore = "1.21, 3"
+OrdinaryDiffEqCore = "3"
 RecipesBase = "1.3.4"
 Reexport = "1.2.2"
 SciMLBase = "2.78"

From 9256967439c2ca6e9e248274be134c9614063714 Mon Sep 17 00:00:00 2001
From: jClugstor <jadonclugston@gmail.com>
Date: Thu, 29 Jan 2026 11:42:36 -0500
Subject: [PATCH 04/18] bump lower compat for SciMLBase to match
 OrdinaryDiffEqCore

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 9ccbd31eb..8ac027870 100644
--- a/Project.toml
+++ b/Project.toml
@@ -32,7 +32,7 @@ MuladdMacro = "0.2.4"
 OrdinaryDiffEqCore = "3"
 RecipesBase = "1.3.4"
 Reexport = "1.2.2"
-SciMLBase = "2.78"
+SciMLBase = "2.116"
 SparseArrays = "1"
 StaticArrays = "1.9.7"
 Statistics = "1"

From 978c1df744724b3d5c98fa472a128ee223aae561 Mon Sep 17 00:00:00 2001
From: jClugstor <jadonclugston@gmail.com>
Date: Thu, 29 Jan 2026 11:42:48 -0500
Subject: [PATCH 05/18] format

---
 ext/JuMPExt.jl       |  2 +-
 src/interpolation.jl |  8 ++--
 src/mpdec.jl         | 31 ++++++++-------
 src/mprk.jl          | 58 +++++++++++++++-------------
 src/sspmprk.jl       | 91 ++++++++++++++++++++++++--------------------
 test/runtests.jl     | 11 +++---
 6 files changed, 111 insertions(+), 90 deletions(-)

diff --git a/ext/JuMPExt.jl b/ext/JuMPExt.jl
index 7c6b051d2..e1c8ab942 100644
--- a/ext/JuMPExt.jl
+++ b/ext/JuMPExt.jl
@@ -74,7 +74,7 @@ function SanduProjection(model, AT, b, eps = nothing; save = true, verbose = fal
     set_string_names_on_creation(model, false)
 
     @variable(model, z[i = 1:s]>=epsv[i])
-    @constraint(model, AT * z.==b)
+    @constraint(model, AT * z .== b)
     # This just initializes the objective. The correct coefficients will be set later
     @expression(model, obj_exp, sum(z .^ 2)+sum(z))
     @objective(model, Min, obj_exp)
diff --git a/src/interpolation.jl b/src/interpolation.jl
index d257009da..40850515e 100644
--- a/src/interpolation.jl
+++ b/src/interpolation.jl
@@ -12,13 +12,13 @@ end
 @muladd @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs::Nothing,
                                              T::Type{Val{0}})
     Θm1 = (1 - Θ)
-    @.. broadcast=false out=Θm1 * u0 + Θ * u1
+    @.. broadcast=false out=Θm1*u0+Θ*u1
     out
 end
 
 @muladd @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs, T::Type{Val{0}})
     Θm1 = (1 - Θ)
-    @views @.. broadcast=false out=Θm1 * u0[idxs] + Θ * u1[idxs]
+    @views @.. broadcast=false out=Θm1*u0[idxs]+Θ*u1[idxs]
     out
 end
 
@@ -31,12 +31,12 @@ end
 end
 
 @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs::Nothing, T::Type{Val{1}})
-    @.. broadcast=false out=(u1 - u0) / dt
+    @.. broadcast=false out=(u1-u0)/dt
     out
 end
 
 @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs, T::Type{Val{1}})
-    @views @.. broadcast=false out=(u1[idxs] - u0[idxs]) / dt
+    @views @.. broadcast=false out=(u1[idxs]-u0[idxs])/dt
     out
 end
 
diff --git a/src/mpdec.jl b/src/mpdec.jl
index 87d6d88b1..5d925586b 100644
--- a/src/mpdec.jl
+++ b/src/mpdec.jl
@@ -242,7 +242,8 @@ end
 function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("MPDeC can only be applied to production-destruction systems"))
     end
@@ -260,14 +261,16 @@ function build_mpdec_matrix_and_rhs_oop(uprev, m, f, C, p, t, dt, nodes, theta,
                                         small_constant)
     if f isa PDSFunction
         # Additional destruction terms
-        Mmat, rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes,
-                                                    theta,
-                                                    small_constant, f.d)
+        Mmat,
+        rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes,
+                                              theta,
+                                              small_constant, f.d)
     else
         # No additional destruction terms
-        Mmat, rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes,
-                                                    theta,
-                                                    small_constant)
+        Mmat,
+        rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes,
+                                              theta,
+                                              small_constant)
     end
 
     return Mmat, rhs
@@ -577,9 +580,10 @@ end
     for _ in 1:K
         C .= C2
         for m in 2:(M + 1)
-            Mmat, rhs = build_mpdec_matrix_and_rhs_oop(uprev, m, f, C, p, t, dt, nodes,
-                                                       theta,
-                                                       small_constant)
+            Mmat,
+            rhs = build_mpdec_matrix_and_rhs_oop(uprev, m, f, C, p, t, dt, nodes,
+                                                 theta,
+                                                 small_constant)
             # solve linear system
             linprob = LinearProblem(Mmat, rhs)
             sol = solve(linprob, alg.linsolve)
@@ -641,7 +645,8 @@ get_tmp_cache(integrator, ::MPDeC, cache::OrdinaryDiffEqMutableCache) = (cache.
 function alg_cache(alg::MPDeC, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{true},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     nodes, theta = get_constant_parameters(alg, uEltypeNoUnits)
     tab = MPDeCConstantCache(alg.K, alg.M, nodes, theta,
                              alg.small_constant_function(uEltypeNoUnits), Val(alg.M))
@@ -732,7 +737,7 @@ end
     σ .= C[:, M + 1] # one order less accurate
 
     # Now σ stores the error estimate
-    @.. broadcast=false σ=u - σ
+    @.. broadcast=false σ=u-σ
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
@@ -776,7 +781,7 @@ end
     σ .= C[:, M + 1] # one order less accurate
 
     # Now σ stores the error estimate
-    @.. broadcast=false σ=u - σ
+    @.. broadcast=false σ=u-σ
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
diff --git a/src/mprk.jl b/src/mprk.jl
index 5eab9b0c9..97e729de4 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -454,7 +454,8 @@ end
 function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("MPE can only be applied to production-destruction systems"))
     end
@@ -503,7 +504,8 @@ get_tmp_cache(integrator, ::MPE, cache::OrdinaryDiffEqMutableCache) = (cache.σ,
 function alg_cache(alg::MPE, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{true},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     P = p_prototype(u, f)
     σ = zero(u)
     tab = MPEConstantCache(alg.small_constant_function(uEltypeNoUnits))
@@ -553,7 +555,7 @@ end
     integrator.stats.nf += 1
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev + small_constant
+    @.. broadcast=false σ=uprev+small_constant
 
     basic_patankar_step!(u, uprev, P, D, σ, dt, linsolve)
     integrator.stats.nsolve += 1
@@ -573,7 +575,7 @@ end
     integrator.stats.nf += 1
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev + small_constant
+    @.. broadcast=false σ=uprev+small_constant
 
     basic_patankar_step_conservative!(u, uprev, P, σ, dt, linsolve)
     integrator.stats.nsolve += 1
@@ -679,7 +681,8 @@ end
 function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("MPRK22 can only be applied to production-destruction systems"))
     end
@@ -761,7 +764,8 @@ get_tmp_cache(integrator, ::MPRK22, cache::OrdinaryDiffEqMutableCache) = (cache.
 function alg_cache(alg::MPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{true},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     a21, b1, b2 = get_constant_parameters(alg)
     tab = MPRK22ConstantCache(a21, b1, b2, alg.small_constant_function(uEltypeNoUnits))
     tmp = zero(u)
@@ -822,16 +826,16 @@ end
     lincomb!(D2, a21, D)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev + small_constant
+    @.. broadcast=false σ=uprev+small_constant
 
     tmp .= uprev
     basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     if isone(a21)
-        @.. broadcast=false σ=u + small_constant
+        @.. broadcast=false σ=u+small_constant
     else
-        @.. broadcast=false σ=σ^(1 - 1 / a21) * u^(1 / a21) + small_constant
+        @.. broadcast=false σ=σ^(1-1/a21)*u^(1/a21)+small_constant
     end
 
     f.p(P2, u, p, t + a21 * dt) # evaluate production terms
@@ -848,7 +852,7 @@ end
     # Now σ stores the error estimate
     # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems.
     # Otherwise, this is not clear.
-    @.. broadcast=false σ=u - σ
+    @.. broadcast=false σ=u-σ
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
@@ -871,16 +875,16 @@ end
     lincomb!(P2, a21, P)
 
     # Avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev + small_constant
+    @.. broadcast=false σ=uprev+small_constant
 
     tmp .= uprev
     basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     if isone(a21)
-        @.. broadcast=false σ=u + small_constant
+        @.. broadcast=false σ=u+small_constant
     else
-        @.. broadcast=false σ=σ^(1 - 1 / a21) * u^(1 / a21) + small_constant
+        @.. broadcast=false σ=σ^(1-1/a21)*u^(1/a21)+small_constant
     end
 
     f.p(P2, u, p, t + a21 * dt) # evaluate production terms
@@ -894,7 +898,7 @@ end
     # Now σ stores the error estimate
     # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems.
     # Otherwise, this is not clear.
-    @.. broadcast=false σ=u - σ
+    @.. broadcast=false σ=u-σ
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
@@ -1121,7 +1125,8 @@ end
 function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("MPRK43 can only be applied to production-destruction systems"))
     end
@@ -1230,7 +1235,8 @@ end
 function alg_cache(alg::Union{MPRK43I, MPRK43II}, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{true},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     a21, a31, a32, b1, b2, b3, c2, c3, beta1, beta2, q1, q2 = get_constant_parameters(alg)
     tab = MPRK43ConstantCache(a21, a31, a32, b1, b2, b3, c2, c3,
                               beta1, beta2, q1, q2,
@@ -1291,7 +1297,7 @@ end
     lincomb!(D3, a21, D)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev + small_constant
+    @.. broadcast=false σ=uprev+small_constant
 
     tmp .= uprev
     basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve)
@@ -1301,7 +1307,7 @@ end
         tmp2 .= u #u2 in out-of-place version
     end
 
-    @.. broadcast=false σ=σ^(1 - q1) * u^q1 + small_constant
+    @.. broadcast=false σ=σ^(1-q1)*u^q1+small_constant
 
     f.p(P2, u, p, t + c2 * dt) # evaluate production terms
     f.d(D2, u, p, t + c2 * dt) # evaluate nonconservative destruction terms
@@ -1315,7 +1321,7 @@ end
     integrator.stats.nsolve += 1
 
     if !(q1 ≈ q2)
-        @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 + small_constant
+        @.. broadcast=false σ=(uprev+small_constant)^(1-q2)*tmp2^q2+small_constant
     end
 
     lincomb!(P3, beta1, P, beta2, P2)
@@ -1326,7 +1332,7 @@ end
     integrator.stats.nsolve += 1
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=σ + small_constant
+    @.. broadcast=false σ=σ+small_constant
 
     f.p(P3, u, p, t + c3 * dt) # evaluate production terms
     f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms
@@ -1340,7 +1346,7 @@ end
     integrator.stats.nsolve += 1
 
     # Now tmp stores the error estimate
-    @.. broadcast=false tmp=u - σ
+    @.. broadcast=false tmp=u-σ
 
     # Now tmp2 stores error residuals
     calculate_residuals!(tmp2, tmp, uprev, u, integrator.opts.abstol,
@@ -1363,7 +1369,7 @@ end
     lincomb!(P3, a21, P)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev + small_constant
+    @.. broadcast=false σ=uprev+small_constant
 
     tmp .= uprev
     basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve)
@@ -1374,7 +1380,7 @@ end
     end
     integrator.stats.nsolve += 1
 
-    @.. broadcast=false σ=σ^(1 - q1) * u^q1 + small_constant
+    @.. broadcast=false σ=σ^(1-q1)*u^q1+small_constant
 
     f.p(P2, u, p, t + c2 * dt) # evaluate production terms
     integrator.stats.nf += 1
@@ -1385,7 +1391,7 @@ end
     integrator.stats.nsolve += 1
 
     if !(q1 ≈ q2)
-        @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 + small_constant
+        @.. broadcast=false σ=(uprev+small_constant)^(1-q2)*tmp2^q2+small_constant
     end
 
     lincomb!(P3, beta1, P, beta2, P2)
@@ -1394,7 +1400,7 @@ end
     integrator.stats.nsolve += 1
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=σ + small_constant
+    @.. broadcast=false σ=σ+small_constant
 
     f.p(P3, u, p, t + c3 * dt) # evaluate production terms
     integrator.stats.nf += 1
@@ -1405,7 +1411,7 @@ end
     integrator.stats.nsolve += 1
 
     # Now tmp stores the error estimate
-    @.. broadcast=false tmp=u - σ
+    @.. broadcast=false tmp=u-σ
 
     # Now tmp2 stores error residuals
     calculate_residuals!(tmp2, tmp, uprev, u, integrator.opts.abstol,
diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index dffa19040..66eecbd21 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -106,7 +106,8 @@ end
 function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("SSPMPRK22 can only be applied to production-destruction systems"))
     end
@@ -199,7 +200,8 @@ get_tmp_cache(integrator, ::SSPMPRK22, cache::OrdinaryDiffEqMutableCache) = (cac
 function alg_cache(alg::SSPMPRK22, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{true},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     a21, a10, a20, b10, b20, b21, s, τ = get_constant_parameters(alg)
     tab = SSPMPRK22ConstantCache(a21, a10, a20, b10, b20, b21, s, τ,
                                  alg.small_constant_function(uEltypeNoUnits))
@@ -257,17 +259,17 @@ end
     lincomb!(D2, b10, D)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev + small_constant
+    @.. broadcast=false σ=uprev+small_constant
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=a10 * uprev
+    @.. broadcast=false tmp=a10*uprev
     basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     if isone(s)
-        @.. broadcast=false σ=u + small_constant
+        @.. broadcast=false σ=u+small_constant
     else
-        @.. broadcast=false σ=σ^(1 - s) * u^s + small_constant
+        @.. broadcast=false σ=σ^(1-s)*u^s+small_constant
     end
 
     f.p(P2, u, p, t + b10 * dt) # evaluate production terms
@@ -278,7 +280,7 @@ end
     lincomb!(D2, b20, D, b21, D2)
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=a20 * uprev + a21 * u
+    @.. broadcast=false tmp=a20*uprev+a21*u
     basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
@@ -289,7 +291,7 @@ end
     # σ2 may become negative, but still can be used for error estimation.
 
     # Now σ stores the error estimate
-    @.. broadcast=false σ=u - (σ - uprev) / τ - uprev
+    @.. broadcast=false σ=u-(σ-uprev)/τ-uprev
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
@@ -312,17 +314,17 @@ end
     lincomb!(P2, b10, P)
 
     # Avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev + small_constant
+    @.. broadcast=false σ=uprev+small_constant
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=a10 * uprev
+    @.. broadcast=false tmp=a10*uprev
     basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     if isone(s)
-        @.. broadcast=false σ=u + small_constant
+        @.. broadcast=false σ=u+small_constant
     else
-        @.. broadcast=false σ=σ^(1 - s) * u^s + small_constant
+        @.. broadcast=false σ=σ^(1-s)*u^s+small_constant
     end
 
     f.p(P2, u, p, t + b10 * dt) # evaluate production terms
@@ -330,7 +332,7 @@ end
 
     lincomb!(P2, b20, P, b21, P2)
 
-    @.. broadcast=false tmp=a20 * uprev + a21 * u
+    @.. broadcast=false tmp=a20*uprev+a21*u
     basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
@@ -341,7 +343,7 @@ end
     # σ2 may become negative, but still can be used for error estimation.
 
     # Now σ stores the error estimate
-    @.. broadcast=false σ=u - (σ - uprev) / τ - uprev
+    @.. broadcast=false σ=u-(σ-uprev)/τ-uprev
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
@@ -491,13 +493,15 @@ end
 function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{false}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+                   ::Val{false},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     if !(f isa PDSFunction || f isa ConservativePDSFunction)
         throw(ArgumentError("SSPMPRK43 can only be applied to production-destruction systems"))
     end
     const_param = get_constant_parameters(alg)
     const_param = convert.(uEltypeNoUnits, const_param)
-    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = const_param
+    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30,
+    β31, β32, c3 = const_param
     small_constant = alg.small_constant_function(uEltypeNoUnits)
     SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31,
                            α32, β10, β20, β21, β30, β31, β32, c3, small_constant)
@@ -509,7 +513,8 @@ end
 @muladd function perform_step!(integrator, cache::SSPMPRK43ConstantCache,
                                repeat_step = false)
     (; alg, t, dt, uprev, f, p) = integrator
-    (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant) = cache
+    (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21,
+     β30, β31, β32, c3, small_constant) = cache
 
     f = integrator.f
 
@@ -611,8 +616,10 @@ get_tmp_cache(integrator, ::SSPMPRK43, cache::OrdinaryDiffEqMutableCache) = (cac
 function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
                    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
-                   ::Val{true}, verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = get_constant_parameters(alg)
+                   ::Val{true},
+                   verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30,
+    β31, β32, c3 = get_constant_parameters(alg)
     tab = SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30,
                                  α31, α32,
                                  β10, β20, β21, β30, β31, β32, c3,
@@ -657,7 +664,8 @@ end
 @muladd function perform_step!(integrator, cache::SSPMPRK43Cache, repeat_step = false)
     (; t, dt, uprev, u, f, p) = integrator
     (; tmp, tmp2, P, P2, P3, D, D2, D3, σ, ρ, linsolve) = cache
-    (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant) = cache.tab
+    (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21,
+     β30, β31, β32, c3, small_constant) = cache.tab
 
     # We use P3 to store the last evaluation of the PDS
     # as well as to store the system matrix of the linear system
@@ -670,17 +678,17 @@ end
     lincomb!(D3, β10, D)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev + small_constant
+    @.. broadcast=false σ=uprev+small_constant
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=α10 * uprev
+    @.. broadcast=false tmp=α10*uprev
     basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     tmp2 .= u
 
-    @.. broadcast=false ρ=n1 * u + n2 * u^2 / σ
-    @.. broadcast=false ρ=ρ + small_constant
+    @.. broadcast=false ρ=n1*u+n2*u^2/σ
+    @.. broadcast=false ρ=ρ+small_constant
 
     f.p(P2, u, p, t + β10 * dt) # evaluate production terms
     f.d(D2, u, p, t + β10 * dt) # evaluate nonconservative destruction terms
@@ -690,11 +698,11 @@ end
     lincomb!(D3, β20, D, β21, D2)
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=α20 * uprev + α21 * tmp2
+    @.. broadcast=false tmp=α20*uprev+α21*tmp2
     basic_patankar_step!(u, tmp, P3, D3, ρ, dt, linsolve)
     integrator.stats.nsolve += 1
 
-    @.. broadcast=false σ=σ^(1 - s) * tmp2^s + small_constant
+    @.. broadcast=false σ=σ^(1-s)*tmp2^s+small_constant
 
     lincomb!(P3, η3, P, η4, P2)
     lincomb!(D3, η3, D, η4, D2)
@@ -702,7 +710,7 @@ end
     # The next stage is the only stage that is not suited
     # for a direct application of basic_patankar_step!
     # We therefore build the system matrix explicitly.
-    @.. broadcast=false tmp=η1 * uprev + η2 * tmp2
+    @.. broadcast=false tmp=η1*uprev+η2*tmp2
     # see (3.25 f) in original paper
     #=
     @inbounds for i in eachindex(tmp)
@@ -721,9 +729,9 @@ end
     σ .= linres
     integrator.stats.nsolve += 1
 
-    @.. broadcast=false σ=σ + z * uprev * u / ρ
+    @.. broadcast=false σ=σ+z*uprev*u/ρ
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=σ + small_constant
+    @.. broadcast=false σ=σ+small_constant
 
     f.p(P3, u, p, t + c3 * dt) # evaluate production terms
     f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms
@@ -733,7 +741,7 @@ end
     lincomb!(D3, β30, D, β31, D2, β32, D3)
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u
+    @.. broadcast=false tmp=α30*uprev+α31*tmp2+α32*u
     basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
@@ -755,7 +763,8 @@ end
                                repeat_step = false)
     (; t, dt, uprev, u, f, p) = integrator
     (; tmp, tmp2, P, P2, P3, σ, ρ, linsolve) = cache
-    (; n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant) = cache.tab
+    (; n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31,
+     β32, c3, small_constant) = cache.tab
 
     # We use P3 to store the last evaluation of the PDS
     # as well as to store the system matrix of the linear system
@@ -765,48 +774,48 @@ end
     lincomb!(P3, β10, P)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev + small_constant
+    @.. broadcast=false σ=uprev+small_constant
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=α10 * uprev
+    @.. broadcast=false tmp=α10*uprev
 
     basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     tmp2 .= u
 
-    @.. broadcast=false ρ=n1 * u + n2 * u^2 / σ
-    @.. broadcast=false ρ=ρ + small_constant
+    @.. broadcast=false ρ=n1*u+n2*u^2/σ
+    @.. broadcast=false ρ=ρ+small_constant
 
     f.p(P2, u, p, t + β10 * dt) # evaluate production terms
     integrator.stats.nf += 1
 
     lincomb!(P3, β20, P, β21, P2)
 
-    @.. broadcast=false tmp=α20 * uprev + α21 * tmp2
+    @.. broadcast=false tmp=α20*uprev+α21*tmp2
 
     basic_patankar_step_conservative!(u, tmp, P3, ρ, dt, linsolve)
     integrator.stats.nsolve += 1
 
-    @.. broadcast=false σ=σ^(1 - s) * tmp2^s + small_constant
+    @.. broadcast=false σ=σ^(1-s)*tmp2^s+small_constant
 
     lincomb!(P3, η3, P, η4, P2)
 
-    @.. broadcast=false tmp=η1 * uprev + η2 * tmp2
+    @.. broadcast=false tmp=η1*uprev+η2*tmp2
 
     basic_patankar_step_conservative!(σ, tmp, P3, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
-    @.. broadcast=false σ=σ + z * uprev * u / ρ
+    @.. broadcast=false σ=σ+z*uprev*u/ρ
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=σ + small_constant
+    @.. broadcast=false σ=σ+small_constant
 
     f.p(P3, u, p, t + c3 * dt) # evaluate production terms
     integrator.stats.nf += 1
 
     lincomb!(P3, β30, P, β31, P2, β32, P3)
 
-    @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u
+    @.. broadcast=false tmp=α30*uprev+α31*tmp2+α32*u
     basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
diff --git a/test/runtests.jl b/test/runtests.jl
index faf5862c8..976a43561 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1581,7 +1581,7 @@ end
 
             @testset "$alg" for alg in algs
                 for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!),
-                                          (dest_1!, dest_2!, dest_3!))
+                        (dest_1!, dest_2!, dest_3!))
                     prod = (u, p, t) -> begin
                         P = similar(u, (length(u), length(u)))
                         prod!(P, u, p, t)
@@ -1711,7 +1711,7 @@ end
             rtol = sqrt(eps(Float32))
             @testset "$alg" for alg in algs
                 for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!),
-                                          (dest_1!, dest_2!, dest_3!))
+                        (dest_1!, dest_2!, dest_3!))
                     prod! = prod_3!
                     dest! = dest_3!
                     prod = (u, p, t) -> begin
@@ -1827,8 +1827,8 @@ end
             #solve and test
             for alg in algs
                 for prob in (prob_default, prob_tridiagonal, prob_dense, prob_sparse,
-                             prob_default2,
-                             prob_tridiagonal2, prob_dense2, prob_sparse2)
+                     prob_default2,
+                     prob_tridiagonal2, prob_dense2, prob_sparse2)
                     sol1 = solve(prob, alg; dt, adaptive = false)
 
                     # test get_tmp_cache and integrator interface - modifying
@@ -2174,7 +2174,8 @@ end
                 dest_noncons!(d, u, p, t)
                 return SVector{2}(d)
             end
-            fct_noncons! = (du, u, p, t) -> begin
+            fct_noncons! = (du, u, p,
+                            t) -> begin
                 du[1] = sin(t)^2 * u[2] - cos(2 * t)^2 * u[1] - sin(2 * t)^2 * u[1]
                 du[2] = -sin(t)^2 * u[2] + cos(2 * t)^2 * u[1] + cos(t)^2 * u[2] -
                         sin(0.5 * t)^2 * u[2]

From 602d3ca6c5e9783f5fc0b6df5d1ebfc36e313b9a Mon Sep 17 00:00:00 2001
From: jClugstor <jadonclugston@gmail.com>
Date: Thu, 29 Jan 2026 11:51:27 -0500
Subject: [PATCH 06/18] reformat with JuliaFormatter v1.0.60

---
 ext/JuMPExt.jl       |  2 +-
 src/interpolation.jl |  8 ++---
 src/mpdec.jl         |  4 +--
 src/mprk.jl          | 40 ++++++++++++-------------
 src/sspmprk.jl       | 70 ++++++++++++++++++++++----------------------
 test/runtests.jl     | 10 +++----
 6 files changed, 67 insertions(+), 67 deletions(-)

diff --git a/ext/JuMPExt.jl b/ext/JuMPExt.jl
index e1c8ab942..7c6b051d2 100644
--- a/ext/JuMPExt.jl
+++ b/ext/JuMPExt.jl
@@ -74,7 +74,7 @@ function SanduProjection(model, AT, b, eps = nothing; save = true, verbose = fal
     set_string_names_on_creation(model, false)
 
     @variable(model, z[i = 1:s]>=epsv[i])
-    @constraint(model, AT * z .== b)
+    @constraint(model, AT * z.==b)
     # This just initializes the objective. The correct coefficients will be set later
     @expression(model, obj_exp, sum(z .^ 2)+sum(z))
     @objective(model, Min, obj_exp)
diff --git a/src/interpolation.jl b/src/interpolation.jl
index 40850515e..d257009da 100644
--- a/src/interpolation.jl
+++ b/src/interpolation.jl
@@ -12,13 +12,13 @@ end
 @muladd @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs::Nothing,
                                              T::Type{Val{0}})
     Θm1 = (1 - Θ)
-    @.. broadcast=false out=Θm1*u0+Θ*u1
+    @.. broadcast=false out=Θm1 * u0 + Θ * u1
     out
 end
 
 @muladd @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs, T::Type{Val{0}})
     Θm1 = (1 - Θ)
-    @views @.. broadcast=false out=Θm1*u0[idxs]+Θ*u1[idxs]
+    @views @.. broadcast=false out=Θm1 * u0[idxs] + Θ * u1[idxs]
     out
 end
 
@@ -31,12 +31,12 @@ end
 end
 
 @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs::Nothing, T::Type{Val{1}})
-    @.. broadcast=false out=(u1-u0)/dt
+    @.. broadcast=false out=(u1 - u0) / dt
     out
 end
 
 @inline function linear_interpolant!(out, Θ, dt, u0, u1, idxs, T::Type{Val{1}})
-    @views @.. broadcast=false out=(u1[idxs]-u0[idxs])/dt
+    @views @.. broadcast=false out=(u1[idxs] - u0[idxs]) / dt
     out
 end
 
diff --git a/src/mpdec.jl b/src/mpdec.jl
index 5d925586b..1ad5dda8a 100644
--- a/src/mpdec.jl
+++ b/src/mpdec.jl
@@ -737,7 +737,7 @@ end
     σ .= C[:, M + 1] # one order less accurate
 
     # Now σ stores the error estimate
-    @.. broadcast=false σ=u-σ
+    @.. broadcast=false σ=u - σ
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
@@ -781,7 +781,7 @@ end
     σ .= C[:, M + 1] # one order less accurate
 
     # Now σ stores the error estimate
-    @.. broadcast=false σ=u-σ
+    @.. broadcast=false σ=u - σ
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
diff --git a/src/mprk.jl b/src/mprk.jl
index 97e729de4..1b7e602b5 100644
--- a/src/mprk.jl
+++ b/src/mprk.jl
@@ -555,7 +555,7 @@ end
     integrator.stats.nf += 1
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev+small_constant
+    @.. broadcast=false σ=uprev + small_constant
 
     basic_patankar_step!(u, uprev, P, D, σ, dt, linsolve)
     integrator.stats.nsolve += 1
@@ -575,7 +575,7 @@ end
     integrator.stats.nf += 1
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev+small_constant
+    @.. broadcast=false σ=uprev + small_constant
 
     basic_patankar_step_conservative!(u, uprev, P, σ, dt, linsolve)
     integrator.stats.nsolve += 1
@@ -826,16 +826,16 @@ end
     lincomb!(D2, a21, D)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev+small_constant
+    @.. broadcast=false σ=uprev + small_constant
 
     tmp .= uprev
     basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     if isone(a21)
-        @.. broadcast=false σ=u+small_constant
+        @.. broadcast=false σ=u + small_constant
     else
-        @.. broadcast=false σ=σ^(1-1/a21)*u^(1/a21)+small_constant
+        @.. broadcast=false σ=σ^(1 - 1 / a21) * u^(1 / a21) + small_constant
     end
 
     f.p(P2, u, p, t + a21 * dt) # evaluate production terms
@@ -852,7 +852,7 @@ end
     # Now σ stores the error estimate
     # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems.
     # Otherwise, this is not clear.
-    @.. broadcast=false σ=u-σ
+    @.. broadcast=false σ=u - σ
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
@@ -875,16 +875,16 @@ end
     lincomb!(P2, a21, P)
 
     # Avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev+small_constant
+    @.. broadcast=false σ=uprev + small_constant
 
     tmp .= uprev
     basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     if isone(a21)
-        @.. broadcast=false σ=u+small_constant
+        @.. broadcast=false σ=u + small_constant
     else
-        @.. broadcast=false σ=σ^(1-1/a21)*u^(1/a21)+small_constant
+        @.. broadcast=false σ=σ^(1 - 1 / a21) * u^(1 / a21) + small_constant
     end
 
     f.p(P2, u, p, t + a21 * dt) # evaluate production terms
@@ -898,7 +898,7 @@ end
     # Now σ stores the error estimate
     # If a21 = 1, then σ is the MPE approximation, i.e. suited for stiff problems.
     # Otherwise, this is not clear.
-    @.. broadcast=false σ=u-σ
+    @.. broadcast=false σ=u - σ
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
@@ -1297,7 +1297,7 @@ end
     lincomb!(D3, a21, D)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev+small_constant
+    @.. broadcast=false σ=uprev + small_constant
 
     tmp .= uprev
     basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve)
@@ -1307,7 +1307,7 @@ end
         tmp2 .= u #u2 in out-of-place version
     end
 
-    @.. broadcast=false σ=σ^(1-q1)*u^q1+small_constant
+    @.. broadcast=false σ=σ^(1 - q1) * u^q1 + small_constant
 
     f.p(P2, u, p, t + c2 * dt) # evaluate production terms
     f.d(D2, u, p, t + c2 * dt) # evaluate nonconservative destruction terms
@@ -1321,7 +1321,7 @@ end
     integrator.stats.nsolve += 1
 
     if !(q1 ≈ q2)
-        @.. broadcast=false σ=(uprev+small_constant)^(1-q2)*tmp2^q2+small_constant
+        @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 + small_constant
     end
 
     lincomb!(P3, beta1, P, beta2, P2)
@@ -1332,7 +1332,7 @@ end
     integrator.stats.nsolve += 1
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=σ+small_constant
+    @.. broadcast=false σ=σ + small_constant
 
     f.p(P3, u, p, t + c3 * dt) # evaluate production terms
     f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms
@@ -1346,7 +1346,7 @@ end
     integrator.stats.nsolve += 1
 
     # Now tmp stores the error estimate
-    @.. broadcast=false tmp=u-σ
+    @.. broadcast=false tmp=u - σ
 
     # Now tmp2 stores error residuals
     calculate_residuals!(tmp2, tmp, uprev, u, integrator.opts.abstol,
@@ -1369,7 +1369,7 @@ end
     lincomb!(P3, a21, P)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev+small_constant
+    @.. broadcast=false σ=uprev + small_constant
 
     tmp .= uprev
     basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve)
@@ -1380,7 +1380,7 @@ end
     end
     integrator.stats.nsolve += 1
 
-    @.. broadcast=false σ=σ^(1-q1)*u^q1+small_constant
+    @.. broadcast=false σ=σ^(1 - q1) * u^q1 + small_constant
 
     f.p(P2, u, p, t + c2 * dt) # evaluate production terms
     integrator.stats.nf += 1
@@ -1391,7 +1391,7 @@ end
     integrator.stats.nsolve += 1
 
     if !(q1 ≈ q2)
-        @.. broadcast=false σ=(uprev+small_constant)^(1-q2)*tmp2^q2+small_constant
+        @.. broadcast=false σ=(uprev + small_constant)^(1 - q2) * tmp2^q2 + small_constant
     end
 
     lincomb!(P3, beta1, P, beta2, P2)
@@ -1400,7 +1400,7 @@ end
     integrator.stats.nsolve += 1
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=σ+small_constant
+    @.. broadcast=false σ=σ + small_constant
 
     f.p(P3, u, p, t + c3 * dt) # evaluate production terms
     integrator.stats.nf += 1
@@ -1411,7 +1411,7 @@ end
     integrator.stats.nsolve += 1
 
     # Now tmp stores the error estimate
-    @.. broadcast=false tmp=u-σ
+    @.. broadcast=false tmp=u - σ
 
     # Now tmp2 stores error residuals
     calculate_residuals!(tmp2, tmp, uprev, u, integrator.opts.abstol,
diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index 66eecbd21..46e936755 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -259,17 +259,17 @@ end
     lincomb!(D2, b10, D)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev+small_constant
+    @.. broadcast=false σ=uprev + small_constant
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=a10*uprev
+    @.. broadcast=false tmp=a10 * uprev
     basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     if isone(s)
-        @.. broadcast=false σ=u+small_constant
+        @.. broadcast=false σ=u + small_constant
     else
-        @.. broadcast=false σ=σ^(1-s)*u^s+small_constant
+        @.. broadcast=false σ=σ^(1 - s) * u^s + small_constant
     end
 
     f.p(P2, u, p, t + b10 * dt) # evaluate production terms
@@ -280,7 +280,7 @@ end
     lincomb!(D2, b20, D, b21, D2)
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=a20*uprev+a21*u
+    @.. broadcast=false tmp=a20 * uprev + a21 * u
     basic_patankar_step!(u, tmp, P2, D2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
@@ -291,7 +291,7 @@ end
     # σ2 may become negative, but still can be used for error estimation.
 
     # Now σ stores the error estimate
-    @.. broadcast=false σ=u-(σ-uprev)/τ-uprev
+    @.. broadcast=false σ=u - (σ - uprev) / τ - uprev
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
@@ -314,17 +314,17 @@ end
     lincomb!(P2, b10, P)
 
     # Avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev+small_constant
+    @.. broadcast=false σ=uprev + small_constant
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=a10*uprev
+    @.. broadcast=false tmp=a10 * uprev
     basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     if isone(s)
-        @.. broadcast=false σ=u+small_constant
+        @.. broadcast=false σ=u + small_constant
     else
-        @.. broadcast=false σ=σ^(1-s)*u^s+small_constant
+        @.. broadcast=false σ=σ^(1 - s) * u^s + small_constant
     end
 
     f.p(P2, u, p, t + b10 * dt) # evaluate production terms
@@ -332,7 +332,7 @@ end
 
     lincomb!(P2, b20, P, b21, P2)
 
-    @.. broadcast=false tmp=a20*uprev+a21*u
+    @.. broadcast=false tmp=a20 * uprev + a21 * u
     basic_patankar_step_conservative!(u, tmp, P2, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
@@ -343,7 +343,7 @@ end
     # σ2 may become negative, but still can be used for error estimation.
 
     # Now σ stores the error estimate
-    @.. broadcast=false σ=u-(σ-uprev)/τ-uprev
+    @.. broadcast=false σ=u - (σ - uprev) / τ - uprev
 
     # Now tmp stores error residuals
     calculate_residuals!(tmp, σ, uprev, u, integrator.opts.abstol,
@@ -514,7 +514,7 @@ end
                                repeat_step = false)
     (; alg, t, dt, uprev, f, p) = integrator
     (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21,
-     β30, β31, β32, c3, small_constant) = cache
+    β30, β31, β32, c3, small_constant) = cache
 
     f = integrator.f
 
@@ -665,7 +665,7 @@ end
     (; t, dt, uprev, u, f, p) = integrator
     (; tmp, tmp2, P, P2, P3, D, D2, D3, σ, ρ, linsolve) = cache
     (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21,
-     β30, β31, β32, c3, small_constant) = cache.tab
+    β30, β31, β32, c3, small_constant) = cache.tab
 
     # We use P3 to store the last evaluation of the PDS
     # as well as to store the system matrix of the linear system
@@ -678,17 +678,17 @@ end
     lincomb!(D3, β10, D)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev+small_constant
+    @.. broadcast=false σ=uprev + small_constant
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=α10*uprev
+    @.. broadcast=false tmp=α10 * uprev
     basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     tmp2 .= u
 
-    @.. broadcast=false ρ=n1*u+n2*u^2/σ
-    @.. broadcast=false ρ=ρ+small_constant
+    @.. broadcast=false ρ=n1 * u + n2 * u^2 / σ
+    @.. broadcast=false ρ=ρ + small_constant
 
     f.p(P2, u, p, t + β10 * dt) # evaluate production terms
     f.d(D2, u, p, t + β10 * dt) # evaluate nonconservative destruction terms
@@ -698,11 +698,11 @@ end
     lincomb!(D3, β20, D, β21, D2)
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=α20*uprev+α21*tmp2
+    @.. broadcast=false tmp=α20 * uprev + α21 * tmp2
     basic_patankar_step!(u, tmp, P3, D3, ρ, dt, linsolve)
     integrator.stats.nsolve += 1
 
-    @.. broadcast=false σ=σ^(1-s)*tmp2^s+small_constant
+    @.. broadcast=false σ=σ^(1 - s) * tmp2^s + small_constant
 
     lincomb!(P3, η3, P, η4, P2)
     lincomb!(D3, η3, D, η4, D2)
@@ -710,7 +710,7 @@ end
     # The next stage is the only stage that is not suited
     # for a direct application of basic_patankar_step!
     # We therefore build the system matrix explicitly.
-    @.. broadcast=false tmp=η1*uprev+η2*tmp2
+    @.. broadcast=false tmp=η1 * uprev + η2 * tmp2
     # see (3.25 f) in original paper
     #=
     @inbounds for i in eachindex(tmp)
@@ -729,9 +729,9 @@ end
     σ .= linres
     integrator.stats.nsolve += 1
 
-    @.. broadcast=false σ=σ+z*uprev*u/ρ
+    @.. broadcast=false σ=σ + z * uprev * u / ρ
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=σ+small_constant
+    @.. broadcast=false σ=σ + small_constant
 
     f.p(P3, u, p, t + c3 * dt) # evaluate production terms
     f.d(D3, u, p, t + c3 * dt) # evaluate nonconservative destruction terms
@@ -741,7 +741,7 @@ end
     lincomb!(D3, β30, D, β31, D2, β32, D3)
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=α30*uprev+α31*tmp2+α32*u
+    @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u
     basic_patankar_step!(u, tmp, P3, D3, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
@@ -764,7 +764,7 @@ end
     (; t, dt, uprev, u, f, p) = integrator
     (; tmp, tmp2, P, P2, P3, σ, ρ, linsolve) = cache
     (; n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31,
-     β32, c3, small_constant) = cache.tab
+    β32, c3, small_constant) = cache.tab
 
     # We use P3 to store the last evaluation of the PDS
     # as well as to store the system matrix of the linear system
@@ -774,48 +774,48 @@ end
     lincomb!(P3, β10, P)
 
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=uprev+small_constant
+    @.. broadcast=false σ=uprev + small_constant
 
     # tmp holds the right hand side of the linear system
-    @.. broadcast=false tmp=α10*uprev
+    @.. broadcast=false tmp=α10 * uprev
 
     basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
     tmp2 .= u
 
-    @.. broadcast=false ρ=n1*u+n2*u^2/σ
-    @.. broadcast=false ρ=ρ+small_constant
+    @.. broadcast=false ρ=n1 * u + n2 * u^2 / σ
+    @.. broadcast=false ρ=ρ + small_constant
 
     f.p(P2, u, p, t + β10 * dt) # evaluate production terms
     integrator.stats.nf += 1
 
     lincomb!(P3, β20, P, β21, P2)
 
-    @.. broadcast=false tmp=α20*uprev+α21*tmp2
+    @.. broadcast=false tmp=α20 * uprev + α21 * tmp2
 
     basic_patankar_step_conservative!(u, tmp, P3, ρ, dt, linsolve)
     integrator.stats.nsolve += 1
 
-    @.. broadcast=false σ=σ^(1-s)*tmp2^s+small_constant
+    @.. broadcast=false σ=σ^(1 - s) * tmp2^s + small_constant
 
     lincomb!(P3, η3, P, η4, P2)
 
-    @.. broadcast=false tmp=η1*uprev+η2*tmp2
+    @.. broadcast=false tmp=η1 * uprev + η2 * tmp2
 
     basic_patankar_step_conservative!(σ, tmp, P3, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
-    @.. broadcast=false σ=σ+z*uprev*u/ρ
+    @.. broadcast=false σ=σ + z * uprev * u / ρ
     # avoid division by zero due to zero Patankar weights
-    @.. broadcast=false σ=σ+small_constant
+    @.. broadcast=false σ=σ + small_constant
 
     f.p(P3, u, p, t + c3 * dt) # evaluate production terms
     integrator.stats.nf += 1
 
     lincomb!(P3, β30, P, β31, P2, β32, P3)
 
-    @.. broadcast=false tmp=α30*uprev+α31*tmp2+α32*u
+    @.. broadcast=false tmp=α30 * uprev + α31 * tmp2 + α32 * u
     basic_patankar_step_conservative!(u, tmp, P3, σ, dt, linsolve)
     integrator.stats.nsolve += 1
 
diff --git a/test/runtests.jl b/test/runtests.jl
index 976a43561..6567617bc 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1581,7 +1581,7 @@ end
 
             @testset "$alg" for alg in algs
                 for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!),
-                        (dest_1!, dest_2!, dest_3!))
+                                          (dest_1!, dest_2!, dest_3!))
                     prod = (u, p, t) -> begin
                         P = similar(u, (length(u), length(u)))
                         prod!(P, u, p, t)
@@ -1711,7 +1711,7 @@ end
             rtol = sqrt(eps(Float32))
             @testset "$alg" for alg in algs
                 for (prod!, dest!) in zip((prod_1!, prod_2!, prod_3!),
-                        (dest_1!, dest_2!, dest_3!))
+                                          (dest_1!, dest_2!, dest_3!))
                     prod! = prod_3!
                     dest! = dest_3!
                     prod = (u, p, t) -> begin
@@ -1827,8 +1827,8 @@ end
             #solve and test
             for alg in algs
                 for prob in (prob_default, prob_tridiagonal, prob_dense, prob_sparse,
-                     prob_default2,
-                     prob_tridiagonal2, prob_dense2, prob_sparse2)
+                             prob_default2,
+                             prob_tridiagonal2, prob_dense2, prob_sparse2)
                     sol1 = solve(prob, alg; dt, adaptive = false)
 
                     # test get_tmp_cache and integrator interface - modifying
@@ -2175,7 +2175,7 @@ end
                 return SVector{2}(d)
             end
             fct_noncons! = (du, u, p,
-                            t) -> begin
+            t) -> begin
                 du[1] = sin(t)^2 * u[2] - cos(2 * t)^2 * u[1] - sin(2 * t)^2 * u[1]
                 du[2] = -sin(t)^2 * u[2] + cos(2 * t)^2 * u[1] + cos(t)^2 * u[2] -
                         sin(0.5 * t)^2 * u[2]

From a6947e1563eea8c8c734f04f702a164d0ca37780 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 13:57:28 +0100
Subject: [PATCH 07/18] format

---
 src/mpdec.jl     | 21 +++++++++------------
 src/sspmprk.jl   | 15 +++++----------
 test/runtests.jl |  3 +--
 3 files changed, 15 insertions(+), 24 deletions(-)

diff --git a/src/mpdec.jl b/src/mpdec.jl
index 1ad5dda8a..bda3fc247 100644
--- a/src/mpdec.jl
+++ b/src/mpdec.jl
@@ -261,16 +261,14 @@ function build_mpdec_matrix_and_rhs_oop(uprev, m, f, C, p, t, dt, nodes, theta,
                                         small_constant)
     if f isa PDSFunction
         # Additional destruction terms
-        Mmat,
-        rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes,
-                                              theta,
-                                              small_constant, f.d)
+        Mmat, rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes,
+                                                    theta,
+                                                    small_constant, f.d)
     else
         # No additional destruction terms
-        Mmat,
-        rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes,
-                                              theta,
-                                              small_constant)
+        Mmat, rhs = _build_mpdec_matrix_and_rhs_oop(uprev, m, f.p, C, p, t, dt, nodes,
+                                                    theta,
+                                                    small_constant)
     end
 
     return Mmat, rhs
@@ -580,10 +578,9 @@ end
     for _ in 1:K
         C .= C2
         for m in 2:(M + 1)
-            Mmat,
-            rhs = build_mpdec_matrix_and_rhs_oop(uprev, m, f, C, p, t, dt, nodes,
-                                                 theta,
-                                                 small_constant)
+            Mmat, rhs = build_mpdec_matrix_and_rhs_oop(uprev, m, f, C, p, t, dt, nodes,
+                                                       theta,
+                                                       small_constant)
             # solve linear system
             linprob = LinearProblem(Mmat, rhs)
             sol = solve(linprob, alg.linsolve)
diff --git a/src/sspmprk.jl b/src/sspmprk.jl
index 46e936755..3cff589f9 100644
--- a/src/sspmprk.jl
+++ b/src/sspmprk.jl
@@ -500,8 +500,7 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
     end
     const_param = get_constant_parameters(alg)
     const_param = convert.(uEltypeNoUnits, const_param)
-    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30,
-    β31, β32, c3 = const_param
+    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = const_param
     small_constant = alg.small_constant_function(uEltypeNoUnits)
     SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31,
                            α32, β10, β20, β21, β30, β31, β32, c3, small_constant)
@@ -513,8 +512,7 @@ end
 @muladd function perform_step!(integrator, cache::SSPMPRK43ConstantCache,
                                repeat_step = false)
     (; alg, t, dt, uprev, f, p) = integrator
-    (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21,
-    β30, β31, β32, c3, small_constant) = cache
+    (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant) = cache
 
     f = integrator.f
 
@@ -618,8 +616,7 @@ function alg_cache(alg::SSPMPRK43, u, rate_prototype, ::Type{uEltypeNoUnits},
                    uprev, uprev2, f, t, dt, reltol, p, calck,
                    ::Val{true},
                    verbose) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
-    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30,
-    β31, β32, c3 = get_constant_parameters(alg)
+    n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3 = get_constant_parameters(alg)
     tab = SSPMPRK43ConstantCache(n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30,
                                  α31, α32,
                                  β10, β20, β21, β30, β31, β32, c3,
@@ -664,8 +661,7 @@ end
 @muladd function perform_step!(integrator, cache::SSPMPRK43Cache, repeat_step = false)
     (; t, dt, uprev, u, f, p) = integrator
     (; tmp, tmp2, P, P2, P3, D, D2, D3, σ, ρ, linsolve) = cache
-    (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21,
-    β30, β31, β32, c3, small_constant) = cache.tab
+    (; n1, n2, z, η1, η2, η3, η4, η5, η6, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant) = cache.tab
 
     # We use P3 to store the last evaluation of the PDS
     # as well as to store the system matrix of the linear system
@@ -763,8 +759,7 @@ end
                                repeat_step = false)
     (; t, dt, uprev, u, f, p) = integrator
     (; tmp, tmp2, P, P2, P3, σ, ρ, linsolve) = cache
-    (; n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31,
-    β32, c3, small_constant) = cache.tab
+    (; n1, n2, z, η1, η2, η3, η4, s, α10, α20, α21, α30, α31, α32, β10, β20, β21, β30, β31, β32, c3, small_constant) = cache.tab
 
     # We use P3 to store the last evaluation of the PDS
     # as well as to store the system matrix of the linear system
diff --git a/test/runtests.jl b/test/runtests.jl
index 6567617bc..faf5862c8 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -2174,8 +2174,7 @@ end
                 dest_noncons!(d, u, p, t)
                 return SVector{2}(d)
             end
-            fct_noncons! = (du, u, p,
-            t) -> begin
+            fct_noncons! = (du, u, p, t) -> begin
                 du[1] = sin(t)^2 * u[2] - cos(2 * t)^2 * u[1] - sin(2 * t)^2 * u[1]
                 du[2] = -sin(t)^2 * u[2] + cos(2 * t)^2 * u[1] + cos(t)^2 * u[2] -
                         sin(0.5 * t)^2 * u[2]

From adf7c6d9624f635b358a746552d2492830f50136 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 14:06:11 +0100
Subject: [PATCH 08/18] bump LinearSolve to v3.39

---
 Project.toml      | 2 +-
 test/Project.toml | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 8ac027870..c4d301f76 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,7 +27,7 @@ JuMPExt = "JuMP"
 FastBroadcast = "0.3.5"
 JuMP = "1.28"
 LinearAlgebra = "1"
-LinearSolve = "3.7.1"
+LinearSolve = "3.39"
 MuladdMacro = "0.2.4"
 OrdinaryDiffEqCore = "3"
 RecipesBase = "1.3.4"
diff --git a/test/Project.toml b/test/Project.toml
index a4c714f20..445d75b28 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -29,7 +29,7 @@ DoubleFloats = "1.4.3"
 ExplicitImports = "1.0.1"
 JuMP = "1.28"
 LinearAlgebra = "1"
-LinearSolve = "3.7.1"
+LinearSolve = "3.39"
 OrdinaryDiffEqLowOrderRK = "1.2"
 OrdinaryDiffEqRosenbrock = "1.9"
 OrdinaryDiffEqSDIRK = "1.3"

From 815f0900056828ae371ef8ed8d3910c6d0f8ad40 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 14:10:59 +0100
Subject: [PATCH 09/18] bump SymbolicIndexingInterface.jl to v0.3.39

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index c4d301f76..143cb1fa4 100644
--- a/Project.toml
+++ b/Project.toml
@@ -36,5 +36,5 @@ SciMLBase = "2.116"
 SparseArrays = "1"
 StaticArrays = "1.9.7"
 Statistics = "1"
-SymbolicIndexingInterface = "0.3.36"
+SymbolicIndexingInterface = "0.3.39"
 julia = "1.10"

From 6f8f11038295603e8e7dd0cf651b8fca8b5b5730 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 14:13:51 +0100
Subject: [PATCH 10/18] bump SymbolicIndexingInterface.jl to v0.3.42

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 143cb1fa4..2b8587c9f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -36,5 +36,5 @@ SciMLBase = "2.116"
 SparseArrays = "1"
 StaticArrays = "1.9.7"
 Statistics = "1"
-SymbolicIndexingInterface = "0.3.39"
+SymbolicIndexingInterface = "0.3.42"
 julia = "1.10"

From c864064d54f3e7f09a6842813e615dcdfc43e281 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 14:18:29 +0100
Subject: [PATCH 11/18] bump sublibraries

---
 test/Project.toml | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/test/Project.toml b/test/Project.toml
index 445d75b28..02cc31cb6 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -30,11 +30,11 @@ ExplicitImports = "1.0.1"
 JuMP = "1.28"
 LinearAlgebra = "1"
 LinearSolve = "3.39"
-OrdinaryDiffEqLowOrderRK = "1.2"
-OrdinaryDiffEqRosenbrock = "1.9"
-OrdinaryDiffEqSDIRK = "1.3"
-OrdinaryDiffEqTsit5 = "1.1"
-OrdinaryDiffEqVerner = "1.1"
+OrdinaryDiffEqLowOrderRK = "1.10"
+OrdinaryDiffEqRosenbrock = "1.22"
+OrdinaryDiffEqSDIRK = "1.11"
+OrdinaryDiffEqTsit5 = "1.9"
+OrdinaryDiffEqVerner = "1.10"
 Plots = "1.25.11"
 RecipesBase = "1.3.4"
 RecursiveFactorization = "0.2.23"

From f82c954128b1aa7ca3c07134c9ca86bcd2996b51 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 14:21:32 +0100
Subject: [PATCH 12/18] bump LinearSolve.jl to v3.46

---
 Project.toml      | 2 +-
 test/Project.toml | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 2b8587c9f..2eb26d8f3 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,7 +27,7 @@ JuMPExt = "JuMP"
 FastBroadcast = "0.3.5"
 JuMP = "1.28"
 LinearAlgebra = "1"
-LinearSolve = "3.39"
+LinearSolve = "3.46"
 MuladdMacro = "0.2.4"
 OrdinaryDiffEqCore = "3"
 RecipesBase = "1.3.4"
diff --git a/test/Project.toml b/test/Project.toml
index 02cc31cb6..70e7eae82 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -29,7 +29,7 @@ DoubleFloats = "1.4.3"
 ExplicitImports = "1.0.1"
 JuMP = "1.28"
 LinearAlgebra = "1"
-LinearSolve = "3.39"
+LinearSolve = "3.46"
 OrdinaryDiffEqLowOrderRK = "1.10"
 OrdinaryDiffEqRosenbrock = "1.22"
 OrdinaryDiffEqSDIRK = "1.11"

From 4743546b9136f04c9d83375732f500127adf442c Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 14:23:12 +0100
Subject: [PATCH 13/18] bump ADTypes.jl to v1.16

---
 test/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/Project.toml b/test/Project.toml
index 70e7eae82..1cbcbf3e2 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -22,7 +22,7 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
-ADTypes = "1.14"
+ADTypes = "1.16"
 Aqua = "0.8"
 Clarabel = "0.11"
 DoubleFloats = "1.4.3"

From e9b22c3dd36d8e0ba103793b1808825c4afee09c Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 14:28:10 +0100
Subject: [PATCH 14/18] bump LinearSolve.jl to v3.54

---
 Project.toml      | 2 +-
 test/Project.toml | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/Project.toml b/Project.toml
index 2eb26d8f3..d58eaf696 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,7 +27,7 @@ JuMPExt = "JuMP"
 FastBroadcast = "0.3.5"
 JuMP = "1.28"
 LinearAlgebra = "1"
-LinearSolve = "3.46"
+LinearSolve = "3.54"
 MuladdMacro = "0.2.4"
 OrdinaryDiffEqCore = "3"
 RecipesBase = "1.3.4"
diff --git a/test/Project.toml b/test/Project.toml
index 1cbcbf3e2..6f26056d2 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -29,7 +29,7 @@ DoubleFloats = "1.4.3"
 ExplicitImports = "1.0.1"
 JuMP = "1.28"
 LinearAlgebra = "1"
-LinearSolve = "3.46"
+LinearSolve = "3.54"
 OrdinaryDiffEqLowOrderRK = "1.10"
 OrdinaryDiffEqRosenbrock = "1.22"
 OrdinaryDiffEqSDIRK = "1.11"

From 7d3ea1eb80e2b1c194817523ece776c80d7c979d Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 14:29:50 +0100
Subject: [PATCH 15/18] bump compat of SciMLBase.jl to v2.128

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index d58eaf696..42365a0d1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -32,7 +32,7 @@ MuladdMacro = "0.2.4"
 OrdinaryDiffEqCore = "3"
 RecipesBase = "1.3.4"
 Reexport = "1.2.2"
-SciMLBase = "2.116"
+SciMLBase = "2.128"
 SparseArrays = "1"
 StaticArrays = "1.9.7"
 Statistics = "1"

From 9e1f789e8d36a76160ba56c8f64e47adc73af060 Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 14:31:18 +0100
Subject: [PATCH 16/18] bump RecursiveFactorization.jl to v0.2.26

---
 test/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/Project.toml b/test/Project.toml
index 6f26056d2..2e1cb27e6 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -37,7 +37,7 @@ OrdinaryDiffEqTsit5 = "1.9"
 OrdinaryDiffEqVerner = "1.10"
 Plots = "1.25.11"
 RecipesBase = "1.3.4"
-RecursiveFactorization = "0.2.23"
+RecursiveFactorization = "0.2.26"
 SparseArrays = "1"
 StaticArrays = "1.9.7"
 Statistics = "1"

From ead044ad30480642ee7e289aa27bcd9756e45f6c Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 14:46:36 +0100
Subject: [PATCH 17/18] bump Plots.jl to v1.33

---
 test/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/Project.toml b/test/Project.toml
index 2e1cb27e6..00d623421 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -35,7 +35,7 @@ OrdinaryDiffEqRosenbrock = "1.22"
 OrdinaryDiffEqSDIRK = "1.11"
 OrdinaryDiffEqTsit5 = "1.9"
 OrdinaryDiffEqVerner = "1.10"
-Plots = "1.25.11"
+Plots = "1.33"
 RecipesBase = "1.3.4"
 RecursiveFactorization = "0.2.26"
 SparseArrays = "1"

From f258f86c75df61722d047af8edeb020d9cd79bcc Mon Sep 17 00:00:00 2001
From: Joshua Lampert <joshua.lampert@uni-hamburg.de>
Date: Sun, 1 Feb 2026 14:49:20 +0100
Subject: [PATCH 18/18] bump SymbolicIndexingInterface.jl to v0.3.43

---
 Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 42365a0d1..3d07b89ee 100644
--- a/Project.toml
+++ b/Project.toml
@@ -36,5 +36,5 @@ SciMLBase = "2.128"
 SparseArrays = "1"
 StaticArrays = "1.9.7"
 Statistics = "1"
-SymbolicIndexingInterface = "0.3.42"
+SymbolicIndexingInterface = "0.3.43"
 julia = "1.10"