diff --git a/Project.toml b/Project.toml
index f0c2129..a83ae67 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,7 +27,7 @@ OrdinaryDiffEqTsit5 = "1.1.0, 2"
 PrecompileTools = "1.1"
 RecursiveArrayTools = "3.39.0, 4"
 SafeTestsets = "0.1.0"
-SciMLBase = "2.77.0, 3.1"
+SciMLBase = "3.36"
 SciMLIterators = "1"
 SciMLTesting = "2.1"
 SymbolicIndexingInterface = "0.3.36"
diff --git a/docs/src/devdocs/index.md b/docs/src/devdocs/index.md
index 82708fc..30387a5 100644
--- a/docs/src/devdocs/index.md
+++ b/docs/src/devdocs/index.md
@@ -108,3 +108,42 @@ end
     # Done :)
 end
 ```
+
+## Dense output
+
+Saving, `saveat` and continuous callback root-finding all go through a single hook,
+so an algorithm only has to describe how to interpolate *within one of its steps*:
+
+```julia
+function OrdinaryDiffEqOperatorSplitting.splitting_interpolant(
+        integrator, cache::MySimpleFirstOrderCache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}}
+) where {D}
+    # Θ = (t - tprev) / dt is the step-local coordinate, y₀ = u(tprev), y₁ = u(t).
+    # ...
+end
+
+function OrdinaryDiffEqOperatorSplitting.splitting_interpolant!(
+        out, integrator, cache::MySimpleFirstOrderCache, Θ, dt, y₀, y₁, idxs,
+        ::Type{Val{D}}
+) where {D}
+    # In-place variant; must return `out`.
+end
+```
+
+Both fall back to linear interpolation for any `AbstractOperatorSplittingCache`, so
+implementing them is optional -- but with the fallback, saved output is first order
+even for a second-order scheme, and continuous callback event times are the exact
+roots of a straight chord across the step. One method improves both.
+
+The fallback is linear for a structural reason, and it is the same reason saving and
+callbacks live on the outer integrator alone: a splitting step advances its children
+sequentially over staggered subintervals, so a child's own interpolant describes a
+different sub-problem over a different interval, and they do not compose into an
+approximation of the split solution. Only the step endpoints, which the outer
+integrator owns, are states of the full split system -- an inner split is a stage,
+not a step.
+
+When implementing this, note that `Θ` is derived from
+`integrator.t - integrator.tprev`, **not** from `integrator.dt`: once a step is
+accepted, `step_accept_controller!` has already replaced `dt` with the step size
+proposed for the *next* step.
diff --git a/docs/src/usage/index.md b/docs/src/usage/index.md
index 94e2bcd..277ed4d 100644
--- a/docs/src/usage/index.md
+++ b/docs/src/usage/index.md
@@ -42,8 +42,9 @@ alg = LieTrotterGodunov(
     (Euler(), Euler())
 )
 
-# Right now OrdinaryDiffEqOperatorSplitting.jl does not implement the SciML solution interface,
-# but we can obtain intermediate solutions via the iterator interface.
+# OrdinaryDiffEqOperatorSplitting.jl implements only part of the SciML solution
+# interface (see "Saving and interpolation" below); intermediate solutions are most
+# directly obtained via the iterator interface.
 integrator = init(prob, alg, dt = 0.1)
 for (u, t) in TimeChoiceIterator(integrator, 0.0:0.5:1.0)
     @show t, u
@@ -65,6 +66,90 @@ for (u, t) in TimeChoiceIterator(integrator, 0.0:0.5:1.0)
 end
 ```
 
+## Saving and interpolation
+
+`solve!` returns a solution whose `t`/`u` are filled as the integration proceeds, so
+the usual saving keywords work and `sol(t)` interpolates between saved points:
+
+```julia
+integrator = init(prob, alg; dt = 0.1, saveat = 0.25)
+sol = solve!(integrator)
+
+sol.t          # [0.0, 0.25, 0.5, 0.75, 1.0]
+sol(0.3)       # interpolated between the saved points
+```
+
+The supported keywords are `saveat`, `save_everystep` (default `false`),
+`save_start`, `save_end` and `save_on`. They apply to the **outer** integrator only:
+the inner splits are stages rather than steps, so their intermediate states do not
+approximate the split solution at any time point.
+
+`saveat` points that fall strictly inside a step are filled from the step's
+interpolant, so requesting output never changes the sequence of steps and therefore
+never changes the splitting error.
+
+!!! note "The interpolant is first order"
+    A splitting step advances its children sequentially over staggered
+    subintervals, so their individual interpolants do not compose into an
+    approximation of the split solution. The dense output available at the outer
+    level is therefore linear between the step endpoints -- exact for the state a
+    `LieTrotterGodunov` step produces, but only first order for the second-order
+    schemes. If you need output that is accurate to the order of the scheme, pass
+    the times as `tstops` so that the integrator lands on them exactly:
+
+    ```julia
+    sol = solve!(init(prob, alg; dt = 0.1, tstops = [0.25, 0.5], saveat = [0.25, 0.5]))
+    ```
+
+    `dense = true` is accepted but has no effect: `sol(t)` already reproduces the
+    integrator's own interpolant from the saved points, so there is no per-step
+    interpolation data left to store. Use `save_everystep` or `saveat` to control
+    how finely that interpolant resolves the trajectory.
+
+## Callbacks and events
+
+The standard SciML callbacks work through the `callback` keyword:
+
+```julia
+# Stop as soon as the first component drops below a threshold.
+cb = ContinuousCallback((u, t, integrator) -> u[1] - 0.5, terminate!)
+sol = solve!(init(prob, alg; dt = 0.1, callback = cb))
+```
+
+The machinery is DiffEqBase's own, so `DiscreteCallback`, `ContinuousCallback`,
+`VectorContinuousCallback`, `CallbackSet` and anything built on top of them (for
+instance DiffEqCallbacks.jl) all behave as they do elsewhere in SciML.
+
+!!! note "Callbacks run on the outer integrator only"
+    A condition is evaluated once per **outer** step, after all the operators of
+    that step have been applied, and never between two inner splits: for the reason
+    given under "Saving and interpolation" above, those intermediate states are
+    stages and approximate the split solution at no time point, so there is nothing
+    meaningful for a condition to test or an `affect!` to modify at that level.
+
+    A consequence worth knowing: an `affect!` that modifies `integrator.u` is
+    propagated into every subintegrator before the next step, so modifying the state
+    from a callback is safe.
+
+### Accuracy of continuous events
+
+Event times are found by root-finding on the step's interpolant, which is linear
+(see "Saving and interpolation" above). Two consequences:
+
+  - The located event time is second order accurate in the step size, and is the
+    *exact* root of the linear interpolant over the step that brackets it. With
+    large steps -- adaptive splittings can grow the step considerably on smooth
+    problems -- the event time degrades accordingly. Cap it with `dtmax` when event
+    accuracy matters.
+  - An event that occurs and reverses **within a single step** cannot be detected,
+    because a linear interpolant has no interior extremum. Raising `interp_points`
+    does not help for the same reason, so setting `interp_points = 0` on the
+    callback avoids a sweep that cannot find anything the endpoints missed.
+
+Once the event time is located, the state there comes from the same interpolant and
+the whole subintegrator tree is re-anchored to it, so integration resumes
+consistently from the event.
+
 ## Configuring individual subintegrators
 
 `init` takes one value per keyword, which is not enough when the operators want
diff --git a/src/OrdinaryDiffEqOperatorSplitting.jl b/src/OrdinaryDiffEqOperatorSplitting.jl
index 1cb9578..19aaf77 100644
--- a/src/OrdinaryDiffEqOperatorSplitting.jl
+++ b/src/OrdinaryDiffEqOperatorSplitting.jl
@@ -42,6 +42,7 @@ abstract type AbstractOperatorSplittingAlgorithm end
 abstract type AbstractOperatorSplittingCache end
 
 @inline SciMLBase.isadaptive(::AbstractOperatorSplittingAlgorithm) = false
+@inline SciMLBase.isdiscrete(::AbstractOperatorSplittingAlgorithm) = false
 @inline isdtchangeable(alg::AbstractOperatorSplittingAlgorithm) = all(isdtchangeable.(alg.inner_algs))
 
 include("function.jl")
diff --git a/src/integrator.jl b/src/integrator.jl
index e6d5499..0d68db0 100644
--- a/src/integrator.jl
+++ b/src/integrator.jl
@@ -22,6 +22,23 @@ Base.@kwdef mutable struct IntegratorOptions{tType, fType, vbType, F3, atolType,
 end
 
 
+"""
+    SaveOptions
+
+The saving settings of the outermost [`OperatorSplittingIntegrator`](@ref).
+
+Consumed by `__init` directly rather than through the [`ConfigTree`](@ref), which is
+what keeps them from travelling down to the leaf integrators. `saveat` lives in the
+integrator's `saveat` heap instead, because it is consumed as the integration proceeds.
+"""
+struct SaveOptions
+    save_on::Bool
+    save_everystep::Bool
+    save_start::Bool
+    save_end::Bool
+end
+
+
 """
     SplitSubIntegratorStatus
 
@@ -103,7 +120,7 @@ mutable struct SplitSubIntegrator{
     controller_cache::controllerType
     force_stepfail::Bool
     last_step_failed::Bool
-    u_modified::Bool # TODO we can probably remove this
+    derivative_discontinuity::Bool
     status::SplitSubIntegratorStatus
     stats::IntegratorStats
     cache::cacheType
@@ -147,6 +164,8 @@ mutable struct OperatorSplittingIntegrator{
         tstopsType,
         saveatType,
         callbackType,
+        callbackCacheType,
+        eventErrType,
         cacheType,
         solType,
         subintTreeType,
@@ -172,14 +191,24 @@ mutable struct OperatorSplittingIntegrator{
     saveat::heapType
     _saveat::saveatType
     callback::callbackType
+    # Scratch buffers a `VectorContinuousCallback` needs; `nothing` when there is none.
+    const callback_cache::callbackCacheType
+    event_last_time::Int
+    vector_event_last_time::Int
+    last_event_error::eventErrType
     advance_to_tstop::Bool
     last_step_failed::Bool
     force_stepfail::Bool
     isout::Bool
-    u_modified::Bool
+    derivative_discontinuity::Bool
     just_hit_tstop::Bool
     cache::cacheType
     sol::solType
+    # `saveiter`, not `length(sol.t)`, is the authoritative length of the saved
+    # prefix of `sol.t`/`sol.u`.
+    const save_opts::SaveOptions
+    saveiter::Int
+    postamble_done::Bool
     # Tuple of SplitSubIntegrator nodes (one per top-level operator).
     child_subintegrators::subintTreeType
     child_solution_indices::childSolidxType # Tuple
@@ -229,7 +258,14 @@ function SciMLBase.__init(
         tstops = (),
         saveat = (),
         d_discontinuities = (),
+        save_on = true,
         save_everystep = false,
+        save_start = true,
+        save_end = true,
+        # Accepted and ignored: `sol(t)` already interpolates over the saved points
+        # with the same (linear) interpolant a dense output would store per step.
+        # Named explicitly so it cannot travel down to the leaf integrators.
+        dense = false,
         callback = nothing,
         advance_to_tstop = false,
         adaptive = nothing,
@@ -241,6 +277,8 @@ function SciMLBase.__init(
     (; u0, p) = prob
     t0, tf = prob.tspan
 
+    save_opts = SaveOptions(save_on, save_everystep, save_start, save_end)
+
     # By default every node adapts exactly if its own algorithm does; a scalar or a
     # TreeOption overrides the whole tree explicitly.
     if adaptive === nothing
@@ -280,9 +318,26 @@ function SciMLBase.__init(
     tmp = setup_u(prob, alg, false)
     uType = typeof(u)
 
-    sol = SciMLBase.build_solution(prob, alg, tType[], uType[])
+    # `build_solution` defaults `interp` to a `LinearInterpolation` over the very
+    # vectors passed in here, so `sol(t)` sees every later `savevalues!` push and
+    # agrees with the integrator's own (linear) interpolant. `stats` has to be a real
+    # `DEStats`: the callback machinery increments `sol.stats.ncondition`.
+    sol = SciMLBase.build_solution(
+        prob, alg, tType[], uType[];
+        stats = SciMLBase.DEStats(0),
+        calculate_error = false,
+    )
     callback = DiffEqBase.CallbackSet(callback)
 
+    # The scratch buffers have to be wide enough for the widest callback in the set.
+    max_len_cb = DiffEqBase.max_vector_callback_length_int(callback)
+    eventErrType = real(eltype(u))
+    callback_cache = if max_len_cb === nothing
+        nothing
+    else
+        DiffEqBase.CallbackCache(u, max_len_cb, eventErrType, eventErrType)
+    end
+
     child_subintegrators = build_subintegrators(
         prob, alg,
         uprev, u,
@@ -315,9 +370,12 @@ function SciMLBase.__init(
         tstops_internal, tstops,
         saveat_internal, saveat,
         callback,
+        callback_cache,
+        0, 1, zero(eventErrType),   # event_last_time, vector_event_last_time, last_event_error
         advance_to_tstop,
         false, false, false, false, false,
         cache, sol,
+        save_opts, 0, false,
         child_subintegrators,
         child_solution_indices,
         child_synchronizers,
@@ -330,10 +388,48 @@ function SciMLBase.__init(
         false,
         config,
     )
-    DiffEqBase.initialize!(callback, u0, t0, integrator)
+    # Initialize first, so that an initializer which modifies `u` is reflected in the
+    # first saved point.
+    initialize_callbacks!(integrator)
+    save_initial_value!(integrator)
     return integrator
 end
 
+# Run the callbacks' `initialize` hooks, mirroring OrdinaryDiffEqCore's
+# `initialize_callbacks!`. The flag starts out set so that a hook can clear it; the
+# default initializer does exactly that.
+function initialize_callbacks!(integrator::OperatorSplittingIntegrator)
+    integrator.derivative_discontinuity = true
+    # `integrator.u`, not `prob.u0`: with `alias_u0 = false` they are different arrays
+    # and a hook has to be able to modify the live state.
+    modified = DiffEqBase.initialize!(
+        integrator.callback, integrator.u, integrator.t, integrator
+    )
+    if modified
+        update_uprev!(integrator)
+        rollback_children!(integrator)
+    end
+    integrator.derivative_discontinuity = false
+    return nothing
+end
+
+function SciMLBase.reeval_internals_due_to_modification!(
+        integrator::OperatorSplittingIntegrator, continuous_modification = true;
+        callback_initializealg = nothing
+    )
+    rollback_children!(integrator)
+    integrator.derivative_discontinuity = false
+    return nothing
+end
+
+function save_initial_value!(integrator::OperatorSplittingIntegrator)
+    integrator.saveiter = 0
+    if integrator.save_opts.save_on && integrator.save_opts.save_start
+        _save_current!(integrator)
+    end
+    return nothing
+end
+
 # ---------------------------------------------------------------------------
 # reinit!
 # ---------------------------------------------------------------------------
@@ -379,12 +475,9 @@ function DiffEqBase.reinit!(
         resize!(integrator.sol.t, 0)
         resize!(integrator.sol.u, 0)
     end
-    if reinit_callbacks
-        DiffEqBase.initialize!(integrator.callback, u0, t0, integrator)
-    else
-        saving_callback = integrator.callback.discrete_callbacks[end]
-        DiffEqBase.initialize!(saving_callback, u0, t0, integrator)
-    end
+    integrator.event_last_time = 0
+    integrator.vector_event_last_time = 1
+    integrator.last_event_error = zero(integrator.last_event_error)
     if reinit_retcode
         integrator.sol = SciMLBase.solution_new_retcode(
             integrator.sol, ReturnCode.Default
@@ -400,6 +493,20 @@ function DiffEqBase.reinit!(
         erase_sol, tstops, saveat,
         reinit_callbacks, reinit_retcode
     )
+    # A leaf's `reinit!` restores it to the `u0` slice captured when it was built, not
+    # to the `u0` handed to this call, so the new state has to be pushed down
+    # explicitly. The forward sync of the next step is not enough: the palindromic
+    # schemes skip it for their first child.
+    rollback_children!(integrator)
+
+    # After the children, so the child reinit cannot undo an initializer's changes.
+    if reinit_callbacks
+        initialize_callbacks!(integrator)
+    end
+    # Redone regardless of `reinit_callbacks`, since saving is built into the
+    # integrator rather than provided by a saving callback.
+    integrator.postamble_done = false
+    save_initial_value!(integrator)
     return nothing
 end
 
@@ -577,10 +684,13 @@ end
 
 # Roll a node's children back to the state of the node itself: the local solution
 # buffers are refilled from the parent's (already restored) `u` and the child clocks
-# are moved back to the parent's time. Leaves get their `u` restored here as well,
-# not only by the forward sync before the next solve: palindromic algorithms skip
-# that sync for the first child (`next_sync_is_continuous`), so a rollback that left
-# the leaf state stale would silently resume from the failed attempt.
+# are moved back to the parent's time. Also the resync path for a callback that
+# modified `u`.
+#
+# Leaves get their `u` restored here as well, not only by the forward sync before the
+# next solve: palindromic algorithms skip that sync for their first child
+# (`next_sync_is_continuous`), so a rollback that left the leaf state stale would
+# silently resume from the failed attempt.
 rollback_children!(parent::AnySplitIntegrator) = _rollback_children!(
     parent.child_subintegrators, parent.child_solution_indices, parent.u, parent.t
 )
@@ -649,7 +759,7 @@ function step_header!(integrator::AnySplitIntegrator)
         else
             reject_step!(integrator)
         end
-    elseif integrator.u_modified
+    elseif integrator.derivative_discontinuity
         update_uprev!(integrator)
     end
     increment_iteration(integrator)
@@ -682,8 +792,10 @@ is_first_iteration(integrator::AnySplitIntegrator) = integrator.iter == 0
 increment_iteration(integrator::AnySplitIntegrator) = integrator.iter += 1
 
 function footer_reset_flags!(integrator)
-    integrator.u_modified = false
+    integrator.derivative_discontinuity = false
     integrator.just_hit_tstop = false
+    # Re-armed, so that stepping on after a `solve!` can close the solution out again.
+    integrator.postamble_done = false
     return
 end
 footer_reset_flags!(::SplitSubIntegrator) = nothing
@@ -692,14 +804,6 @@ function setup_validity_flags!(integrator, t_next)
     return
 end
 setup_validity_flags!(::SplitSubIntegrator, _) = nothing
-function fix_solution_buffer_sizes!(integrator, sol)
-    resize!(integrator.sol.t, integrator.saveiter)
-    resize!(integrator.sol.u, integrator.saveiter)
-    if !(integrator.sol isa SciMLBase.DAESolution)
-        resize!(integrator.sol.k, integrator.saveiter_dense)
-    end
-    return
-end
 
 # Window for absorbing floating point drift when landing on a time point. Scaled
 # by the local time scale *and* the step size: near t = 0 (e.g. integrating
@@ -753,6 +857,7 @@ function step_footer!(integrator::AnySplitIntegrator)
         try_snap_children_to_tstop!.(integrator.child_subintegrators, integrator.t)
         step_accept_controller!(integrator)
         validate_time_point(integrator)
+        handle_callbacks!(integrator)   # also does the saving
     elseif integrator.force_stepfail
         # Failure escalation protocol: the failing node's own adaptivity decides.
         fatal_rc = _fatal_child_retcode(integrator.child_subintegrators)
@@ -928,6 +1033,19 @@ function SciMLBase.check_error!(integrator::SplitSubIntegrator)
     return code
 end
 
+# Only a real failure finalizes here. The generic method finalizes on anything but
+# `Success`, but a mid-solve node legitimately reports `Default`, and this is called
+# before every step -- so the postamble would run once per step.
+function SciMLBase.check_error!(integrator::OperatorSplittingIntegrator)
+    code = SciMLBase.check_error(integrator)
+    # Rebuilding the solution allocates, and this runs before every step.
+    integrator.sol.retcode !== code && _set_retcode!(integrator, code)
+    if code ∉ (ReturnCode.Success, ReturnCode.Default)
+        SciMLBase.postamble!(integrator)
+    end
+    return code
+end
+
 @unroll function _check_error_children(current_retcode, children::Tuple)
     @unroll for child in children
         rc = _child_retcode(child)
@@ -949,15 +1067,244 @@ end
     return (integrator.tmp,)
 end
 
-function linear_interpolation!(y, t, y1, y2, t1, t2)
-    return y .= y1 + (t - t1) * (y2 - y1) / (t2 - t1)
+# ---------------------------------------------------------------------------
+# Interpolation
+#
+# See the "Dense output" section of the developer documentation for why the only
+# dense output well defined here is the one built from the endpoints the outer
+# integrator owns, and hence why the generic fallback is linear.
+# ---------------------------------------------------------------------------
+
+"""
+    splitting_interpolant(integrator, cache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}})
+    splitting_interpolant!(out, integrator, cache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}})
+
+Dense output of a single operator splitting step, evaluated in the step-local
+coordinate `Θ = (t - tprev) / dt`, where `y₀`/`y₁` are the state at `tprev`/`t` and
+`dt` is the step that was taken. `D` is the requested derivative order and `idxs`
+selects components (`nothing` for all of them).
+
+This is the extension point for higher order dense output: dispatch on the
+algorithm's cache type. The fallback for any
+[`AbstractOperatorSplittingCache`](@ref) is linear interpolation.
+"""
+splitting_interpolant, splitting_interpolant!
+
+splitting_interpolant(
+    integrator, ::AbstractOperatorSplittingCache, Θ, dt, y₀, y₁, idxs, ::Type{Val{D}}
+) where {D} = SciMLBase.linear_interpolant(Θ, dt, y₀, y₁, idxs, Val{D})
+
+function splitting_interpolant!(
+        out, integrator, ::AbstractOperatorSplittingCache, Θ, dt, y₀, y₁, idxs,
+        ::Type{Val{D}}
+    ) where {D}
+    SciMLBase.linear_interpolant!(out, Θ, dt, y₀, y₁, idxs, Val{D})
+    return out
+end
+
+# Step-local coordinate `Θ` and the length of the step that was just taken. The latter
+# must not come from `integrator.dt`: once a step is accepted `step_accept_controller!`
+# has overwritten it with the *next* proposal. A zero-length interval only happens
+# before the first step, where `uprev == u` and every Θ gives the same value.
+function _interp_coords(integrator::OperatorSplittingIntegrator, t)
+    dt = integrator.t - integrator.tprev
+    Θ = iszero(dt) ? zero(t / oneunit(dt)) : (t - integrator.tprev) / dt
+    return Θ, dt
+end
+
+function (integrator::OperatorSplittingIntegrator)(
+        t::Number, ::Type{deriv} = Val{0}; idxs = nothing
+    ) where {deriv}
+    Θ, dt = _interp_coords(integrator, t)
+    return splitting_interpolant(
+        integrator, integrator.cache, Θ, dt,
+        integrator.uprev, integrator.u, idxs, deriv
+    )
 end
-function (integrator::OperatorSplittingIntegrator)(tmp, t)
-    return linear_interpolation!(
-        tmp, t, integrator.uprev, integrator.u, integrator.tprev, integrator.t
+
+function (integrator::OperatorSplittingIntegrator)(
+        val::AbstractArray, t::Number, ::Type{deriv} = Val{0}; idxs = nothing
+    ) where {deriv}
+    Θ, dt = _interp_coords(integrator, t)
+    return splitting_interpolant!(
+        val, integrator, integrator.cache, Θ, dt,
+        integrator.uprev, integrator.u, idxs, deriv
     )
 end
 
+# ---------------------------------------------------------------------------
+# change_t_via_interpolation!
+# ---------------------------------------------------------------------------
+
+"""
+    change_t_via_interpolation!(integrator::OperatorSplittingIntegrator, t, modify_save_endpoint = Val{false}, reinitialize_alg = nothing)
+
+Move the integrator back to a time `t` inside the step that was just taken.
+
+`u` is refilled from the step's interpolant and the whole subintegrator tree is
+re-anchored to `t` through [`rollback_children!`](@ref), so every child's state and
+clock stay consistent with the parent (`validate_time_point` asserts the latter).
+`integrator.dt` is deliberately left alone: it already holds the next step's proposal.
+"""
+function SciMLBase.change_t_via_interpolation!(
+        integrator::OperatorSplittingIntegrator, t,
+        modify_save_endpoint::Type{Val{T}} = Val{false},
+        reinitialize_alg = nothing
+    ) where {T}
+    if integrator.tdir * t < integrator.tdir * integrator.tprev
+        error("Current interpolant only works between tprev and t")
+    elseif t != integrator.t
+        integrator(integrator.u, t)
+        integrator.t = t
+        rollback_children!(integrator)
+        # Compared by identity, not `if T`: DiffEqBase's continuous callback path
+        # passes `Val{:false}`, a `Val` of the *Symbol*, which `if T` would throw on.
+        if modify_save_endpoint === Val{true}
+            save_endpoint!(integrator)
+        end
+    end
+    return nothing
+end
+
+# ---------------------------------------------------------------------------
+# Saving
+#
+# Outer integrator only, for the reason given in the "Dense output" section of the
+# developer documentation.
+# ---------------------------------------------------------------------------
+
+_saved_at_current_t(integrator::OperatorSplittingIntegrator) =
+    integrator.saveiter > 0 &&
+    integrator.sol.t[integrator.saveiter] == integrator.t
+
+# `copy_u = false` hands over ownership of a `u` the caller freshly allocated.
+function _save_at!(integrator::OperatorSplittingIntegrator, t, u, copy_u = true)
+    integrator.saveiter += 1
+    RecursiveArrayTools.copyat_or_push!(integrator.sol.t, integrator.saveiter, t)
+    RecursiveArrayTools.copyat_or_push!(integrator.sol.u, integrator.saveiter, u, copy_u)
+    return nothing
+end
+
+_save_current!(integrator::OperatorSplittingIntegrator) =
+    _save_at!(integrator, integrator.t, integrator.u)
+
+"""
+    savevalues!(integrator::OperatorSplittingIntegrator, force_save = false, reduce_size = true)
+
+Save the state at the current time point and at any pending `saveat` points that the
+step just taken has passed, returning `(saved, savedexactly)`.
+
+`saveat` points strictly inside the step are filled from the step's interpolant
+(see [`splitting_interpolant`](@ref)), so asking for output never changes the sequence
+of steps and therefore never changes the splitting error. The generic interpolant is
+linear, so such points are only first order accurate; pass the times as `tstops`
+instead to have the integrator land on them exactly.
+"""
+function SciMLBase.savevalues!(
+        integrator::OperatorSplittingIntegrator,
+        force_save = false, reduce_size = true
+    )::Tuple{Bool, Bool}
+    saved = false
+    savedexactly = false
+    integrator.save_opts.save_on || return saved, savedexactly
+    (; save_everystep, save_end) = integrator.save_opts
+    (force_save || save_everystep || !isempty(integrator.saveat)) ||
+        return saved, savedexactly
+    tf = integrator.sol.prob.tspan[2]
+
+    # The heaps store raw times and carry the direction in their ordering, so the
+    # comparison is scaled by tdir rather than the stored value.
+    tdir_t = integrator.tdir * integrator.t
+    while !isempty(integrator.saveat) &&
+            integrator.tdir * first(integrator.saveat) <= tdir_t
+        curt = pop!(integrator.saveat)
+        if curt == integrator.t
+            # `save_end` owns the final point; leave it to the postamble.
+            (!save_end && curt == tf) && continue
+            savedexactly = true
+            _save_current!(integrator)
+        else
+            _save_at!(integrator, curt, integrator(curt), false)
+        end
+        saved = true
+    end
+
+    if force_save || (
+            save_everystep && !_saved_at_current_t(integrator) &&
+                (save_end || integrator.t != tf)
+        )
+        saved = true
+        savedexactly = true
+        _save_current!(integrator)
+    end
+
+    return saved, savedexactly
+end
+
+function save_endpoint!(integrator::OperatorSplittingIntegrator)
+    integrator.save_opts.save_on || return nothing
+    integrator.save_opts.save_end || return nothing
+    _saved_at_current_t(integrator) && return nothing
+    _save_current!(integrator)
+    return nothing
+end
+
+# ---------------------------------------------------------------------------
+# Callbacks
+#
+# Outer integrator only, for the reason given in the "Dense output" section of the
+# developer documentation. Root finding, `save_positions` handling and `affect!`
+# dispatch are DiffEqBase's; this is just the per-step driver, modelled on
+# OrdinaryDiffEqCore's `handle_callbacks!`, and it also does the saving when no
+# callback saved.
+# ---------------------------------------------------------------------------
+
+handle_callbacks!(::SplitSubIntegrator) = nothing
+
+function handle_callbacks!(integrator::OperatorSplittingIntegrator)
+    discrete_callbacks = integrator.callback.discrete_callbacks
+    continuous_callbacks = integrator.callback.continuous_callbacks
+
+    continuous_modified = false
+    discrete_modified = false
+    saved_in_cb = false
+
+    if !(continuous_callbacks isa Tuple{})
+        time, upcrossing, event_occurred, event_idx, idx, _counter =
+            DiffEqBase.find_first_continuous_callback(
+            integrator, continuous_callbacks...
+        )
+        if event_occurred
+            integrator.event_last_time = idx
+            integrator.vector_event_last_time = event_idx
+            continuous_modified, saved_in_cb = OrdinaryDiffEqCore.apply_ith_callback!(
+                integrator, time, upcrossing, event_idx, idx, continuous_callbacks
+            )
+            # The step was cut at the event and the state may have jumped, so the
+            # controller's error history no longer describes what happens next.
+            reinit_node_controller!(integrator)
+        else
+            # Clearing these stops the root finder from nudging `tprev` and
+            # re-detecting the event handled on an earlier step.
+            integrator.event_last_time = 0
+            integrator.vector_event_last_time = 1
+        end
+    end
+
+    if !integrator.force_stepfail && !(discrete_callbacks isa Tuple{})
+        discrete_modified, saved_in_cb = DiffEqBase.apply_discrete_callback!(
+            integrator, discrete_callbacks...
+        )
+    end
+
+    if !saved_in_cb
+        SciMLBase.savevalues!(integrator)
+    end
+
+    integrator.derivative_discontinuity = continuous_modified | discrete_modified
+    return nothing
+end
+
 # ---------------------------------------------------------------------------
 # Step size control
 #
@@ -1116,7 +1463,16 @@ function SciMLBase.done(integrator::OperatorSplittingIntegrator)
 end
 
 function SciMLBase.postamble!(integrator::OperatorSplittingIntegrator)
-    return DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator)
+    # `solve!`, `done` and `check_error!` all call this; only the first one after a
+    # step may run the finalizers and close out the solution.
+    integrator.postamble_done && return nothing
+    integrator.postamble_done = true
+    DiffEqBase.finalize!(integrator.callback, integrator.u, integrator.t, integrator)
+    save_endpoint!(integrator)
+    # Drop whatever a previous, longer run left behind (`reinit!` without `erase_sol`).
+    resize!(integrator.sol.t, integrator.saveiter)
+    resize!(integrator.sol.u, integrator.saveiter)
+    return nothing
 end
 
 function __step!(integrator::AnySplitIntegrator)
@@ -1325,7 +1681,7 @@ function _build_child(
         0, 0,           # iter, success_iter
         EEst_val,
         controller_cache,
-        false, false, false,  # force_stepfail, last_step_failed, u_modified
+        false, false, false,  # force_stepfail, last_step_failed, derivative_discontinuity
         SplitSubIntegratorStatus(),
         IntegratorStats(),
         level_cache,
@@ -1412,6 +1768,31 @@ SciMLBase.first_tstop(i::AnySplitIntegrator) = first(i.tstops)
 SciMLBase.pop_tstop!(i::AnySplitIntegrator) = pop!(i.tstops)
 
 DiffEqBase.get_dt(i::AnySplitIntegrator) = i.dt
+
+function SciMLBase.set_proposed_dt!(integrator::OperatorSplittingIntegrator, dt)
+    if integrator.dtcache != abs(dt)
+        integrator.dtcache = abs(dt)
+        if !isadaptive(integrator)
+            integrator.dt = dt
+        end
+    end
+    return nothing
+end
+
+"""
+    terminate!(integrator::OperatorSplittingIntegrator, retcode = ReturnCode.Terminated)
+
+Stop the integration from inside a callback's `affect!`. Emptying the outermost
+`tstops` heap is what ends the `solve!` loop and makes `done` report completion.
+"""
+function SciMLBase.terminate!(
+        integrator::OperatorSplittingIntegrator, retcode = ReturnCode.Terminated
+    )
+    _set_retcode!(integrator, retcode)
+    empty!(integrator.tstops)
+    return nothing
+end
+
 function set_dt!(i::DEIntegrator, dt)
     iszero(dt) && error("dt must be nonzero")
     return i.dt = dt
@@ -1434,12 +1815,5 @@ function DiffEqBase.add_saveat!(i::OperatorSplittingIntegrator, t)
     return nothing
 end
 
-# SciMLBase v3 renamed `u_modified!` → `derivative_discontinuity!`.
-@static if isdefined(DiffEqBase, :u_modified!)
-    DiffEqBase.u_modified!(i::OperatorSplittingIntegrator, bool) = i.u_modified = bool
-    DiffEqBase.u_modified!(i::SplitSubIntegrator, bool) = i.u_modified = bool
-end
-@static if isdefined(SciMLBase, :derivative_discontinuity!)
-    SciMLBase.derivative_discontinuity!(i::OperatorSplittingIntegrator, bool) = i.u_modified = bool
-    SciMLBase.derivative_discontinuity!(i::SplitSubIntegrator, bool) = i.u_modified = bool
-end
+SciMLBase.derivative_discontinuity!(i::AnySplitIntegrator, bool) =
+    i.derivative_discontinuity = bool
diff --git a/src/utils.jl b/src/utils.jl
index 642fd79..ab3e99d 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -7,16 +7,17 @@ function tstops_and_saveat_heaps(t0, tf, tstops, saveat)
     tstops = [filter(t -> t0 < t < tf || tf < t < t0, tstops)..., tf]
     tstops = BinaryHeaps.BinaryHeap{FT, ordering}(tstops)
 
-    if isnothing(saveat)
-        saveat = [t0, tf]
+    # Keep `t0 < t <= tf` in tdir-space: `save_start` owns the initial point and
+    # `save_end` the final one, so leaving either in the heap would duplicate it.
+    tdir = tf > t0 ? one(FT) : -one(FT)
+    saveat = if isnothing(saveat)
+        FT[]
     elseif saveat isa Number
         saveat > zero(saveat) || error("saveat value must be positive")
-        saveat = tf > t0 ? saveat : -saveat
-        saveat = [t0:saveat:tf..., tf]
+        step = tdir * saveat
+        collect((t0 + step):step:tf)
     else
-        # We do not need to filter saveat like tstops because the saving
-        # callback will ignore any times that are not between t0 and tf.
-        saveat = collect(saveat)
+        filter(t -> tdir * t0 < tdir * t <= tdir * tf, collect(FT, saveat))
     end
     saveat = BinaryHeaps.BinaryHeap{FT, ordering}(saveat)
 
@@ -75,14 +76,9 @@ function forward_sync_subintegrator!(
 end
 
 # Tell a leaf integrator that its state was changed from the outside so it discards
-# FSAL information. SciMLBase v3 renamed `u_modified!` → `derivative_discontinuity!`;
-# call the appropriate name based on which SciMLBase is loaded.
+# FSAL information.
 function mark_state_modified!(child::DEIntegrator)
-    @static if isdefined(SciMLBase, :derivative_discontinuity!)
-        SciMLBase.derivative_discontinuity!(child, true)
-    else
-        SciMLBase.u_modified!(child, true)
-    end
+    SciMLBase.derivative_discontinuity!(child, true)
     return nothing
 end
 
diff --git a/test/callbacks.jl b/test/callbacks.jl
new file mode 100644
index 0000000..3ef3fe9
--- /dev/null
+++ b/test/callbacks.jl
@@ -0,0 +1,497 @@
+using OrdinaryDiffEqOperatorSplitting
+using Test
+
+import SciMLBase
+import SciMLBase: ReturnCode
+import DiffEqBase: DiffEqBase, ODEFunction, DiscreteCallback, ContinuousCallback,
+    VectorContinuousCallback, CallbackSet
+using OrdinaryDiffEqLowOrderRK
+using OrdinaryDiffEqTsit5
+
+const OS = OrdinaryDiffEqOperatorSplitting
+
+# ---------------------------------------------------------------------------
+# Reference problem. The third component is decoupled and solved by the first
+# operator alone, so `u[3](t) = 3exp(-t/10)` exactly -- which gives continuous
+# callback tests an analytic event time to compare against.
+# ---------------------------------------------------------------------------
+const U0 = [1.0, 2.0, 3.0]
+const TSPAN = (0.0, 1.0)
+
+f1 = ODEFunction((du, u, p, t) -> (@. du = -0.1u))
+f2 = ODEFunction(
+    function (du, u, p, t)
+        du[1] = -0.01u[2]
+        du[2] = -0.01u[1]
+        return nothing
+    end
+)
+fsplit = GenericSplitFunction((f1, f2), ([1, 2, 3], [1, 2]))
+
+make_prob(tspan = TSPAN) = OperatorSplittingProblem(fsplit, copy(U0), tspan)
+ltg() = LieTrotterGodunov((Euler(), Euler()))
+# Accurate inner solvers, so that an event time's error is the interpolant's alone.
+ltg_exact() = LieTrotterGodunov((Tsit5(), Tsit5()))
+
+# Exact crossing time of u[3] = level.
+exact_crossing(level) = -10 * log(level / 3)
+
+@testset "discrete callbacks" begin
+    prob = make_prob()
+
+    @testset "condition is evaluated once per outer step" begin
+        calls = Ref(0)
+        cb = DiscreteCallback((u, t, integrator) -> (calls[] += 1; false), integrator -> nothing)
+        DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb))
+        @test calls[] == 10
+    end
+
+    @testset "nested splitting does not run callbacks per inner stage" begin
+        inner = GenericSplitFunction((f1, f1), ([1, 2], [1, 2]))
+        outer = GenericSplitFunction((f1, inner), ([1, 2, 3], [1, 2]))
+        nprob = OperatorSplittingProblem(outer, copy(U0), TSPAN)
+        alg = LieTrotterGodunov((Euler(), LieTrotterGodunov((Euler(), Euler()))))
+        calls = Ref(0)
+        cb = DiscreteCallback((u, t, integrator) -> (calls[] += 1; false), integrator -> nothing)
+        DiffEqBase.solve!(DiffEqBase.init(nprob, alg; dt = 0.1, callback = cb))
+        @test calls[] == 10   # outer steps only, whatever the inner node does
+    end
+
+    @testset "affect! fires and sees the integrator" begin
+        seen = Float64[]
+        cb = DiscreteCallback(
+            (u, t, integrator) -> t >= 0.45,
+            integrator -> push!(seen, integrator.t)
+        )
+        sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb))
+        @test seen ≈ [0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
+        @test sol.retcode == ReturnCode.Success
+    end
+
+    @testset "save_positions brackets the affect!" begin
+        cb = DiscreteCallback(
+            (u, t, integrator) -> isapprox(t, 0.5),
+            integrator -> (integrator.u[1] += 10.0)
+        )
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb, tstops = [0.5])
+        )
+        # Two points at the event time: the state before and after the affect!.
+        idxs = findall(t -> isapprox(t, 0.5), sol.t)
+        @test length(idxs) == 2
+        @test sol.u[idxs[2]][1] - sol.u[idxs[1]][1] ≈ 10.0
+
+        # save_positions = (false, false) records nothing extra.
+        cb_quiet = DiscreteCallback(
+            (u, t, integrator) -> isapprox(t, 0.5),
+            integrator -> (integrator.u[1] += 10.0);
+            save_positions = (false, false)
+        )
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb_quiet, tstops = [0.5])
+        )
+        @test count(t -> isapprox(t, 0.5), sol.t) == 0
+    end
+
+    @testset "a modified u reaches the whole subintegrator tree" begin
+        # Every child holds its own copy of its slice, and the palindromic schemes
+        # skip the forward sync of their first child.
+        for alg in (
+                ltg(),
+                StrangMarchuk((Euler(), Euler())),
+                PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())),
+            )
+            cb = DiscreteCallback(
+                (u, t, integrator) -> isapprox(t, 0.5),
+                integrator -> (integrator.u .= [100.0, 200.0, 300.0])
+            )
+            integrator = DiffEqBase.init(
+                prob, alg; dt = 0.1, callback = cb, tstops = [0.5]
+            )
+            while integrator.t < 0.5
+                DiffEqBase.step!(integrator)
+            end
+            @test integrator.u == [100.0, 200.0, 300.0]
+            @test integrator.child_subintegrators[1].u == [100.0, 200.0, 300.0]
+            @test integrator.child_subintegrators[2].u == [100.0, 200.0]
+            OS.validate_time_point(integrator)
+        end
+    end
+
+    @testset "a modified u actually changes the trajectory" begin
+        plain = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, tstops = [0.5])
+        )
+        cb = DiscreteCallback(
+            (u, t, integrator) -> isapprox(t, 0.5),
+            integrator -> (integrator.u .*= 2)
+        )
+        bumped = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb, tstops = [0.5])
+        )
+        # Doubling halfway through must roughly double the endpoint of a linear ODE.
+        @test bumped.u[end] ≈ 2 .* plain.u[end] rtol = 1.0e-8
+    end
+
+    @testset "multiple callbacks in a CallbackSet" begin
+        # Thresholds sit between grid points, and times are compared with isapprox,
+        # because of ulp drift (0.1 * 8 is 0.7999999999999999, not 0.8).
+        early = Float64[]
+        late = Float64[]
+        cb1 = DiscreteCallback((u, t, integrator) -> t >= 0.35, integrator -> push!(early, integrator.t))
+        cb2 = DiscreteCallback((u, t, integrator) -> t >= 0.85, integrator -> push!(late, integrator.t))
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, callback = CallbackSet(cb1, cb2))
+        )
+        @test early ≈ [0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
+        @test late ≈ [0.9, 1.0]
+        @test sol.retcode == ReturnCode.Success
+    end
+
+    @testset "sol.stats.ncondition counts condition evaluations" begin
+        cb = DiscreteCallback((u, t, integrator) -> false, integrator -> nothing)
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)
+        sol = DiffEqBase.solve!(integrator)
+        @test sol.stats.ncondition == 10
+    end
+end
+
+@testset "terminate!" begin
+    prob = make_prob()
+
+    @testset "stops the solve and sets the retcode" begin
+        cb = DiscreteCallback(
+            (u, t, integrator) -> t >= 0.5,
+            integrator -> SciMLBase.terminate!(integrator)
+        )
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)
+        sol = DiffEqBase.solve!(integrator)
+        @test sol.retcode == ReturnCode.Terminated
+        @test sol.t[end] ≈ 0.5
+        @test integrator.t ≈ 0.5
+        @test SciMLBase.done(integrator)
+    end
+
+    @testset "accepts an explicit retcode" begin
+        cb = DiscreteCallback(
+            (u, t, integrator) -> t >= 0.5,
+            integrator -> SciMLBase.terminate!(integrator, ReturnCode.Success)
+        )
+        sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb))
+        @test sol.retcode == ReturnCode.Success
+    end
+end
+
+@testset "continuous callbacks" begin
+    prob = make_prob()
+    level = 2.9
+    exact = exact_crossing(level)
+
+    @testset "locates the crossing and fires exactly once" begin
+        hits = Float64[]
+        cb = ContinuousCallback(
+            (u, t, integrator) -> u[3] - level,
+            integrator -> push!(hits, integrator.t)
+        )
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg_exact(); dt = 0.05, callback = cb)
+        )
+        @test length(hits) == 1
+        @test hits[1] ≈ exact atol = 1.0e-3
+        @test hits[1] in sol.t
+        @test issorted(sol.t)
+        @test sol.retcode == ReturnCode.Success
+    end
+
+    @testset "event time converges under refinement" begin
+        # The located root is at least second order accurate. The step grid moves
+        # with dt, so the constant is noisy: assert only that a factor-of-four
+        # refinement buys a factor of four.
+        errs = map((0.1, 0.025)) do dt
+            hits = Float64[]
+            cb = ContinuousCallback(
+                (u, t, integrator) -> u[3] - level,
+                integrator -> push!(hits, integrator.t)
+            )
+            DiffEqBase.solve!(
+                DiffEqBase.init(prob, ltg_exact(); dt = dt, callback = cb)
+            )
+            abs(hits[1] - exact)
+        end
+        @test errs[2] < errs[1] / 4
+    end
+
+    @testset "affect_neg! handles the downcrossing" begin
+        ups = Ref(0)
+        downs = Ref(0)
+        cb = ContinuousCallback(
+            (u, t, integrator) -> u[3] - level,
+            integrator -> (ups[] += 1);
+            affect_neg! = integrator -> (downs[] += 1)
+        )
+        DiffEqBase.solve!(DiffEqBase.init(prob, ltg_exact(); dt = 0.05, callback = cb))
+        # u[3] decays through the level, so this is a downcrossing only.
+        @test ups[] == 0
+        @test downs[] == 1
+    end
+
+    @testset "the tree is re-anchored to the event time" begin
+        cb = ContinuousCallback(
+            (u, t, integrator) -> u[3] - level,
+            integrator -> nothing
+        )
+        integrator = DiffEqBase.init(prob, ltg_exact(); dt = 0.05, callback = cb)
+        while integrator.t < exact
+            DiffEqBase.step!(integrator)
+        end
+        @test integrator.t ≈ exact atol = 1.0e-3
+        OS.validate_time_point(integrator)
+        @test integrator.child_subintegrators[1].u ≈ integrator.u
+        @test integrator.child_subintegrators[2].u ≈ integrator.u[1:2]
+    end
+
+    @testset "an affect! that modifies u reaches the children" begin
+        cb = ContinuousCallback(
+            (u, t, integrator) -> u[3] - level,
+            integrator -> (integrator.u .= [7.0, 8.0, 9.0])
+        )
+        integrator = DiffEqBase.init(prob, ltg_exact(); dt = 0.05, callback = cb)
+        while integrator.t < exact
+            DiffEqBase.step!(integrator)
+        end
+        @test integrator.u == [7.0, 8.0, 9.0]
+        @test integrator.child_subintegrators[1].u == [7.0, 8.0, 9.0]
+        @test integrator.child_subintegrators[2].u == [7.0, 8.0]
+    end
+
+    @testset "NoRootFind fires at the step endpoint" begin
+        hits = Float64[]
+        cb = ContinuousCallback(
+            (u, t, integrator) -> u[3] - level,
+            integrator -> push!(hits, integrator.t);
+            rootfind = SciMLBase.NoRootFind
+        )
+        DiffEqBase.solve!(DiffEqBase.init(prob, ltg_exact(); dt = 0.1, callback = cb))
+        @test length(hits) == 1
+        # No root finding: the event is reported at the end of the bracketing step.
+        @test hits[1] ≈ 0.4
+    end
+
+    @testset "works with an adaptive splitting" begin
+        hits = Float64[]
+        cb = ContinuousCallback(
+            (u, t, integrator) -> u[3] - level,
+            integrator -> push!(hits, integrator.t)
+        )
+        # `dtmax` matters: this problem is so nearly linear that the controller
+        # otherwise spans the rest of the tspan in one step, and the event is only
+        # as accurate as the (linear) interpolant over the bracketing step.
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(
+                prob, PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5()));
+                dt = 0.05, dtmax = 0.05, callback = cb
+            )
+        )
+        @test length(hits) == 1
+        @test hits[1] ≈ exact atol = 1.0e-3
+        @test sol.retcode == ReturnCode.Success
+    end
+
+    @testset "a long adaptive step degrades the event time" begin
+        # An event is only as accurate as the interpolant over the step that brackets
+        # it, so capping the step with dtmax improves it.
+        function locate(; kwargs...)
+            hits = Float64[]
+            cb = ContinuousCallback(
+                (u, t, integrator) -> u[3] - level,
+                integrator -> push!(hits, integrator.t)
+            )
+            DiffEqBase.solve!(
+                DiffEqBase.init(
+                    prob, PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5()));
+                    dt = 0.05, callback = cb, kwargs...
+                )
+            )
+            return only(hits)
+        end
+
+        @test abs(locate(dtmax = 0.05) - exact) < abs(locate() - exact)
+    end
+
+    @testset "interp_points = 0 skips the safety sweep" begin
+        calls = Ref(0)
+        cb = ContinuousCallback(
+            (u, t, integrator) -> (calls[] += 1; u[3] - level),
+            integrator -> nothing;
+            interp_points = 0
+        )
+        DiffEqBase.solve!(DiffEqBase.init(prob, ltg_exact(); dt = 0.25, callback = cb))
+        sweeping = Ref(0)
+        cb2 = ContinuousCallback(
+            (u, t, integrator) -> (sweeping[] += 1; u[3] - level),
+            integrator -> nothing
+        )
+        DiffEqBase.solve!(DiffEqBase.init(prob, ltg_exact(); dt = 0.25, callback = cb2))
+        # The default `interp_points = 11` sweeps the interpolant for sign changes
+        # the endpoints missed; with `0` the condition is only seen at the endpoints.
+        @test calls[] < sweeping[]
+    end
+end
+
+@testset "VectorContinuousCallback" begin
+    prob = make_prob()
+    levels = (2.9, 2.85)
+
+    # Current DiffEqBase hands `affect!` a mask of simultaneous events rather than a
+    # single index; accept either shape.
+    fired_index(idx::Integer) = Int(idx)
+    fired_index(mask) = only(findall(!iszero, mask))
+
+    fired = Tuple{Int, Float64}[]
+    cb = VectorContinuousCallback(
+        function (out, u, t, integrator)
+            out[1] = u[3] - levels[1]
+            out[2] = u[3] - levels[2]
+            return nothing
+        end,
+        (integrator, idx) -> push!(fired, (fired_index(idx), integrator.t)),
+        2
+    )
+    integrator = DiffEqBase.init(prob, ltg_exact(); dt = 0.05, callback = cb)
+    @test integrator.callback_cache !== nothing
+    sol = DiffEqBase.solve!(integrator)
+
+    @test length(fired) == 2
+    @test first.(fired) == [1, 2]                       # the higher level is crossed first
+    @test last(fired[1]) ≈ exact_crossing(levels[1]) atol = 1.0e-3
+    @test last(fired[2]) ≈ exact_crossing(levels[2]) atol = 1.0e-3
+    @test sol.retcode == ReturnCode.Success
+
+    @testset "no cache is allocated without a vector callback" begin
+        scalar_cb = ContinuousCallback((u, t, integrator) -> u[3] - 2.9, integrator -> nothing)
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = scalar_cb)
+        @test integrator.callback_cache === nothing
+    end
+end
+
+@testset "callback initialization" begin
+    prob = make_prob()
+
+    @testset "an initializer can add tstops" begin
+        # The pattern PresetTimeCallback uses.
+        hits = Float64[]
+        cb = DiscreteCallback(
+            (u, t, integrator) -> t in (0.35, 0.65),
+            integrator -> push!(hits, integrator.t);
+            initialize = function (c, u, t, integrator)
+                DiffEqBase.add_tstop!(integrator, 0.35)
+                DiffEqBase.add_tstop!(integrator, 0.65)
+                return SciMLBase.derivative_discontinuity!(integrator, false)
+            end
+        )
+        sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb))
+        @test hits == [0.35, 0.65]
+        @test sol.retcode == ReturnCode.Success
+    end
+
+    @testset "an initializer that modifies u is reflected everywhere" begin
+        cb = DiscreteCallback(
+            (u, t, integrator) -> false,
+            integrator -> nothing;
+            initialize = function (c, u, t, integrator)
+                integrator.u .= [5.0, 6.0, 7.0]
+                return SciMLBase.derivative_discontinuity!(integrator, true)
+            end
+        )
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)
+        @test integrator.u == [5.0, 6.0, 7.0]
+        @test integrator.uprev == [5.0, 6.0, 7.0]
+        @test integrator.child_subintegrators[1].u == [5.0, 6.0, 7.0]
+        @test integrator.child_subintegrators[2].u == [5.0, 6.0]
+        # The first saved point is the modified state, not the problem's u0.
+        @test integrator.sol.u[1] == [5.0, 6.0, 7.0]
+    end
+
+    @testset "finalize runs exactly once" begin
+        finals = Ref(0)
+        cb = DiscreteCallback(
+            (u, t, integrator) -> false,
+            integrator -> nothing;
+            finalize = (c, u, t, integrator) -> (finals[] += 1)
+        )
+        DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb))
+        @test finals[] == 1
+    end
+
+    @testset "reinit! keeps an initializer's modification" begin
+        # The child reinit must not undo what the initializer did, so the callbacks
+        # have to be initialized after the subintegrator tree is restored.
+        cb = DiscreteCallback(
+            (u, t, integrator) -> false,
+            integrator -> nothing;
+            initialize = function (c, u, t, integrator)
+                integrator.u .= [5.0, 6.0, 7.0]
+                return SciMLBase.derivative_discontinuity!(integrator, true)
+            end
+        )
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)
+        DiffEqBase.solve!(integrator)
+        DiffEqBase.reinit!(integrator)
+        @test integrator.u == [5.0, 6.0, 7.0]
+        @test integrator.uprev == [5.0, 6.0, 7.0]
+        @test integrator.child_subintegrators[1].u == [5.0, 6.0, 7.0]
+        @test integrator.child_subintegrators[2].u == [5.0, 6.0]
+        @test integrator.sol.u[1] == [5.0, 6.0, 7.0]
+    end
+
+    @testset "reinit! re-runs the initializers" begin
+        inits = Ref(0)
+        cb = DiscreteCallback(
+            (u, t, integrator) -> false,
+            integrator -> nothing;
+            initialize = function (c, u, t, integrator)
+                inits[] += 1
+                return SciMLBase.derivative_discontinuity!(integrator, false)
+            end
+        )
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)
+        @test inits[] == 1
+        DiffEqBase.solve!(integrator)
+        DiffEqBase.reinit!(integrator)
+        @test inits[] == 2
+        DiffEqBase.reinit!(integrator; reinit_callbacks = false)
+        @test inits[] == 2
+        @test DiffEqBase.solve!(integrator).retcode == ReturnCode.Success
+    end
+end
+
+@testset "callbacks combined with saving" begin
+    prob = make_prob()
+
+    @testset "saveat and an event coexist" begin
+        cb = ContinuousCallback(
+            (u, t, integrator) -> u[3] - 2.9,
+            integrator -> nothing
+        )
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(
+                prob, ltg_exact(); dt = 0.1, saveat = [0.2, 0.8], callback = cb
+            )
+        )
+        @test issorted(sol.t)
+        @test 0.2 in sol.t
+        @test 0.8 in sol.t
+        @test sol.t[1] == 0.0
+        @test sol.t[end] == 1.0
+    end
+
+    @testset "a discrete callback still leaves a usable interpolation" begin
+        cb = DiscreteCallback((u, t, integrator) -> false, integrator -> nothing)
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(
+                prob, ltg(); dt = 0.1, callback = cb, save_everystep = true
+            )
+        )
+        @test sol(0.35) ≈ (sol.u[4] .+ sol.u[5]) ./ 2
+    end
+end
diff --git a/test/qa/qa.jl b/test/qa/qa.jl
index cd55cef..42f9f80 100644
--- a/test/qa/qa.jl
+++ b/test/qa/qa.jl
@@ -31,7 +31,12 @@ run_qa(
                 :setup_controller_cache, :reinit_controller!, :post_newton_controller!,  # OrdinaryDiffEqCore
                 :get_EEst, :set_EEst!, :get_current_adaptive_order,                      # OrdinaryDiffEqCore
                 :gamma_default, :failfactor_default, :AbstractControllerCache,           # OrdinaryDiffEqCore
+                :apply_ith_callback!,                                                    # OrdinaryDiffEqCore
                 :promote_tspan,                                                          # SciMLBase
+                # Shared linear interpolation kernels, reused so that the
+                # integrator's own dense output and `sol(t)` (which goes through
+                # SciMLBase's `LinearInterpolation`) cannot drift apart.
+                :linear_interpolant, :linear_interpolant!,                                # SciMLBase
                 # Broadcast extension points a TreeOption implements (src/config_tree.jl).
                 :Broadcasted, :broadcastable, :dotview, :materialize!,                   # Base
             ),
diff --git a/test/saving.jl b/test/saving.jl
new file mode 100644
index 0000000..bc545f3
--- /dev/null
+++ b/test/saving.jl
@@ -0,0 +1,404 @@
+using OrdinaryDiffEqOperatorSplitting
+using Test
+
+import SciMLBase
+import SciMLBase: ReturnCode
+import DiffEqBase: DiffEqBase, ODEFunction, DiscreteCallback, CallbackSet
+using OrdinaryDiffEqLowOrderRK
+using OrdinaryDiffEqTsit5
+
+const OS = OrdinaryDiffEqOperatorSplitting
+
+# ---------------------------------------------------------------------------
+# Reference problem: two linear operators sharing the first two components.
+# ---------------------------------------------------------------------------
+const U0 = [1.0, 2.0, 3.0]
+const TSPAN = (0.0, 1.0)
+
+f1 = ODEFunction((du, u, p, t) -> (@. du = -0.1u))
+f2 = ODEFunction(
+    function (du, u, p, t)
+        du[1] = -0.01u[2]
+        du[2] = -0.01u[1]
+        return nothing
+    end
+)
+fsplit = GenericSplitFunction((f1, f2), ([1, 2, 3], [1, 2]))
+
+make_prob(tspan = TSPAN) = OperatorSplittingProblem(fsplit, copy(U0), tspan)
+ltg() = LieTrotterGodunov((Euler(), Euler()))
+
+# Swept by every test that is not about one particular scheme.
+algs() = (
+    LieTrotterGodunov((Tsit5(), Tsit5())),
+    StrangMarchuk((Tsit5(), Tsit5())),
+    PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())),
+)
+
+function exact(t)
+    p = 1.5 * exp(-0.11t)     # (1, 1) mode
+    m = -0.5 * exp(-0.09t)    # (1, -1) mode
+    return [p + m, p - m, 3 * exp(-0.1t)]
+end
+rhs(u) = [-0.1u[1] - 0.01u[2], -0.1u[2] - 0.01u[1], -0.1u[3]]
+
+@testset "saved time points" begin
+    prob = make_prob()
+
+    @testset "default saves only the interval endpoints" begin
+        sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1))
+        @test sol.t == [0.0, 1.0]
+        @test length(sol.u) == 2
+        @test sol.u[1] == U0
+        @test sol.retcode == ReturnCode.Success
+    end
+
+    @testset "save_everystep saves every accepted step" begin
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, save_everystep = true)
+        )
+        @test length(sol.t) == 11          # t0 plus ten steps of dt = 0.1
+        @test sol.t[1] == 0.0
+        @test sol.t[end] == 1.0
+        @test issorted(sol.t)
+        # No duplicate of the final point from `save_end`.
+        @test allunique(sol.t)
+    end
+
+    @testset "saveat vector" begin
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = [0.25, 0.5])
+        )
+        @test sol.t == [0.0, 0.25, 0.5, 1.0]
+    end
+
+    @testset "saveat number excludes t0 and does not duplicate tf" begin
+        sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = 0.25))
+        @test sol.t == [0.0, 0.25, 0.5, 0.75, 1.0]
+        # A step that does not divide the span evenly leaves tf to `save_end`.
+        sol = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = 0.3))
+        @test sol.t ≈ [0.0, 0.3, 0.6, 0.9, 1.0]
+    end
+
+    @testset "saveat points outside the tspan are dropped" begin
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = [-1.0, 0.5, 7.0])
+        )
+        @test sol.t == [0.0, 0.5, 1.0]
+    end
+
+    @testset "save_start / save_end / save_on" begin
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = [0.5], save_start = false)
+        )
+        @test sol.t == [0.5, 1.0]
+
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = [0.5], save_end = false)
+        )
+        @test sol.t == [0.0, 0.5]
+
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(
+                prob, ltg(); dt = 0.1, saveat = [0.5],
+                save_on = false, save_everystep = true
+            )
+        )
+        @test isempty(sol.t)
+        @test isempty(sol.u)
+    end
+
+    @testset "backward integration keeps the direction ordering" begin
+        bprob = make_prob((1.0, 0.0))
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(bprob, ltg(); dt = 0.1, saveat = [0.25, 0.75])
+        )
+        @test sol.t == [1.0, 0.75, 0.25, 0.0]
+
+        sol = DiffEqBase.solve!(DiffEqBase.init(bprob, ltg(); dt = 0.1, saveat = 0.25))
+        @test sol.t == [1.0, 0.75, 0.5, 0.25, 0.0]
+    end
+
+    @testset "asking for output does not change the trajectory" begin
+        # saveat points are interpolated, never stepped onto.
+        plain = DiffEqBase.solve!(DiffEqBase.init(prob, ltg(); dt = 0.1))
+        dense = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = 0.017, save_everystep = true)
+        )
+        @test dense.u[end] == plain.u[end]
+    end
+
+    @testset "saveat times passed as tstops are landed on, not interpolated" begin
+        # The documented way to get output accurate to the order of the scheme.
+        ts = [0.25, 0.55]
+        landed = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, tstops = ts, saveat = ts)
+        )
+        interpolated = DiffEqBase.solve!(
+            DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = ts)
+        )
+        @test issorted(landed.t)
+        @test allunique(landed.t)
+        for t in ts
+            @test t in landed.t
+            @test t in interpolated.t
+        end
+        # Same times, different values: stepped state vs. interpolated one.
+        for t in ts
+            a = landed.u[findfirst(==(t), landed.t)]
+            b = interpolated.u[findfirst(==(t), interpolated.t)]
+            @test a != b
+            # ... but both approximate the same solution.
+            @test isapprox(a, b; rtol = 1.0e-3)
+        end
+    end
+
+    @testset "add_saveat!" begin
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1)
+        DiffEqBase.step!(integrator)
+        DiffEqBase.add_saveat!(integrator, 0.65)
+        sol = DiffEqBase.solve!(integrator)
+        @test sol.t == [0.0, 0.65, 1.0]
+        # Cannot save behind the current time.
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1)
+        DiffEqBase.step!(integrator)
+        DiffEqBase.step!(integrator)
+        @test_throws ErrorException DiffEqBase.add_saveat!(integrator, 0.05)
+    end
+
+    @testset "other algorithms" begin
+        for alg in algs()
+            sol = DiffEqBase.solve!(
+                DiffEqBase.init(prob, alg; dt = 0.1, saveat = [0.5])
+            )
+            @test sol.t[1] == 0.0
+            @test sol.t[end] == 1.0
+            @test 0.5 in sol.t
+            @test sol.retcode == ReturnCode.Success
+        end
+    end
+
+    @testset "nested splitting saves only at the outer level" begin
+        inner = GenericSplitFunction((f1, f1), ([1, 2], [1, 2]))
+        outer = GenericSplitFunction((f1, inner), ([1, 2, 3], [1, 2]))
+        nprob = OperatorSplittingProblem(outer, copy(U0), TSPAN)
+        alg = LieTrotterGodunov((Euler(), LieTrotterGodunov((Euler(), Euler()))))
+        integrator = DiffEqBase.init(nprob, alg; dt = 0.1, saveat = [0.5])
+        sol = DiffEqBase.solve!(integrator)
+        @test sol.t == [0.0, 0.5, 1.0]
+        # The inner node owns no solution storage of its own.
+        @test !hasproperty(integrator.child_subintegrators[2], :saveiter)
+    end
+end
+
+@testset "interpolation" begin
+    prob = make_prob()
+
+    @testset "$(nameof(typeof(alg)))" for alg in algs()
+        integrator = DiffEqBase.init(prob, alg; dt = 0.1)
+        DiffEqBase.step!(integrator)
+        DiffEqBase.step!(integrator)
+        (; tprev, t) = integrator
+        tmid = (tprev + t) / 2
+
+        @test integrator(tprev) ≈ integrator.uprev
+        @test integrator(t) ≈ integrator.u
+        @test integrator(tmid) ≈ exact(tmid) rtol = 1.0e-2
+        @test integrator(tmid, Val{1}) ≈ rhs(exact(tmid)) rtol = 1.0e-2
+        # For an adaptive scheme `dt` is already the next proposal, so the endpoints
+        # above are only reproduced if the interpolant does not key off it.
+        SciMLBase.isadaptive(alg) && @test integrator.dt != t - tprev
+
+        out = similar(integrator.u)
+        integrator(out, tmid)
+        @test out ≈ integrator(tmid)
+        @test integrator(tmid, Val{0}; idxs = 2) ≈ integrator(tmid)[2]
+        sub = similar(integrator.u, 2)
+        integrator(sub, tmid, Val{0}; idxs = [1, 3])
+        @test sub ≈ integrator(tmid)[[1, 3]]
+
+        sol = DiffEqBase.solve!(
+            DiffEqBase.init(
+                prob, alg; dt = 0.1, saveat = [0.05], save_everystep = true
+            )
+        )
+        @test sol.u[findfirst(==(0.05), sol.t)] ≈ exact(0.05) rtol = 1.0e-2
+        for τ in (0.35, 0.62, 0.97)
+            @test sol(τ) ≈ exact(τ) rtol = 1.0e-2
+        end
+    end
+end
+
+@testset "change_t_via_interpolation!" begin
+    prob = make_prob()
+
+    @testset "moves the whole tree back to the interpolated state" begin
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1)
+        DiffEqBase.step!(integrator)
+        (; tprev, t) = integrator
+        tmid = (tprev + t) / 2
+        expected = integrator(tmid)
+
+        SciMLBase.change_t_via_interpolation!(integrator, tmid)
+        @test integrator.t == tmid
+        @test integrator.u ≈ expected
+        OS.validate_time_point(integrator)
+        for child in integrator.child_subintegrators
+            @test child.t == tmid
+        end
+        @test integrator.child_subintegrators[1].u ≈ expected[1:3]
+        @test integrator.child_subintegrators[2].u ≈ expected[1:2]
+    end
+
+    @testset "refuses to move before tprev" begin
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1)
+        DiffEqBase.step!(integrator)
+        DiffEqBase.step!(integrator)
+        @test_throws ErrorException SciMLBase.change_t_via_interpolation!(
+            integrator, integrator.tprev - 0.01
+        )
+    end
+
+    @testset "moving to the current t is a no-op" begin
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1)
+        DiffEqBase.step!(integrator)
+        u_before = copy(integrator.u)
+        SciMLBase.change_t_via_interpolation!(integrator, integrator.t)
+        @test integrator.u == u_before
+    end
+
+    @testset "accepts the Val{:false} the callback path passes" begin
+        # A `Val` of the *Symbol*, so the flag has to be decoded by dispatch.
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1)
+        DiffEqBase.step!(integrator)
+        tmid = (integrator.tprev + integrator.t) / 2
+        SciMLBase.change_t_via_interpolation!(integrator, tmid, Val{:false}, nothing)
+        @test integrator.t == tmid
+    end
+
+    @testset "modify_save_endpoint rewrites the saved endpoint" begin
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, save_everystep = true)
+        DiffEqBase.step!(integrator)
+        nsaved = integrator.saveiter
+        tmid = (integrator.tprev + integrator.t) / 2
+        SciMLBase.change_t_via_interpolation!(integrator, tmid, Val{true})
+        @test integrator.saveiter == nsaved + 1
+        @test integrator.sol.t[integrator.saveiter] == tmid
+    end
+end
+
+@testset "solution object plumbing" begin
+    prob = make_prob()
+
+    @testset "dense is accepted and does not reach the leaves" begin
+        integrator = DiffEqBase.init(
+            prob, ltg(); dt = 0.1, dense = true, save_everystep = true
+        )
+        sol = DiffEqBase.solve!(integrator)
+        @test sol.retcode == ReturnCode.Success
+        @test sol(0.35) ≈ (sol.u[4] .+ sol.u[5]) ./ 2
+        # `dense` must be swallowed by `__init`; a leaf handed `dense = true` would
+        # store interpolation data for stages that are not solution points.
+        @test integrator.child_subintegrators[1].sol.dense == false
+    end
+
+    @testset "sol.stats is a real DEStats" begin
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1)
+        @test integrator.sol.stats !== nothing
+        @test integrator.sol.stats.ncondition == 0
+    end
+
+    @testset "postamble! is idempotent" begin
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1)
+        sol = DiffEqBase.solve!(integrator)
+        n = length(sol.t)
+        SciMLBase.postamble!(integrator)
+        SciMLBase.postamble!(integrator)
+        @test length(integrator.sol.t) == n
+    end
+
+    @testset "postamble! does not run once per step" begin
+        finalized = Ref(0)
+        cb = DiscreteCallback(
+            (u, t, integrator) -> false,
+            integrator -> nothing;
+            finalize = (c, u, t, integrator) -> (finalized[] += 1)
+        )
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)
+        DiffEqBase.solve!(integrator)
+        @test finalized[] == 1
+    end
+end
+
+@testset "reinit!" begin
+    prob = make_prob()
+
+    @testset "restores the saved solution" begin
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, saveat = [0.5])
+        first_run = copy(DiffEqBase.solve!(integrator).t)
+        DiffEqBase.reinit!(integrator)
+        @test integrator.saveiter == 1
+        @test integrator.sol.t[1] == 0.0
+        DiffEqBase.solve!(integrator)
+        @test integrator.sol.t == first_run
+    end
+
+    @testset "reinit_callbacks = false works without a saving callback" begin
+        # Regression: this used to index an empty `discrete_callbacks` tuple.
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1)
+        DiffEqBase.solve!(integrator)
+        DiffEqBase.reinit!(integrator; reinit_callbacks = false)
+        @test integrator.t == 0.0
+        @test integrator.saveiter == 1
+        @test DiffEqBase.solve!(integrator).retcode == ReturnCode.Success
+    end
+
+    @testset "a new u0 reaches the whole tree" begin
+        # Regression: a leaf's `reinit!` restores the `u0` slice captured when it was
+        # built, and the palindromic schemes skip the forward sync of their first
+        # child, so this used to integrate a stale child state.
+        u0b = [10.0, 20.0, 30.0]
+        probb = OperatorSplittingProblem(fsplit, copy(u0b), TSPAN)
+        for alg in (
+                ltg(),
+                StrangMarchuk((Euler(), Euler())),
+                StrangMarchuk((Tsit5(), Euler())),
+                PalindromicPairLieTrotterGodunov((Tsit5(), Tsit5())),
+            )
+            integrator = DiffEqBase.init(make_prob(), alg; dt = 0.1)
+            DiffEqBase.solve!(integrator)
+            DiffEqBase.reinit!(integrator, copy(u0b))
+            @test integrator.child_subintegrators[1].u == u0b
+            @test integrator.child_subintegrators[2].u == u0b[1:2]
+            got = DiffEqBase.solve!(integrator).u[end]
+            ref = DiffEqBase.solve!(DiffEqBase.init(probb, alg; dt = 0.1)).u[end]
+            @test got ≈ ref rtol = 1.0e-12
+        end
+    end
+
+    @testset "erase_sol clears the stored solution" begin
+        integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, save_everystep = true)
+        DiffEqBase.solve!(integrator)
+        DiffEqBase.reinit!(integrator; erase_sol = true)
+        @test integrator.saveiter == 1
+        @test length(integrator.sol.t) == 1
+    end
+end
+
+@testset "callbacks can be attached" begin
+    # Regression: `init` used to throw as soon as any callback was passed.
+    prob = make_prob()
+    calls = Ref(0)
+    cb = DiscreteCallback((u, t, integrator) -> false, integrator -> (calls[] += 1))
+    integrator = DiffEqBase.init(prob, ltg(); dt = 0.1, callback = cb)
+    @test integrator.derivative_discontinuity == false
+    sol = DiffEqBase.solve!(integrator)
+    @test sol.retcode == ReturnCode.Success
+    @test calls[] == 0
+
+    integrator = DiffEqBase.init(
+        prob, ltg(); dt = 0.1,
+        callback = CallbackSet(cb, cb)
+    )
+    @test DiffEqBase.solve!(integrator).retcode == ReturnCode.Success
+end