From 87bd4b06d0e008aac27e0650c9457630ca895e25 Mon Sep 17 00:00:00 2001
From: Orjan Ameye <orjan.ameye@hotmail.com>
Date: Sat, 9 May 2026 08:19:09 +0200
Subject: [PATCH 1/8] Add StrangMarchuk second-order symmetric operator
 splitting
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

Implement the classical ABA scheme for two operators A and B:
A(dt/2) → B(dt) → A(dt/2), achieving second-order accuracy through
symmetry.

Closes #10
---
 src/OrdinaryDiffEqOperatorSplitting.jl |  2 +-
 src/solver.jl                          | 98 ++++++++++++++++++++++++++
 src/utils.jl                           |  8 +++
 test/operator_splitting_api.jl         | 71 ++++++++++++++++++-
 4 files changed, 177 insertions(+), 2 deletions(-)

diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl
index d752233..8c69eb5 100644
--- a/src/OrdinaryDiffEqOperatorSplitting.jl
+++ b/src/OrdinaryDiffEqOperatorSplitting.jl
@@ -35,7 +35,7 @@ include("integrator.jl")
 include("solver.jl")
 include("utils.jl")
 
-export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov
+export GenericSplitFunction, OperatorSplittingProblem, LieTrotterGodunov, StrangMarchuk
 
 include("precompilation.jl")
 
diff --git a/src/solver.jl b/src/solver.jl
index 3163474..1431a2a 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -56,3 +56,101 @@ end
         backward_sync_subintegrator!(parent, child, idxs, sync)
     end
 end
+
+# ---------------------------------------------------------------------------
+# Strang-Marchuk operator splitting
+# ---------------------------------------------------------------------------
+"""
+    StrangMarchuk <: AbstractOperatorSplittingAlgorithm
+
+Second-order symmetric operator splitting algorithm attributed to
+[Str:1968:ccd,Mar:1971:tsm](@cite).
+
+For two operators ``A`` and ``B`` the scheme performs
+``A(\\Delta t/2) \\to B(\\Delta t) \\to A(\\Delta t/2)``,
+achieving second-order accuracy through symmetry.
+"""
+struct StrangMarchuk{AlgTupleType} <: AbstractOperatorSplittingAlgorithm
+    inner_algs::AlgTupleType
+end
+
+function Base.show(io::IO, alg::StrangMarchuk)
+    print(io, "SM (")
+    for inner_alg in alg.inner_algs[1:(end - 1)]
+        Base.show(io, inner_alg)
+        print(io, " -> ")
+    end
+    length(alg.inner_algs) > 0 && Base.show(io, alg.inner_algs[end])
+    return print(io, ")")
+end
+
+struct StrangMarchukCache{uType, uprevType} <: AbstractOperatorSplittingCache
+    u::uType
+    uprev::uprevType
+end
+
+function init_cache(
+        f::GenericSplitFunction, alg::StrangMarchuk;
+        uprev::AbstractArray, u::AbstractVector,
+    )
+    return StrangMarchukCache(u, uprev)
+end
+
+function _perform_step!(
+        parent,
+        children::Tuple,
+        cache::StrangMarchukCache,
+        dt
+    )
+    half_dt = dt / 2
+
+    # A(dt/2)
+    let child = children[1]
+        idxs = parent.child_solution_indices[1]
+        sync = parent.child_synchronizers[1]
+        @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync)
+        @timeit_debug "time solve" advance_solution_by!(parent, child, half_dt)
+        if _child_failed(child)
+            parent.force_stepfail = true
+            return
+        end
+        backward_sync_subintegrator!(parent, child, idxs, sync)
+    end
+
+    # B(dt)
+    let child = children[2]
+        idxs = parent.child_solution_indices[2]
+        sync = parent.child_synchronizers[2]
+        @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync)
+        @timeit_debug "time solve" advance_solution_by!(parent, child, dt)
+        if _child_failed(child)
+            parent.force_stepfail = true
+            return
+        end
+        backward_sync_subintegrator!(parent, child, idxs, sync)
+    end
+
+    # If B contaminated the solution (e.g. NaN), skip the second A step
+    # and let the outer check_error! detect instability on the next iteration.
+    if !all(isfinite, parent.u)
+        _force_set_time!(children[1], children[2].t)
+        return
+    end
+
+    # A(dt/2)
+    let child = children[1]
+        idxs = parent.child_solution_indices[1]
+        sync = parent.child_synchronizers[1]
+        @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync)
+        @timeit_debug "time solve" advance_solution_by!(parent, child, half_dt)
+        if _child_failed(child)
+            parent.force_stepfail = true
+            return
+        end
+        backward_sync_subintegrator!(parent, child, idxs, sync)
+    end
+
+    # Snap child 1's accumulated time (two half-steps) to match child 2's
+    # (one full step) to prevent floating-point drift.
+    try_snap_children_to_tstop!(children[1], children[2].t)
+end
diff --git a/src/utils.jl b/src/utils.jl
index a103e69..d06e44e 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -162,6 +162,14 @@ function OrdinaryDiffEqCore.fix_dt_at_bounds!(integrator::AnySplitIntegrator)
     return nothing
 end
 
+_force_set_time!(child::DEIntegrator, t) = (child.t = t)
+function _force_set_time!(child::SplitSubIntegrator, t)
+    child.t = t
+    for sub in child.child_subintegrators
+        _force_set_time!(sub, t)
+    end
+end
+
 # Check time-step information consistency
 validate_time_point(integrator::AnySplitIntegrator) = validate_time_point(integrator, integrator.child_subintegrators)
 function validate_time_point(parent, child::SplitSubIntegrator)
diff --git a/test/operator_splitting_api.jl b/test/operator_splitting_api.jl
index 2f11a30..1682cd1 100644
--- a/test/operator_splitting_api.jl
+++ b/test/operator_splitting_api.jl
@@ -281,6 +281,75 @@ end
         end
     end
 
+    @testset "StrangMarchuk reinit and convergence" begin
+        @testset "$tstepper" for (prob, tstepper) in (
+                (prob1a, StrangMarchuk((Euler(), Euler()))),
+                (prob1a, StrangMarchuk((Tsit5(), Euler()))),
+                (prob1a, StrangMarchuk((Euler(), Tsit5()))),
+                (prob1a, StrangMarchuk((Tsit5(), Tsit5()))),
+                (prob1b, StrangMarchuk((Euler(), Euler()))),
+                (prob1b, StrangMarchuk((Tsit5(), Euler()))),
+                (prob1b, StrangMarchuk((Euler(), Tsit5()))),
+                (prob1b, StrangMarchuk((Tsit5(), Tsit5()))),
+                (prob2, StrangMarchuk((Euler(), StrangMarchuk((Euler(), Euler()))))),
+                (prob2, StrangMarchuk((Tsit5(), StrangMarchuk((Tsit5(), Tsit5()))))),
+            )
+            integrator = DiffEqBase.init(
+                prob, tstepper, dt = dt, verbose = true, alias_u0 = false, adaptive = false
+            )
+            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default
+
+            sub1 = integrator.child_subintegrators[1]
+            sub2 = integrator.child_subintegrators[2]
+
+            DiffEqBase.solve!(integrator)
+            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success
+            ufinal = copy(integrator.u)
+            @test isapprox(ufinal, trueu, atol = 1.0e-6)
+            @test integrator.t ≈ tspan[2]
+            @test integrator.dtcache ≈ dt
+            @test integrator.iter == nsteps
+
+            @test sub1.t ≈ tspan[2]
+            @test sub1.iter == 2 * nsteps
+
+            @test sub2.t ≈ tspan[2]
+            @test sub2.iter == nsteps
+
+            DiffEqBase.reinit!(integrator; dt = dt)
+            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default
+            for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, tspan[1]:5.0:tspan[2])
+            end
+            @test isapprox(ufinal, integrator.u, atol = 1.0e-12)
+            @test integrator.t ≈ tspan[2]
+            @test integrator.dtcache ≈ dt
+            @test integrator.iter == nsteps
+
+            DiffEqBase.reinit!(integrator; dt = dt)
+            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default
+            for (uprev, tprev, u, t) in DiffEqBase.intervals(integrator)
+            end
+            @test isapprox(ufinal, integrator.u, atol = 1.0e-12)
+            @test integrator.t ≈ tspan[2]
+            @test integrator.dtcache ≈ dt
+            @test integrator.iter == nsteps
+
+            DiffEqBase.reinit!(integrator; dt = dt)
+            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default
+            DiffEqBase.solve!(integrator)
+            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success
+            @test integrator.t ≈ tspan[2]
+            @test integrator.dtcache ≈ dt
+            @test integrator.iter == nsteps
+
+            @test sub1.t ≈ tspan[2]
+            @test sub1.iter == 2 * nsteps
+
+            @test sub2.t ≈ tspan[2]
+            @test sub2.iter == nsteps
+        end
+    end
+
     @testset "Instability detection" begin
         dt = 0.01π
 
@@ -296,7 +365,7 @@ end
         fsplit_NaN = GenericSplitFunction((f1, f_NaN), (f1dofs, f3dofs))
         prob_NaN = OperatorSplittingProblem(fsplit_NaN, u0, tspan)
 
-        for TimeStepperType in (LieTrotterGodunov,)
+        for TimeStepperType in (LieTrotterGodunov, StrangMarchuk)
             @testset "Solver type $TimeStepperType | $tstepper" for tstepper in (
                     TimeStepperType((Euler(), Euler())),
                     TimeStepperType((Tsit5(), Euler())),

From 5c97963b351471a6e06c3b480264fda1871b12ee Mon Sep 17 00:00:00 2001
From: Orjan Ameye <orjan.ameye@hotmail.com>
Date: Sat, 9 May 2026 08:30:21 +0200
Subject: [PATCH 2/8] Integrate StrangMarchuk into existing test loops

Remove separate StrangMarchuk testset and instead add it to the main
convergence and adaptive test loops alongside LieTrotterGodunov. This
gives StrangMarchuk the same coverage: all nested splitting
combinations, FakeAdaptive wrapper, and adaptive mode tests.
---
 test/operator_splitting_api.jl | 84 +++++-----------------------------
 1 file changed, 11 insertions(+), 73 deletions(-)

diff --git a/test/operator_splitting_api.jl b/test/operator_splitting_api.jl
index 1682cd1..8eb45cf 100644
--- a/test/operator_splitting_api.jl
+++ b/test/operator_splitting_api.jl
@@ -126,6 +126,7 @@ end
 end
 
 FakeAdaptiveLTG(inner) = FakeAdaptiveAlgorithm(LieTrotterGodunov(inner))
+FakeAdaptiveSM(inner) = FakeAdaptiveAlgorithm(StrangMarchuk(inner))
 
 function Base.show(io::IO, alg::FakeAdaptiveAlgorithm)
     print(io, "FAKE (")
@@ -133,6 +134,11 @@ function Base.show(io::IO, alg::FakeAdaptiveAlgorithm)
     return print(io, ")")
 end
 
+# StrangMarchuk steps child 1 twice per outer step (two half-steps).
+_sub1_iter_factor(::LieTrotterGodunov) = 1
+_sub1_iter_factor(::StrangMarchuk) = 2
+_sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg)
+
 
 # ---------------------------------------------------------------------------
 # Tests
@@ -161,7 +167,7 @@ end
 
     nsteps = ceil(Int, (tspan[2] - tspan[1]) / dt)
 
-    for TimeStepperType in (LieTrotterGodunov, FakeAdaptiveLTG)
+    for TimeStepperType in (LieTrotterGodunov, FakeAdaptiveLTG, StrangMarchuk, FakeAdaptiveSM)
         @testset "$tstepper" for (prob, tstepper) in (
                 (prob1a, TimeStepperType((Euler(), Euler()))),
                 (prob1a, TimeStepperType((Tsit5(), Euler()))),
@@ -185,6 +191,7 @@ end
 
             sub1 = integrator.child_subintegrators[1]
             sub2 = integrator.child_subintegrators[2]
+            expected_sub1_iters = _sub1_iter_factor(tstepper) * nsteps
 
             DiffEqBase.solve!(integrator)
             @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success
@@ -195,7 +202,7 @@ end
             @test integrator.iter == nsteps
 
             @test sub1.t ≈ tspan[2]
-            @test sub1.iter == nsteps
+            @test sub1.iter == expected_sub1_iters
 
             @test sub2.t ≈ tspan[2]
             @test sub2.iter == nsteps
@@ -227,14 +234,14 @@ end
             @test integrator.iter == nsteps
 
             @test sub1.t ≈ tspan[2]
-            @test sub1.iter == nsteps
+            @test sub1.iter == expected_sub1_iters
 
             @test sub2.t ≈ tspan[2]
             @test sub2.iter == nsteps
         end
     end
 
-    for TimeStepperType in (FakeAdaptiveLTG,)
+    for TimeStepperType in (FakeAdaptiveLTG, FakeAdaptiveSM)
         @testset "Adaptive solver type $TimeStepperType | $tstepper" for (prob, tstepper) in (
                 (prob1a, TimeStepperType((Tsit5(), Tsit5()))),
                 (prob2, TimeStepperType((Tsit5(), TimeStepperType((Tsit5(), Tsit5()))))),
@@ -281,75 +288,6 @@ end
         end
     end
 
-    @testset "StrangMarchuk reinit and convergence" begin
-        @testset "$tstepper" for (prob, tstepper) in (
-                (prob1a, StrangMarchuk((Euler(), Euler()))),
-                (prob1a, StrangMarchuk((Tsit5(), Euler()))),
-                (prob1a, StrangMarchuk((Euler(), Tsit5()))),
-                (prob1a, StrangMarchuk((Tsit5(), Tsit5()))),
-                (prob1b, StrangMarchuk((Euler(), Euler()))),
-                (prob1b, StrangMarchuk((Tsit5(), Euler()))),
-                (prob1b, StrangMarchuk((Euler(), Tsit5()))),
-                (prob1b, StrangMarchuk((Tsit5(), Tsit5()))),
-                (prob2, StrangMarchuk((Euler(), StrangMarchuk((Euler(), Euler()))))),
-                (prob2, StrangMarchuk((Tsit5(), StrangMarchuk((Tsit5(), Tsit5()))))),
-            )
-            integrator = DiffEqBase.init(
-                prob, tstepper, dt = dt, verbose = true, alias_u0 = false, adaptive = false
-            )
-            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default
-
-            sub1 = integrator.child_subintegrators[1]
-            sub2 = integrator.child_subintegrators[2]
-
-            DiffEqBase.solve!(integrator)
-            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success
-            ufinal = copy(integrator.u)
-            @test isapprox(ufinal, trueu, atol = 1.0e-6)
-            @test integrator.t ≈ tspan[2]
-            @test integrator.dtcache ≈ dt
-            @test integrator.iter == nsteps
-
-            @test sub1.t ≈ tspan[2]
-            @test sub1.iter == 2 * nsteps
-
-            @test sub2.t ≈ tspan[2]
-            @test sub2.iter == nsteps
-
-            DiffEqBase.reinit!(integrator; dt = dt)
-            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default
-            for (u, t) in DiffEqBase.TimeChoiceIterator(integrator, tspan[1]:5.0:tspan[2])
-            end
-            @test isapprox(ufinal, integrator.u, atol = 1.0e-12)
-            @test integrator.t ≈ tspan[2]
-            @test integrator.dtcache ≈ dt
-            @test integrator.iter == nsteps
-
-            DiffEqBase.reinit!(integrator; dt = dt)
-            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default
-            for (uprev, tprev, u, t) in DiffEqBase.intervals(integrator)
-            end
-            @test isapprox(ufinal, integrator.u, atol = 1.0e-12)
-            @test integrator.t ≈ tspan[2]
-            @test integrator.dtcache ≈ dt
-            @test integrator.iter == nsteps
-
-            DiffEqBase.reinit!(integrator; dt = dt)
-            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Default
-            DiffEqBase.solve!(integrator)
-            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success
-            @test integrator.t ≈ tspan[2]
-            @test integrator.dtcache ≈ dt
-            @test integrator.iter == nsteps
-
-            @test sub1.t ≈ tspan[2]
-            @test sub1.iter == 2 * nsteps
-
-            @test sub2.t ≈ tspan[2]
-            @test sub2.iter == nsteps
-        end
-    end
-
     @testset "Instability detection" begin
         dt = 0.01π
 

From b9f0c6face7e7da90f8787666ced8fca87cf990a Mon Sep 17 00:00:00 2001
From: Orjan Ameye <orjan.ameye@hotmail.com>
Date: Sat, 9 May 2026 08:49:57 +0200
Subject: [PATCH 3/8] Add review fixes: constructor guard, convergence test,
 precompilation, bibliography

- Validate StrangMarchuk requires exactly 2 inner algorithms
- Add convergence order test with non-commuting operators to verify
  first-order (LTG) vs second-order (SM) splitting accuracy
- Add StrangMarchuk precompilation workload
- Add Strang 1968 and Marchuk 1971 bibliography entries
- Remove accidental OrdinaryDiffEqTsit5 from [deps]
---
 docs/src/assets/references.bib | 20 +++++++++++++++
 src/precompilation.jl          | 14 ++++++++---
 src/solver.jl                  |  6 +++++
 test/operator_splitting_api.jl | 45 ++++++++++++++++++++++++++++++++++
 4 files changed, 82 insertions(+), 3 deletions(-)

diff --git a/docs/src/assets/references.bib b/docs/src/assets/references.bib
index 7184cca..3cb2732 100644
--- a/docs/src/assets/references.bib
+++ b/docs/src/assets/references.bib
@@ -32,3 +32,23 @@ @article{God:1959:dmn
   year={1959},
   publisher={Russian Academy of Sciences, Steklov Mathematical Institute}
 }
+
+@article{Str:1968:ccd,
+  title={On the construction and comparison of difference schemes},
+  author={Strang, Gilbert},
+  journal={SIAM Journal on Numerical Analysis},
+  volume={5},
+  number={3},
+  pages={506--517},
+  year={1968},
+  publisher={SIAM}
+}
+
+@incollection{Mar:1971:tsm,
+  title={On the theory of the splitting-up method},
+  author={Marchuk, Guri Ivanovich},
+  booktitle={Numerical Solution of Partial Differential Equations-{II}},
+  pages={469--500},
+  year={1971},
+  publisher={Academic Press}
+}
diff --git a/src/precompilation.jl b/src/precompilation.jl
index f638e74..df96482 100644
--- a/src/precompilation.jl
+++ b/src/precompilation.jl
@@ -34,10 +34,18 @@ end
     fsplit = GenericSplitFunction((f1, fsplitinner), (f1dofs, [1, 2, 3]))
 
     prob = OperatorSplittingProblem(fsplit, u0, tspan)
-    tstepper = LieTrotterGodunov((Euler(), LieTrotterGodunov((Euler(), Euler()))))
 
-    # Precompile init and a few steps
-    integrator = DiffEqBase.init(prob, tstepper, dt = 0.01, verbose = false)
+    # Precompile LieTrotterGodunov
+    tstepper_ltg = LieTrotterGodunov((Euler(), LieTrotterGodunov((Euler(), Euler()))))
+    integrator = DiffEqBase.init(prob, tstepper_ltg, dt = 0.01, verbose = false)
     step!(integrator)
     solve!(integrator)
+
+    # Precompile StrangMarchuk
+    fsplit_sm = GenericSplitFunction((f1, f2), (f1dofs, f2dofs))
+    prob_sm = OperatorSplittingProblem(fsplit_sm, u0, tspan)
+    tstepper_sm = StrangMarchuk((Euler(), Euler()))
+    integrator_sm = DiffEqBase.init(prob_sm, tstepper_sm, dt = 0.01, verbose = false)
+    step!(integrator_sm)
+    solve!(integrator_sm)
 end
diff --git a/src/solver.jl b/src/solver.jl
index 1431a2a..7dde197 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -72,6 +72,12 @@ achieving second-order accuracy through symmetry.
 """
 struct StrangMarchuk{AlgTupleType} <: AbstractOperatorSplittingAlgorithm
     inner_algs::AlgTupleType
+    function StrangMarchuk(inner_algs::T) where {T <: Tuple}
+        length(inner_algs) == 2 || throw(
+            ArgumentError("StrangMarchuk requires exactly 2 inner algorithms, got $(length(inner_algs))")
+        )
+        return new{T}(inner_algs)
+    end
 end
 
 function Base.show(io::IO, alg::StrangMarchuk)
diff --git a/test/operator_splitting_api.jl b/test/operator_splitting_api.jl
index 8eb45cf..a0506d7 100644
--- a/test/operator_splitting_api.jl
+++ b/test/operator_splitting_api.jl
@@ -288,6 +288,51 @@ _sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg)
         end
     end
 
+    @testset "Convergence order" begin
+        # Use non-commuting operators so splitting error is non-zero.
+        # A = diag(-1,-2), B = [0 0.5; 0.5 0] have [A,B] ≠ 0.
+        function ode_conv_A(du, u, p, t)
+            du[1] = -u[1]
+            return du[2] = -2 * u[2]
+        end
+        function ode_conv_B(du, u, p, t)
+            du[1] = 0.5 * u[2]
+            return du[2] = 0.5 * u[1]
+        end
+        fA = ODEFunction(ode_conv_A)
+        fB = ODEFunction(ode_conv_B)
+
+        conv_tspan = (0.0, 1.0)
+        conv_u0 = [1.0, 1.0]
+        conv_trueu = exp(conv_tspan[2] * [-1.0 0.5; 0.5 -2.0]) * conv_u0
+
+        conv_dofs = [1, 2]
+        fsplit_conv = GenericSplitFunction((fA, fB), (conv_dofs, conv_dofs))
+        prob_conv = OperatorSplittingProblem(fsplit_conv, conv_u0, conv_tspan)
+
+        dts = [0.1, 0.05, 0.025]
+        for (TimeStepperType, expected_order) in (
+                (LieTrotterGodunov, 1),
+                (StrangMarchuk, 2),
+            )
+            @testset "$TimeStepperType (order $expected_order)" begin
+                errors = map(dts) do dt_i
+                    tstepper = TimeStepperType((Tsit5(), Tsit5()))
+                    integrator = DiffEqBase.init(
+                        prob_conv, tstepper, dt = dt_i, verbose = false,
+                        alias_u0 = false, adaptive = false
+                    )
+                    DiffEqBase.solve!(integrator)
+                    maximum(abs, integrator.u .- conv_trueu)
+                end
+                for i in 1:(length(errors) - 1)
+                    rate = log2(errors[i] / errors[i + 1])
+                    @test rate ≈ expected_order atol = 0.3
+                end
+            end
+        end
+    end
+
     @testset "Instability detection" begin
         dt = 0.01π
 

From 8cc1c0f0269436dbd59869afe6a904e4ce4c3733 Mon Sep 17 00:00:00 2001
From: Orjan Ameye <orjan.ameye@hotmail.com>
Date: Sat, 9 May 2026 09:04:28 +0200
Subject: [PATCH 4/8] Add StrangMarchuk to docs: API reference, usage example,
 theory section

---
 docs/src/api-reference/index.md     |  3 ++-
 docs/src/topics/time-integration.md | 42 +++++++++++++++++++++++++++++
 docs/src/usage/index.md             | 15 +++++++++++
 3 files changed, 59 insertions(+), 1 deletion(-)

diff --git a/docs/src/api-reference/index.md b/docs/src/api-reference/index.md
index 1e84973..7c630e8 100644
--- a/docs/src/api-reference/index.md
+++ b/docs/src/api-reference/index.md
@@ -12,8 +12,9 @@ OperatorSplittingProblem
 GenericSplitFunction
 ```
 
-## Solver
+## Solvers
 
 ```@docs
 LieTrotterGodunov
+StrangMarchuk
 ```
diff --git a/docs/src/topics/time-integration.md b/docs/src/topics/time-integration.md
index f6fccf1..4af9039 100644
--- a/docs/src/topics/time-integration.md
+++ b/docs/src/topics/time-integration.md
@@ -116,6 +116,48 @@ $n \in \mathbb{N}$ the following bound
 
 which implies stability of the scheme.
 
+### Strang-Marchuk Splitting
+
+A natural way to improve the accuracy of operator splitting is to symmetrize the
+scheme. The Strang-Marchuk splitting [Str:1968:ccd,Mar:1971:tsm](@cite) achieves
+second-order accuracy for two operators $F_1$ and $F_2$ by performing
+
+```math
+\begin{aligned}
+    \text{Solve} \quad d_t u^1(t) &= F_1(u^1(t), p, t) & & \quad \text{on} \; [t_0, t_0 + \Delta t/2] \; \text{with} \; u^1(t_0) = u_0 \\
+    \text{Solve} \quad d_t u^2(t) &= F_2(u^2(t), p, t) & & \quad \text{on} \; [t_0, t_0 + \Delta t] \; \text{with} \; u^2(t_0) = u^1(t_0 + \Delta t/2) \\
+    \text{Solve} \quad d_t u^3(t) &= F_1(u^3(t), p, t) & & \quad \text{on} \; [t_0 + \Delta t/2, t_0 + \Delta t] \; \text{with} \; u^3(t_0 + \Delta t/2) = u^2(t_0 + \Delta t)
+\end{aligned}
+```
+
+yielding $u(t_0 + \Delta t) \approx u^3(t_0 + \Delta t)$.
+
+### Analysis of Strang-Marchuk
+
+For two bounded linear operators $L_1$ and $L_2$ the Strang-Marchuk approximation
+reads
+
+```math
+\tilde{u}(t) = e^{L_1 t/2} \, e^{L_2 t} \, e^{L_1 t/2} \, u_0 \, .
+```
+
+Expanding the exponentials:
+
+```math
+\begin{aligned}
+e^{L_1 t/2} \, e^{L_2 t} \, e^{L_1 t/2}
+&= \bigl(I + \tfrac{t}{2}L_1 + \tfrac{t^2}{8}L_1^2 + \cdots\bigr)
+   \bigl(I + t L_2 + \tfrac{t^2}{2}L_2^2 + \cdots\bigr)
+   \bigl(I + \tfrac{t}{2}L_1 + \tfrac{t^2}{8}L_1^2 + \cdots\bigr) \\
+&= I + t(L_1 + L_2) + \tfrac{t^2}{2}(L_1 + L_2)^2 + O(t^3)
+\end{aligned}
+```
+
+which matches the Taylor expansion of $e^{(L_1+L_2)t}$ through the $t^2$ term.
+The symmetry of the scheme causes the first-order commutator term
+$[L_1, L_2] = L_1 L_2 - L_2 L_1$ to cancel, leaving a local truncation error
+of $O(t^3)$ and hence second-order global accuracy.
+
 ## References
 
 ```@bibliography
diff --git a/docs/src/usage/index.md b/docs/src/usage/index.md
index 9dc456c..e95cecd 100644
--- a/docs/src/usage/index.md
+++ b/docs/src/usage/index.md
@@ -49,3 +49,18 @@ for (u, t) in TimeChoiceIterator(integrator, 0.0:0.5:1.0)
     @show t, u
 end
 ```
+
+For second-order accuracy, use the `StrangMarchuk` algorithm instead.
+It requires exactly two subproblems and performs the symmetric
+A(Δt/2) → B(Δt) → A(Δt/2) splitting:
+
+```julia
+alg = StrangMarchuk(
+    (Euler(), Euler())
+)
+
+integrator = init(prob, alg, dt = 0.1)
+for (u, t) in TimeChoiceIterator(integrator, 0.0:0.5:1.0)
+    @show t, u
+end
+```

From d3bfd76e10028bcfc7d18e644071d0f1360558d5 Mon Sep 17 00:00:00 2001
From: Orjan Ameye <orjan.ameye@hotmail.com>
Date: Sun, 10 May 2026 07:40:46 +0200
Subject: [PATCH 5/8] Generalize StrangMarchuk to N operators with palindromic
 scheme

Replace the hardcoded 2-operator implementation with a general
N-operator palindromic sweep using @unroll, matching LTG's style.
Add 3-operator test verifying iteration counts and solution accuracy.
Update docs and docstrings for the general scheme.
---
 docs/src/topics/time-integration.md | 15 +++--
 docs/src/usage/index.md             |  4 +-
 src/solver.jl                       | 97 ++++++++++++++++-------------
 test/operator_splitting_api.jl      | 33 ++++++++++
 4 files changed, 98 insertions(+), 51 deletions(-)

diff --git a/docs/src/topics/time-integration.md b/docs/src/topics/time-integration.md
index 4af9039..9ecb7f5 100644
--- a/docs/src/topics/time-integration.md
+++ b/docs/src/topics/time-integration.md
@@ -120,7 +120,13 @@ which implies stability of the scheme.
 
 A natural way to improve the accuracy of operator splitting is to symmetrize the
 scheme. The Strang-Marchuk splitting [Str:1968:ccd,Mar:1971:tsm](@cite) achieves
-second-order accuracy for two operators $F_1$ and $F_2$ by performing
+second-order accuracy for $N$ operators by performing a palindromic sweep
+
+```math
+F_1(\Delta t/2) \to \cdots \to F_{N-1}(\Delta t/2) \to F_N(\Delta t) \to F_{N-1}(\Delta t/2) \to \cdots \to F_1(\Delta t/2)
+```
+
+More formally, for the simplest case of two operators $F_1$ and $F_2$
 
 ```math
 \begin{aligned}
@@ -134,8 +140,8 @@ yielding $u(t_0 + \Delta t) \approx u^3(t_0 + \Delta t)$.
 
 ### Analysis of Strang-Marchuk
 
-For two bounded linear operators $L_1$ and $L_2$ the Strang-Marchuk approximation
-reads
+We show the second-order accuracy for two bounded linear operators $L_1$ and
+$L_2$. The Strang-Marchuk approximation reads
 
 ```math
 \tilde{u}(t) = e^{L_1 t/2} \, e^{L_2 t} \, e^{L_1 t/2} \, u_0 \, .
@@ -156,7 +162,8 @@ e^{L_1 t/2} \, e^{L_2 t} \, e^{L_1 t/2}
 which matches the Taylor expansion of $e^{(L_1+L_2)t}$ through the $t^2$ term.
 The symmetry of the scheme causes the first-order commutator term
 $[L_1, L_2] = L_1 L_2 - L_2 L_1$ to cancel, leaving a local truncation error
-of $O(t^3)$ and hence second-order global accuracy.
+of $O(t^3)$ and hence second-order global accuracy. The same argument extends to
+the general $N$-operator palindromic scheme.
 
 ## References
 
diff --git a/docs/src/usage/index.md b/docs/src/usage/index.md
index e95cecd..4e3d0e4 100644
--- a/docs/src/usage/index.md
+++ b/docs/src/usage/index.md
@@ -51,8 +51,8 @@ end
 ```
 
 For second-order accuracy, use the `StrangMarchuk` algorithm instead.
-It requires exactly two subproblems and performs the symmetric
-A(Δt/2) → B(Δt) → A(Δt/2) splitting:
+It performs the symmetric palindromic splitting
+A₁(Δt/2) → … → Aₙ(Δt) → … → A₁(Δt/2):
 
 ```julia
 alg = StrangMarchuk(
diff --git a/src/solver.jl b/src/solver.jl
index 7dde197..c4dd938 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -63,21 +63,17 @@ end
 """
     StrangMarchuk <: AbstractOperatorSplittingAlgorithm
 
-Second-order symmetric operator splitting algorithm attributed to
+Second-order symmetric (palindromic) operator splitting algorithm attributed to
 [Str:1968:ccd,Mar:1971:tsm](@cite).
 
-For two operators ``A`` and ``B`` the scheme performs
-``A(\\Delta t/2) \\to B(\\Delta t) \\to A(\\Delta t/2)``,
+For ``N`` operators the scheme performs
+
+``A_1(\\Delta t/2) \\to \\cdots \\to A_{N-1}(\\Delta t/2) \\to A_N(\\Delta t) \\to A_{N-1}(\\Delta t/2) \\to \\cdots \\to A_1(\\Delta t/2)``
+
 achieving second-order accuracy through symmetry.
 """
 struct StrangMarchuk{AlgTupleType} <: AbstractOperatorSplittingAlgorithm
-    inner_algs::AlgTupleType
-    function StrangMarchuk(inner_algs::T) where {T <: Tuple}
-        length(inner_algs) == 2 || throw(
-            ArgumentError("StrangMarchuk requires exactly 2 inner algorithms, got $(length(inner_algs))")
-        )
-        return new{T}(inner_algs)
-    end
+    inner_algs::AlgTupleType # Tuple of timesteppers for inner problems
 end
 
 function Base.show(io::IO, alg::StrangMarchuk)
@@ -102,61 +98,72 @@ function init_cache(
     return StrangMarchukCache(u, uprev)
 end
 
-function _perform_step!(
-        parent,
-        children::Tuple,
-        cache::StrangMarchukCache,
-        dt
-    )
-    half_dt = dt / 2
+# Forward pass: A₁(dt/2) → … → Aₙ₋₁(dt/2) → Aₙ(dt)
+@unroll function _sm_forward_pass!(parent, children::Tuple, half_dt, dt)
+    N = length(children)
+    i = 0
+    @unroll for child in children
+        i += 1
+        step_dt = i < N ? half_dt : dt
+
+        idxs = parent.child_solution_indices[i]
+        sync = parent.child_synchronizers[i]
 
-    # A(dt/2)
-    let child = children[1]
-        idxs = parent.child_solution_indices[1]
-        sync = parent.child_synchronizers[1]
         @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync)
-        @timeit_debug "time solve" advance_solution_by!(parent, child, half_dt)
+        @timeit_debug "time solve" advance_solution_by!(parent, child, step_dt)
         if _child_failed(child)
             parent.force_stepfail = true
             return
         end
+
         backward_sync_subintegrator!(parent, child, idxs, sync)
     end
+end
+
+# Reverse pass: Aₙ₋₁(dt/2) → … → A₁(dt/2)
+@unroll function _sm_reverse_pass!(parent, rev_front::Tuple, half_dt, N)
+    j = 0
+    @unroll for child in rev_front
+        j += 1
+        i = N - j
+
+        idxs = parent.child_solution_indices[i]
+        sync = parent.child_synchronizers[i]
 
-    # B(dt)
-    let child = children[2]
-        idxs = parent.child_solution_indices[2]
-        sync = parent.child_synchronizers[2]
         @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync)
-        @timeit_debug "time solve" advance_solution_by!(parent, child, dt)
+        @timeit_debug "time solve" advance_solution_by!(parent, child, half_dt)
         if _child_failed(child)
             parent.force_stepfail = true
             return
         end
+
         backward_sync_subintegrator!(parent, child, idxs, sync)
     end
+end
+
+function _perform_step!(
+        parent,
+        children::Tuple,
+        cache::StrangMarchukCache,
+        dt
+    )
+    half_dt = dt / 2
+    N = length(children)
+
+    _sm_forward_pass!(parent, children, half_dt, dt)
+    parent.force_stepfail && return
 
-    # If B contaminated the solution (e.g. NaN), skip the second A step
-    # and let the outer check_error! detect instability on the next iteration.
     if !all(isfinite, parent.u)
-        _force_set_time!(children[1], children[2].t)
+        for i in 1:(N - 1)
+            _force_set_time!(children[i], children[N].t)
+        end
         return
     end
 
-    # A(dt/2)
-    let child = children[1]
-        idxs = parent.child_solution_indices[1]
-        sync = parent.child_synchronizers[1]
-        @timeit_debug "sync ->" forward_sync_subintegrator!(parent, child, idxs, sync)
-        @timeit_debug "time solve" advance_solution_by!(parent, child, half_dt)
-        if _child_failed(child)
-            parent.force_stepfail = true
-            return
-        end
-        backward_sync_subintegrator!(parent, child, idxs, sync)
-    end
+    _sm_reverse_pass!(parent, reverse(children[1:(end - 1)]), half_dt, N)
+    parent.force_stepfail && return
 
-    # Snap child 1's accumulated time (two half-steps) to match child 2's
-    # (one full step) to prevent floating-point drift.
-    try_snap_children_to_tstop!(children[1], children[2].t)
+    for i in 1:(N - 1)
+        try_snap_children_to_tstop!(children[i], children[N].t)
+    end
 end
diff --git a/test/operator_splitting_api.jl b/test/operator_splitting_api.jl
index a0506d7..7cabb90 100644
--- a/test/operator_splitting_api.jl
+++ b/test/operator_splitting_api.jl
@@ -288,6 +288,39 @@ _sub1_iter_factor(alg::FakeAdaptiveAlgorithm) = _sub1_iter_factor(alg.alg)
         end
     end
 
+    @testset "StrangMarchuk with 3 operators" begin
+        dt = 0.01π
+        # f1 + f3 + f3 = f1 + f2, so the reference solution is the same trueu.
+        f1dofs = [1, 2, 3]
+        f3dofs = [1, 3]
+        fsplit3 = GenericSplitFunction((f1, f3, f3), (f1dofs, f3dofs, f3dofs))
+        prob3 = OperatorSplittingProblem(fsplit3, u0, tspan)
+        nsteps = ceil(Int, (tspan[2] - tspan[1]) / dt)
+
+        @testset "$tstepper" for tstepper in (
+                StrangMarchuk((Euler(), Euler(), Euler())),
+                StrangMarchuk((Tsit5(), Euler(), Tsit5())),
+                StrangMarchuk((Tsit5(), Tsit5(), Tsit5())),
+            )
+            integrator = DiffEqBase.init(
+                prob3, tstepper, dt = dt, verbose = true, alias_u0 = false, adaptive = false
+            )
+            DiffEqBase.solve!(integrator)
+            @test integrator.sol.retcode == DiffEqBase.ReturnCode.Success
+            @test isapprox(integrator.u, trueu, atol = 1.0e-6)
+            @test integrator.t ≈ tspan[2]
+            @test integrator.iter == nsteps
+
+            sub1 = integrator.child_subintegrators[1]
+            sub2 = integrator.child_subintegrators[2]
+            sub3 = integrator.child_subintegrators[3]
+            # Palindromic: children 1 & 2 get two half-steps, child 3 gets one full step
+            @test sub1.iter == 2 * nsteps
+            @test sub2.iter == 2 * nsteps
+            @test sub3.iter == nsteps
+        end
+    end
+
     @testset "Convergence order" begin
         # Use non-commuting operators so splitting error is non-zero.
         # A = diag(-1,-2), B = [0 0.5; 0.5 0] have [A,B] ≠ 0.

From 41875ac6b51950a43e6033c66981f5eb201c13cc Mon Sep 17 00:00:00 2001
From: Orjan Ameye <orjan.ameye@hotmail.com>
Date: Sun, 10 May 2026 07:41:18 +0200
Subject: [PATCH 6/8] format

---
 src/solver.jl | 1 +
 src/utils.jl  | 1 +
 2 files changed, 2 insertions(+)

diff --git a/src/solver.jl b/src/solver.jl
index c4dd938..43b4aef 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -166,4 +166,5 @@ function _perform_step!(
     for i in 1:(N - 1)
         try_snap_children_to_tstop!(children[i], children[N].t)
     end
+    return
 end
diff --git a/src/utils.jl b/src/utils.jl
index d06e44e..ebf381b 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -168,6 +168,7 @@ function _force_set_time!(child::SplitSubIntegrator, t)
     for sub in child.child_subintegrators
         _force_set_time!(sub, t)
     end
+    return
 end
 
 # Check time-step information consistency

From 45132b4f67bfdc0396f8ff1e20ff38768e693f31 Mon Sep 17 00:00:00 2001
From: Orjan Ameye <orjan.ameye@hotmail.com>
Date: Sun, 10 May 2026 21:54:03 +0200
Subject: [PATCH 7/8] adress review

---
 src/integrator.jl |  5 +++++
 src/solver.jl     | 11 ++++++++---
 src/utils.jl      |  6 ++++++
 3 files changed, 19 insertions(+), 3 deletions(-)

diff --git a/src/integrator.jl b/src/integrator.jl
index 01c5e23..37e9993 100644
--- a/src/integrator.jl
+++ b/src/integrator.jl
@@ -605,6 +605,11 @@ function step_footer!(integrator::AnySplitIntegrator)
         integrator.last_step_failed = false
         integrator.tprev = integrator.t
         integrator.t = fixed_t_for_floatingpoint_error!(integrator, ttmp)
+        # Children that step with subdivided dt (e.g. StrangMarchuk's `dt/2`
+        # halves) accumulate ulp-level drift from the parent's exact `t`.
+        # Re-anchor children to the parent's canonical time here so the
+        # drift cannot accumulate across outer steps.
+        try_snap_children_to_tstop!.(integrator.child_subintegrators, integrator.t)
         step_accept_controller!(integrator)
     elseif integrator.force_stepfail
         if isadaptive(integrator)
diff --git a/src/solver.jl b/src/solver.jl
index 43b4aef..3ff3acd 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -153,6 +153,14 @@ function _perform_step!(
     _sm_forward_pass!(parent, children, half_dt, dt)
     parent.force_stepfail && return
 
+    # If the forward pass produced NaN/Inf in `parent.u`, skip the reverse
+    # pass: the next reverse-pass `forward_sync_subintegrator!` would copy
+    # the unstable state into a child whose inner integrator then detects
+    # it during `step!` and hard-throws via the leaf `advance_solution_by!`.
+    # Returning here lets the outer `check_error!` see the unstable `u` and
+    # finish with `Unstable` / `DtNaN` retcode on the next iteration. The
+    # non-master children only got `dt/2`; force their `t` to the master's
+    # so `validate_time_point` doesn't assert on the half-step gap.
     if !all(isfinite, parent.u)
         for i in 1:(N - 1)
             _force_set_time!(children[i], children[N].t)
@@ -163,8 +171,5 @@ function _perform_step!(
     _sm_reverse_pass!(parent, reverse(children[1:(end - 1)]), half_dt, N)
     parent.force_stepfail && return
 
-    for i in 1:(N - 1)
-        try_snap_children_to_tstop!(children[i], children[N].t)
-    end
     return
 end
diff --git a/src/utils.jl b/src/utils.jl
index ebf381b..f1f2867 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -162,6 +162,12 @@ function OrdinaryDiffEqCore.fix_dt_at_bounds!(integrator::AnySplitIntegrator)
     return nothing
 end
 
+# Force-set the local time of an integrator (and recursively its children).
+# Used on the StrangMarchuk NaN-bailout path where a non-master child has only
+# finished half the outer step and we need to skip its remaining sub-steps:
+# we still have to advance its `t` to match the master so the next outer
+# iteration's `check_error!` can see the unstable `u` and report the failure
+# without `validate_time_point` asserting on the half-step gap.
 _force_set_time!(child::DEIntegrator, t) = (child.t = t)
 function _force_set_time!(child::SplitSubIntegrator, t)
     child.t = t

From 25fd68305101a9b94ca95ce24e76934c8e1e8505 Mon Sep 17 00:00:00 2001
From: Orjan Ameye <orjan.ameye@hotmail.com>
Date: Wed, 13 May 2026 20:47:25 +0200
Subject: [PATCH 8/8] validate_time_point in should_accept_step check

---
 src/integrator.jl | 16 +---------------
 src/solver.jl     | 18 +-----------------
 src/utils.jl      | 15 ---------------
 3 files changed, 2 insertions(+), 47 deletions(-)

diff --git a/src/integrator.jl b/src/integrator.jl
index 37e9993..2f9a288 100644
--- a/src/integrator.jl
+++ b/src/integrator.jl
@@ -611,6 +611,7 @@ function step_footer!(integrator::AnySplitIntegrator)
         # drift cannot accumulate across outer steps.
         try_snap_children_to_tstop!.(integrator.child_subintegrators, integrator.t)
         step_accept_controller!(integrator)
+        validate_time_point(integrator)
     elseif integrator.force_stepfail
         if isadaptive(integrator)
             step_reject_controller!(integrator)
@@ -622,7 +623,6 @@ function step_footer!(integrator::AnySplitIntegrator)
         end
         integrator.last_step_failed = true
     end
-    validate_time_point(integrator)
     return nothing
 end
 
@@ -912,26 +912,12 @@ function advance_solution_by!(
         dt
     )
     SciMLBase.step!(sub, dt, true)
-
-    # Unrecoverable failure: error immediately regardless of adaptive/non-adaptive
-    if !SciMLBase.successful_retcode(sub.status.retcode) &&
-            sub.status.retcode != ReturnCode.Default
-        error("Inner integrator failed unrecoverably with retcode \
-               $(sub.status.retcode) at t=$(child.t). Aborting.")
-    end
     return nothing
 end
 
 # Leaf disptach
 function advance_solution_by!(outer::AnySplitIntegrator, child::DEIntegrator, dt)
     SciMLBase.step!(child, dt, true)
-
-    # Unrecoverable failure: error immediately regardless of adaptive/non-adaptive
-    if !SciMLBase.successful_retcode(child.sol.retcode) &&
-            child.sol.retcode != ReturnCode.Default
-        error("Inner integrator failed unrecoverably with retcode \
-               $(child.sol.retcode) at t=$(child.t). Aborting.")
-    end
     return nothing
 end
 
diff --git a/src/solver.jl b/src/solver.jl
index 3ff3acd..a72c62f 100644
--- a/src/solver.jl
+++ b/src/solver.jl
@@ -148,27 +148,11 @@ function _perform_step!(
         dt
     )
     half_dt = dt / 2
-    N = length(children)
 
     _sm_forward_pass!(parent, children, half_dt, dt)
     parent.force_stepfail && return
 
-    # If the forward pass produced NaN/Inf in `parent.u`, skip the reverse
-    # pass: the next reverse-pass `forward_sync_subintegrator!` would copy
-    # the unstable state into a child whose inner integrator then detects
-    # it during `step!` and hard-throws via the leaf `advance_solution_by!`.
-    # Returning here lets the outer `check_error!` see the unstable `u` and
-    # finish with `Unstable` / `DtNaN` retcode on the next iteration. The
-    # non-master children only got `dt/2`; force their `t` to the master's
-    # so `validate_time_point` doesn't assert on the half-step gap.
-    if !all(isfinite, parent.u)
-        for i in 1:(N - 1)
-            _force_set_time!(children[i], children[N].t)
-        end
-        return
-    end
-
-    _sm_reverse_pass!(parent, reverse(children[1:(end - 1)]), half_dt, N)
+    _sm_reverse_pass!(parent, reverse(children[1:(end - 1)]), half_dt, length(children))
     parent.force_stepfail && return
 
     return
diff --git a/src/utils.jl b/src/utils.jl
index f1f2867..a103e69 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -162,21 +162,6 @@ function OrdinaryDiffEqCore.fix_dt_at_bounds!(integrator::AnySplitIntegrator)
     return nothing
 end
 
-# Force-set the local time of an integrator (and recursively its children).
-# Used on the StrangMarchuk NaN-bailout path where a non-master child has only
-# finished half the outer step and we need to skip its remaining sub-steps:
-# we still have to advance its `t` to match the master so the next outer
-# iteration's `check_error!` can see the unstable `u` and report the failure
-# without `validate_time_point` asserting on the half-step gap.
-_force_set_time!(child::DEIntegrator, t) = (child.t = t)
-function _force_set_time!(child::SplitSubIntegrator, t)
-    child.t = t
-    for sub in child.child_subintegrators
-        _force_set_time!(sub, t)
-    end
-    return
-end
-
 # Check time-step information consistency
 validate_time_point(integrator::AnySplitIntegrator) = validate_time_point(integrator, integrator.child_subintegrators)
 function validate_time_point(parent, child::SplitSubIntegrator)