Base.@constprop :aggressive function SciMLBase.__solve(
        prob::Union{
            SciMLBase.AbstractODEProblem,
            SciMLBase.AbstractDAEProblem,
        },
        alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}, args...;
        kwargs...
    )
    integrator = SciMLBase.__init(prob, alg, args...; kwargs...)
    solve!(integrator)
    return integrator.sol
end

determine_controller_datatype(u::AbstractVector{<:Number}, internalnorm, ts::Tuple{<:Number, <:Number}) = promote_type(typeof(SciMLBase.value(internalnorm(u, ts[1]))), typeof(SciMLBase.value(internalnorm(u, ts[2]))), eltype(SciMLBase.value.(ts)))
determine_controller_datatype(u, internalnorm, ts::Tuple{<:Number, <:Number}) = promote_type(typeof(SciMLBase.value(ts[1])), typeof(SciMLBase.value(ts[2]))) # This seems to be an assumption implicitly taken somewhere
determine_controller_datatype(u::AbstractVector{<:Number}, internalnorm, ts::Tuple{<:Integer, <:Integer}) = promote_type(typeof(SciMLBase.value(internalnorm(u, ts[1]))), typeof(SciMLBase.value(internalnorm(u, ts[2]))), eltype(float.(SciMLBase.value(ts))))
determine_controller_datatype(u, internalnorm, ts::Tuple{<:Integer, <:Integer}) = promote_type(typeof(float(SciMLBase.value(ts[1]))), typeof(float(SciMLBase.value(ts[2])))) # This seems to be an assumption implicitly taken somewhere

mutable struct zero_func_struct{u1Type, uType, tType, kType, CacheType, idxsType, varsType, callbackType, outType, outCacheType, FunctionType, tType2, ParameterType}
    u₁::u1Type
    callback::callbackType
    dt::tType
    uprev::uType
    u::uType
    k::kType
    cache::CacheType
    idxs::idxsType
    differential_vars::varsType
    ind::Int
    out::outType
    out_low::outCacheType
    out_high::outCacheType
    f::FunctionType
    tprev::tType2
    p::ParameterType
end

parameter_values(z::zero_func_struct) = z.p

function (z::zero_func_struct)(θ, p)
    iszero(θ) && return z.out_low[z.ind]
    isone(θ) && return z.out_high[z.ind]
    ode_interpolant!(z.u₁, θ, z.dt, z.uprev, z.u, z.k, z.cache, z.idxs, Val{0}, z.differential_vars)
    return zero_condition(z.callback, z.out, z.u₁, z.tprev + θ * z.dt, z, z.ind)
end

@inline zero_condition(cb::ContinuousCallback, out::Nothing, u, t, z, ind) = cb.condition(u, t, z)
@inline function zero_condition(cb::VectorContinuousCallback, out, u, t, z, ind)
    cb.condition(out, u, t, z)
    return out[ind]
end

Base.@constprop :aggressive function SciMLBase.__init(
        prob::Union{
            SciMLBase.AbstractODEProblem,
            SciMLBase.AbstractDAEProblem,
            SciMLBase.AbstractRODEProblem,
            SciMLBase.AbstractJumpProblem,
        },
        alg::Union{
            OrdinaryDiffEqAlgorithm,
            DAEAlgorithm,
            StochasticDiffEqAlgorithm,
            StochasticDiffEqRODEAlgorithm,
        },
        timeseries_init = (),
        ts_init = (),
        ks_init = ();
        kwargs...
    )
    return _ode_init(prob, alg, timeseries_init, ts_init, ks_init; kwargs...)
end

"""
    resolve_stage_step_limiters(alg, stage_limiter, step_limiter, verbose_spec)

Resolve the effective `(stage_limiter!, step_limiter!)` from the solve-level
`stage_limiter`/`step_limiter` keywords. A non-trivial per-algorithm
`stage_limiter!`/`step_limiter!` field is still honored (with a deprecation
warning) when the matching keyword is not supplied. `alg` should be the concrete
method (e.g. unwrapped from `MethodOfSteps`). Supplying `stage_limiter` to a
method for which [`has_stage_limiter`](@ref) is `false` triggers the
`stage_limiter_unused` verbosity toggle, which defaults to `ErrorLevel` (it
errors) and can be lowered to `WarnLevel`/`Silent` to allow it.
"""
function resolve_stage_step_limiters(alg, stage_limiter, step_limiter, verbose_spec)
    stage_limiter! = stage_limiter
    step_limiter! = step_limiter
    if stage_limiter === trivial_limiter!
        if has_stage_limiter(alg) && hasproperty(alg, :stage_limiter!) &&
                alg.stage_limiter! !== trivial_limiter!
            Base.depwarn(
                "Passing `stage_limiter!` to the algorithm constructor is deprecated; " *
                    "pass `stage_limiter` as a keyword argument to `solve`/`init` instead.",
                :stage_limiter,
            )
            stage_limiter! = alg.stage_limiter!
        end
    elseif !has_stage_limiter(alg)
        @SciMLMessage(
            "`stage_limiter` was supplied, but $(nameof(typeof(alg))) does not apply stage " *
                "limiters (it does not opt into `OrdinaryDiffEqCore.has_stage_limiter`). This " *
                "errors by default; set the `stage_limiter_unused` verbosity toggle to " *
                "`WarnLevel()` or `Silent()` to allow it. Use `step_limiter` for a limiter " *
                "applied once per accepted step, which every method supports.",
            verbose_spec, :stage_limiter_unused
        )
    end
    if step_limiter === trivial_limiter! && hasproperty(alg, :step_limiter!) &&
            alg.step_limiter! !== trivial_limiter!
        Base.depwarn(
            "Passing `step_limiter!` to the algorithm constructor is deprecated; " *
                "pass `step_limiter` as a keyword argument to `solve`/`init` instead.",
            :step_limiter,
        )
        step_limiter! = alg.step_limiter!
    end
    return stage_limiter!, step_limiter!
end

"""
    _ode_init(prob, alg, timeseries_init=(), ts_init=(), ks_init=(); kwargs...)

Internal implementation of `__init` for ODE/DAE/SDE/RODE problems. This is
separated from `__init` so that SDE packages can call it directly, bypassing
method dispatch (which would otherwise re-enter SDE's more specific `__init`).
"""
Base.@constprop :aggressive function _ode_init(
        prob,
        alg,
        timeseries_init = (),
        ts_init = (),
        ks_init = ();
        saveat = (),
        tstops = (),
        d_discontinuities = (),
        save_idxs = nothing,
        save_everystep = isempty(saveat),
        save_on = true,
        save_discretes = true,
        save_start = save_everystep || isempty(saveat) ||
            saveat isa Number || prob.tspan[1] in saveat,
        save_end = nothing,
        callback = nothing,
        dense = save_everystep && isempty(saveat) &&
            !default_linear_interpolation(prob, alg),
        calck = (callback !== nothing && callback !== CallbackSet()) ||
            (dense) || !isempty(saveat), # and no dense output
        dt = nothing,
        # For runtime-unit quantities (DynamicQuantities), eltype(prob.tspan)(0) would
        # drop units; use a value-based zero to preserve units.
        dtmin = zero(prob.tspan[1]),
        dtmax = (prob.tspan[end] - prob.tspan[1]),
        force_dtmin = false,
        adaptive = anyadaptive(alg),
        abstol = nothing,
        reltol = nothing,
        controller = nothing,
        fullnormalize = true,
        failfactor = 2,
        maxiters = anyadaptive(alg) ? 1000000 : typemax(Int),
        internalnorm = ODE_DEFAULT_NORM,
        internalopnorm = opnorm,
        stage_limiter = trivial_limiter!,
        step_limiter = trivial_limiter!,
        isoutofdomain = ODE_DEFAULT_ISOUTOFDOMAIN,
        unstable_check = ODE_DEFAULT_UNSTABLE_CHECK,
        verbose = Standard(),
        timeseries_errors = true,
        dense_errors = false,
        advance_to_tstop = false,
        stop_at_next_tstop = false,
        initialize_save = true,
        progress = false,
        progress_steps = 1000,
        progress_name = "ODE",
        progress_message = ODE_DEFAULT_PROG_MESSAGE,
        progress_id = :OrdinaryDiffEq,
        userdata = nothing,
        allow_extrapolation = alg_extrapolates(alg),
        initialize_integrator = true,
        alias = ODEAliasSpecifier(),
        initializealg = DefaultInit(),
        rng = nothing,
        # SDE/RODE fields: accepted here so that SDE packages can delegate to
        # _ode_init and construct an ODEIntegrator with noise populated.
        save_noise = false,
        delta = nothing,
        W = nothing,
        P = nothing,
        sqdt = nothing,
        noise = nothing,
        c = nothing,
        rate_constants = nothing,
        # Pre-built cache for SDE delegation (skip alg_cache call)
        _cache = nothing,
        # Pre-built u/uprev for SDE delegation (cache holds references to these)
        _u = nothing,
        _uprev = nothing,
        seed = UInt64(0),
        kwargs...
    )
    # ODE/DAE-specific validation (skip for RODE/SDE problems)
    if !(prob isa SciMLBase.AbstractRODEProblem)
        if prob isa SciMLBase.AbstractDAEProblem && alg isa OrdinaryDiffEqAlgorithm
            error("You cannot use an ODE Algorithm with a DAEProblem")
        end

        if prob isa SciMLBase.AbstractODEProblem && alg isa DAEAlgorithm
            error("You cannot use an DAE Algorithm with a ODEProblem")
        end

        if prob isa SciMLBase.ODEProblem
            if !(prob.f isa SciMLBase.DynamicalODEFunction) && alg isa PartitionedAlgorithm
                error("You can not use a solver designed for partitioned ODE with this problem. Please choose a solver suitable for your problem")
            end
        end

        if prob.f isa DynamicalODEFunction && prob.f.mass_matrix isa Tuple
            if any(mm != I for mm in prob.f.mass_matrix)
                error("This solver is not able to use mass matrices. For compatible solvers see https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/")
            end
        elseif !(prob isa SciMLBase.AbstractDiscreteProblem) &&
                !(prob isa SciMLBase.AbstractDAEProblem) &&
                !is_mass_matrix_alg(alg) &&
                prob.f.mass_matrix != I
            error("This solver is not able to use mass matrices. For compatible solvers see https://docs.sciml.ai/DiffEqDocs/stable/solvers/dae_solve/")
        end
    end

    verbose_spec = _process_verbose_param(verbose)

    # `stage_limiter`/`step_limiter` are solver-level options stored in `opts`: the
    # stage limiter is threaded through `opts` and applied per stage inside
    # `perform_step!`, while the step limiter is applied centrally on each accepted
    # step from the loop footer. The per-algorithm `stage_limiter!`/`step_limiter!`
    # constructor fields are deprecated in favor of these keywords.
    stage_limiter!, step_limiter! = resolve_stage_step_limiters(
        alg, stage_limiter, step_limiter, verbose_spec
    )

    if alg isa OrdinaryDiffEqRosenbrockAdaptiveAlgorithm &&
            # https://github.com/SciML/OrdinaryDiffEq.jl/pull/2079 fixes this for Rosenbrock23 and 32
            !only_diagonal_mass_matrix(alg) &&
            prob.f.mass_matrix isa AbstractMatrix &&
            all(isequal(0), prob.f.mass_matrix)
        # technically this should also warn for zero operators but those are hard to check for
        if (dense || !isempty(saveat))
            @SciMLMessage(
                "Rosenbrock methods on equations without differential states do not bound the error on interpolations.",
                verbose_spec, :rosenbrock_no_differential_states
            )
        end
    end

    if only_diagonal_mass_matrix(alg) &&
            prob.f.mass_matrix isa AbstractMatrix &&
            !isdiag(prob.f.mass_matrix)
        throw(ArgumentError("$(typeof(alg).name.name) only works with diagonal mass matrices. Please choose a solver suitable for your problem (e.g. Rodas5P)"))
    end

    if !isempty(saveat) && dense
        @SciMLMessage(
            "Dense output is incompatible with saveat. Please use the SavingCallback from the Callback Library to mix the two behaviors.",
            verbose_spec, :dense_output_saveat
        )
    end

    progress && @logmsg(LogLevel(-1), progress_name, _id = progress_id, progress = 0)

    tType = eltype(prob.tspan)
    tspan = prob.tspan
    tdir = sign(tspan[end] - tspan[1])

    t = tspan[1]

    if (!isadaptive(alg) || !adaptive) &&
            isnothing(dt) && isempty(tstops) && dt_required(alg)
        throw(ArgumentError("Fixed timestep methods require a choice of dt or choosing the tstops"))
    end
    if !isadaptive(alg) && adaptive
        throw(ArgumentError("Fixed timestep methods can not be run with adaptive=true"))
    end

    isdae = if !isnothing(W)
        false # RODE/SDE — never DAE
    else
        alg isa DAEAlgorithm || (
            !(prob isa SciMLBase.AbstractDiscreteProblem) &&
                prob.f.mass_matrix != I &&
                !(prob.f.mass_matrix isa Tuple) &&
                ArrayInterface.issingular(prob.f.mass_matrix)
        )
    end
    if alg isa CompositeAlgorithm && alg.choice_function isa AutoSwitch
        auto = alg.choice_function
        _alg = CompositeAlgorithm(
            alg.algs,
            AutoSwitchCache(
                0, 0,
                auto.nonstiffalg,
                auto.stiffalg,
                auto.stiffalgfirst,
                auto.maxstiffstep,
                auto.maxnonstiffstep,
                auto.nonstifftol,
                auto.stifftol,
                auto.dtfac,
                auto.stiffalgfirst,
                auto.switch_max, 0
            )
        )
    else
        _alg = alg
    end

    alias
    if !(alias isa ODEAliasSpecifier)
        throw(ArgumentError("alias kwarg must be an ODEAliasSpecifier"))
    end

    if isnothing(alias.alias_f) || alias.alias_f
        f = prob.f
    else
        f = deepcopy(prob.f)
    end

    if isnothing(alias.alias_p) || alias.alias_p
        p = prob.p
    else
        p = recursivecopy(prob.p)
    end

    if !isnothing(alias.alias_u0) && alias.alias_u0
        u = prob.u0
    else
        u = recursivecopy(prob.u0)
    end

    # Null u0 (e.g., MTK systems with only callbacks and no state variables).
    # Most solvers' caches are built with zero(u)/similar(u) and can't ingest `nothing`,
    # so coerce to Float64[] to let them limp through. Solvers that can handle the
    # `nothing` signal directly (e.g., IDSolve bypasses NonlinearSolve entirely)
    # opt in via allows_null_u0 and keep the richer signal.
    if u === nothing && !allows_null_u0(_alg)
        u = Float64[]
    end

    # SDE delegation: use pre-built u if provided (cache holds references to it)
    if _u !== nothing
        u = _u
    end

    if _alg isa DAEAlgorithm
        if !isnothing(alias.alias_du0) && alias.alias_du0
            du = prob.du0
        else
            du = recursivecopy(prob.du0)
        end
        duprev = recursivecopy(du)
    else
        du = nothing
        duprev = nothing
    end

    uType = typeof(u)
    uBottomEltype = recursive_bottom_eltype(u)
    uBottomEltypeNoUnits = recursive_unitless_bottom_eltype(u)

    uEltypeNoUnits = recursive_unitless_eltype(u)
    tTypeNoUnits = typeof(DiffEqBase.stripunits(oneunit(first(tspan))))

    scalar_type_tol =
        uBottomEltypeNoUnits == uBottomEltype &&
        uBottomEltype <: Union{Real, Complex}

    if prob isa SciMLBase.AbstractDiscreteProblem && abstol === nothing
        abstol_internal = false
    elseif abstol === nothing
        if scalar_type_tol
            abstol_internal = unitfulvalue(real(convert(uBottomEltype, oneunit(uBottomEltype) * 1 // 10^6)))
        else
            abstol_internal = unitfulvalue.(real.(oneunit.(u) .* 1 // 10^6))
        end
    else
        abstol_internal = real.(abstol)
    end

    if prob isa SciMLBase.AbstractDiscreteProblem && reltol === nothing
        reltol_internal = false
    elseif reltol === nothing
        if scalar_type_tol
            reltol_internal = unitfulvalue(real(convert(uBottomEltype, oneunit(uBottomEltype) * 1 // 10^3)))
        else
            reltol_internal = unitfulvalue.(real.(oneunit.(u) .* 1 // 10^3))
        end
    else
        reltol_internal = real.(reltol)
    end

    dtmax > zero(dtmax) && tdir < 0 && (dtmax *= tdir) # Allow positive dtmax, but auto-convert
    # dtmin is all abs => does not care about sign already.

    if !isdae && isinplace(prob) && u isa AbstractArray && eltype(u) <: Number &&
            uBottomEltypeNoUnits == uBottomEltype && tType == tTypeNoUnits # Could this be more efficient for other arrays?
        # SDE delegation provides `_cache` and a non-nothing `W`. In that case
        # `rate_prototype`'s value is unused: only its type is taken downstream,
        # the alg_cache calls are skipped (cache is provided), get_fsalfirstlast
        # is short-circuited, and the [rate_prototype] push is gated by
        # `isnothing(W)` which is false for SDE. Alias instead of allocating.
        rate_prototype = (_cache !== nothing && W !== nothing) ? u : recursivecopy(u)
    elseif prob isa DAEProblem
        rate_prototype = prob.du0
    else
        if (uBottomEltypeNoUnits == uBottomEltype && tType == tTypeNoUnits) ||
                eltype(u) <: Enum
            rate_prototype = u
        else # has units!
            rate_prototype = u / oneunit(first(tspan))
        end
    end
    rateType = typeof(rate_prototype) ## Can be different if united

    res_prototype = nothing
    if isdae
        if uBottomEltype == uBottomEltypeNoUnits
            res_prototype = u
        else
            res_prototype = one(u)
        end
        resType = typeof(res_prototype)
    end

    if isnothing(alias.alias_tstops) || alias.alias_tstops
        tstops = tstops
    else
        tstops = recursivecopy(tstops)
    end

    if tstops isa AbstractArray || tstops isa Tuple || tstops isa Number
        _tstops_cache = tstops
    else
        _tstops_cache = tstops
        tstops = ()
    end
    tstops_internal = initialize_tstops(tType, tstops, d_discontinuities, tspan)
    saveat_internal = initialize_saveat(tType, saveat, tspan)
    d_discontinuities_internal = initialize_d_discontinuities(
        tType, d_discontinuities,
        tspan
    )

    callbacks_internal = CallbackSet(callback)

    max_len_cb = DiffEqBase.max_vector_callback_length_int(callbacks_internal)
    if max_len_cb !== nothing
        uBottomEltypeReal = real(uBottomEltype)
        if isinplace(prob)
            callback_cache = DiffEqBase.CallbackCache(
                u, max_len_cb, uBottomEltypeReal,
                uBottomEltypeReal
            )
        else
            callback_cache = DiffEqBase.CallbackCache(
                max_len_cb, uBottomEltypeReal,
                uBottomEltypeReal
            )
        end
    else
        callback_cache = nothing
    end

    ### Algorithm-specific defaults ###
    save_idxs,
        saved_subsystem = SciMLBase.get_save_idxs_and_saved_subsystem(
        prob, save_idxs
    )

    ks_prototype = nothing
    if save_idxs === nothing
        ksEltype = Vector{rateType}
    else
        ks_prototype = rate_prototype[save_idxs]
        ksEltype = Vector{typeof(ks_prototype)}
    end

    # Have to convert in case passed in wrong.
    if save_idxs === nothing
        timeseries = timeseries_init === () ? uType[] :
            convert(Vector{uType}, timeseries_init)
    else
        u_initial = u[save_idxs]
        timeseries = timeseries_init === () ? typeof(u_initial)[] :
            convert(Vector{uType}, timeseries_init)
    end

    ts = ts_init === () ? tType[] : convert(Vector{tType}, ts_init)
    ks = ks_init === () ? ksEltype[] : convert(Vector{ksEltype}, ks_init)
    alg_choice = (
            _alg isa CompositeAlgorithm ||
            _alg isa StochasticDiffEqCompositeAlgorithm ||
            _alg isa StochasticDiffEqRODECompositeAlgorithm
        ) ? Int[] : nothing

    if (!adaptive || !isadaptive(_alg)) && save_everystep && tspan[2] - tspan[1] != Inf
        if isnothing(dt)
            steps = length(tstops)
        else
            # For fixed dt, the only time dtmin makes sense is if it's smaller than eps().
            # Therefore user specified dtmin doesn't matter, but we need to ensure dt>=eps()
            # to prevent infinite loops.
            abs(dt) < dtmin &&
                throw(ArgumentError("Supplied dt is smaller than dtmin"))
            steps = ceil(Int, internalnorm((tspan[2] - tspan[1]) / dt, tspan[1]))
        end
        sizehint!(timeseries, steps + 1)
        sizehint!(ts, steps + 1)
        sizehint!(ks, steps + 1)
    elseif save_everystep
        sizehint!(timeseries, 50)
        sizehint!(ts, 50)
        sizehint!(ks, 50)
    elseif !isempty(saveat_internal)
        savelength = length(saveat_internal) + 1
        if save_start == false
            savelength -= 1
        end
        if save_end == false && prob.tspan[2] in saveat_internal.valtree
            savelength -= 1
        end
        sizehint!(timeseries, savelength)
        sizehint!(ts, savelength)
        sizehint!(ks, savelength)
    else
        sizehint!(timeseries, 2)
        sizehint!(ts, 2)
        sizehint!(ks, 2)
    end

    k = rateType[]

    if _uprev !== nothing
        # SDE delegation provides a pre-built uprev (cache holds references to it).
        # Use it directly instead of allocating a recursivecopy(u) that would
        # immediately be discarded.
        uprev = _uprev
    elseif uses_uprev(_alg, adaptive) || calck
        uprev = recursivecopy(u)
    else
        # Some algorithms do not use `uprev` explicitly. In that case, we can save
        # some memory by aliasing `uprev = u`, e.g. for "2N" low storage methods.
        uprev = u
    end
    if allow_extrapolation
        uprev2 = recursivecopy(u)
    else
        uprev2 = uprev
    end

    _dt = if isnothing(dt)
        isdiscretealg(alg) && isempty(tstops) ?
            eltype(prob.tspan)(1) : eltype(prob.tspan)(0)
    else
        dt
    end
    if _cache !== nothing
        cache = _cache
    elseif prob isa DAEProblem
        cache = alg_cache(
            _alg, du, u, res_prototype, rate_prototype, uEltypeNoUnits,
            uBottomEltypeNoUnits, tTypeNoUnits, uprev, uprev2, f, t, _dt,
            reltol_internal, p, calck, Val(isinplace(prob)), verbose_spec
        )
    else
        cache = alg_cache(
            _alg, u, rate_prototype, uEltypeNoUnits, uBottomEltypeNoUnits,
            tTypeNoUnits, uprev, uprev2, f, t, _dt, reltol_internal, p, calck,
            Val(isinplace(prob)), verbose_spec
        )
    end

    QT = determine_controller_datatype(u, internalnorm, tspan)

    if controller === nothing
        controller = default_controller(QT, alg)
    end

    EEstT = if tTypeNoUnits <: Integer
        QT
    elseif prob isa SciMLBase.AbstractDiscreteProblem
        constvalue(tTypeNoUnits)
    else
        typeof(internalnorm(u, t))
    end

    save_end_user = save_end
    save_end = save_end === nothing ?
        save_everystep || isempty(saveat) || saveat isa Number ||
        prob.tspan[2] in saveat : save_end

    opts = DEOptions{
        typeof(abstol_internal), typeof(reltol_internal),
        QT, tType,
        typeof(internalnorm), typeof(internalopnorm),
        typeof(save_end_user),
        typeof(callbacks_internal),
        typeof(isoutofdomain),
        typeof(progress_message), typeof(unstable_check),
        typeof(tstops_internal),
        typeof(d_discontinuities_internal), typeof(userdata),
        typeof(save_idxs),
        typeof(maxiters), typeof(_tstops_cache),
        typeof(saveat), typeof(d_discontinuities), typeof(verbose_spec),
        typeof(delta), typeof(stage_limiter!), typeof(step_limiter!),
    }(
        maxiters, save_everystep,
        adaptive, abstol_internal,
        reltol_internal,
        QT(failfactor),
        tType(dtmax), tType(dtmin),
        internalnorm,
        internalopnorm,
        save_idxs, tstops_internal,
        saveat_internal,
        d_discontinuities_internal,
        _tstops_cache, saveat,
        d_discontinuities,
        userdata, progress,
        progress_steps,
        progress_name,
        progress_message,
        progress_id,
        timeseries_errors,
        dense_errors, delta, stage_limiter!, step_limiter!, dense,
        save_on, save_start,
        save_end, save_noise, save_discretes, save_end_user,
        callbacks_internal,
        isoutofdomain,
        unstable_check,
        verbose_spec, calck, force_dtmin,
        advance_to_tstop,
        stop_at_next_tstop
    )

    stats = SciMLBase.DEStats(0)
    differential_vars = if prob isa DAEProblem
        prob.differential_vars
    elseif !isnothing(W)
        nothing
    else
        get_differential_vars(f, u)
    end

    # SDE/RODE solvers never populate derivative stages (ks) — RODESolution has
    # no .k field, so _has_ks returns false and save_dense_at_t! is a no-op.
    # Force dense=false so ode_interpolation uses linear interpolation instead
    # of trying to access the empty ks vector.
    if !isnothing(W)
        dense = false
    end

    # Select the interval-search strategy for scalar interpolation from what
    # the solve options already guarantee about the saved time grid: fixed-dt
    # stepping and pure range-`saveat` saving produce (near-)uniform grids
    # where interpolation search is optimal and no re-probing is needed.
    # Fully adaptive grids start on the robust hinted gallop and re-select on
    # grid growth and at the ending phase (`_finalize_ts_hint!`).
    ts_hint = if !adaptive || (saveat isa AbstractRange && length(saveat) > 1 && !save_everystep)
        TsSearchHint(ts, 0, KIND_INTERPOLATION_SEARCH, typemax(Int))
    else
        TsSearchHint(ts)
    end
    id = InterpolationData(
        f, timeseries, ts, ks, alg_choice, dense, cache, differential_vars, false,
        ts_hint
    )

    _sol_kwargs = if !isnothing(W)
        (;
            W = W, seed = seed, interp = id, dense = dense, alg_choice = alg_choice,
            calculate_error = false, stats = stats, saved_subsystem = saved_subsystem,
        )
    else
        (;
            dense = dense, k = ks, interp = id, alg_choice = alg_choice,
            calculate_error = false, stats = stats, saved_subsystem = saved_subsystem,
        )
    end
    sol = SciMLBase.build_solution(prob, _alg, ts, timeseries; _sol_kwargs...)

    FType = typeof(f)
    SolType = typeof(sol)
    cacheType = typeof(cache)

    # rate/state = (state/time)/state = 1/t units, internalnorm drops units
    # we don't want to differentiate through eigenvalue estimation
    eigen_est = inv(one(tType))
    tprev = t
    dtcache = tType(_dt)
    dtpropose = tType(_dt)
    iter = 0
    kshortsize = 0
    reeval_fsal = false
    derivative_discontinuity = false
    EEst = oneunit(EEstT) # https://github.com/JuliaPhysics/Measurements.jl/pull/135
    just_hit_tstop = false
    next_step_tstop = false
    tstop_target = zero(t)
    isout = false
    accept_step = false
    force_stepfail = false
    last_stepfail = false
    do_error_check = true
    event_last_time = 0
    vector_event_last_time = 1
    last_event_error = prob isa SciMLBase.AbstractDiscreteProblem ? false :
        (
            Base.isbitstype(uBottomEltypeNoUnits) ? zero(uBottomEltypeNoUnits) :
            0.0
        )
    dtchangeable = isdtchangeable(_alg)
    success_iter = 0
    reinitialize = true
    saveiter = 0 # Starts at 0 so first save is at 1
    saveiter_dense = 0
    fsalfirst, fsallast = _cache !== nothing ? (nothing, nothing) :
        get_fsalfirstlast(cache, rate_prototype)

    _rng = rng === nothing ? Random.default_rng() : rng

    num_probs = 0
    for i in callbacks_internal.continuous_callbacks
        if i.maybe_discontinuity
            num_probs += 1
        end
    end
    disco_probs = Vector{IntervalNonlinearProblem}(undef, num_probs)
    idx = 1
    for i in callbacks_internal.continuous_callbacks
        if i.maybe_discontinuity
            u₁ = (u isa AbstractArray) ? similar(u) : zero(u)
            out, out_low, out_high = if i isa VectorContinuousCallback
                arr = (u isa AbstractArray) ? similar(u, i.len) : zeros(typeof(u), i.len)
                arr, similar(arr), similar(arr)
            else
                nothing, zeros(Float64, 1), zeros(Float64, 1)
            end
            zero_func = zero_func_struct(u₁, i, _dt, uprev, u, k, cache, save_idxs, differential_vars, 1, out, out_low, out_high, f, tprev, p)
            zero_func_wrapped = FunctionWrapper{Float64, Tuple{Float64, typeof(p)}}(zero_func)
            disco_prob = IntervalNonlinearProblem{false}(zero_func_wrapped, [zero(tType), one(tType)], p)
            disco_probs[idx] = disco_prob
            idx += 1
        end
    end

    if num_probs > 0
        disco_probs = convert(Vector{typeof(disco_probs[1])}, disco_probs)
    end

    controller_cache = setup_controller_cache(_alg, cache, controller, EEstT, disco_probs)

    is_disco_step = false
    disco_checkpoint = zero(t)
    # Seed the initial EEst on the controller cache (was previously
    # `integrator.EEst = oneunit(EEstT)`).
    set_EEst!(controller_cache, EEst)

    integrator = ODEIntegrator{
        typeof(_alg), isinplace(prob), uType, typeof(du),
        tType, typeof(p), typeof(eigen_est),
        typeof(tdir), typeof(k), SolType,
        FType, cacheType,
        typeof(opts), typeof(fsalfirst),
        typeof(last_event_error), typeof(callback_cache),
        typeof(initializealg), typeof(differential_vars),
        typeof(controller_cache), typeof(_rng),
        typeof(W), typeof(P), typeof(sqdt),
        typeof(noise), typeof(c), typeof(rate_constants),
    }(
        sol, u, du, k, t, tType(_dt), f, p,
        uprev, uprev2, duprev, tprev,
        _alg, dtcache, dtchangeable,
        dtpropose, tdir, eigen_est,
        controller_cache,
        success_iter,
        iter, saveiter, saveiter_dense, cache,
        callback_cache,
        kshortsize, force_stepfail,
        last_stepfail,
        just_hit_tstop, next_step_tstop, tstop_target, do_error_check,
        event_last_time,
        vector_event_last_time,
        last_event_error, accept_step,
        isout, reeval_fsal,
        derivative_discontinuity, reinitialize, isdae,
        opts, stats, initializealg, differential_vars,
        fsalfirst, fsallast, _rng, is_disco_step, disco_checkpoint,
        W, P, sqdt,
        noise, c, rate_constants
    )

    if initialize_integrator
        if isdae || SciMLBase.has_initializeprob(prob.f) ||
                prob isa SciMLBase.ImplicitDiscreteProblem
            SciMLBase.initialize_dae!(integrator)
            !isnothing(integrator.u) && update_uprev!(integrator)
        end

        if save_start
            integrator.saveiter += 1 # Starts at 1 so first save is at 2
            integrator.saveiter_dense += 1
            copyat_or_push!(ts, 1, t)
            # N.B.: integrator.u can be modified by initialized_dae!
            if save_idxs === nothing
                copyat_or_push!(timeseries, 1, integrator.u)
                if isnothing(W)
                    copyat_or_push!(ks, 1, [rate_prototype])
                end
            else
                copyat_or_push!(timeseries, 1, integrator.u[save_idxs], Val{false})
                if isnothing(W)
                    copyat_or_push!(ks, 1, [ks_prototype])
                end
            end
        else
            integrator.saveiter = 0 # Starts at 0 so first save is at 1
            integrator.saveiter_dense = 0
        end

        initialize_callbacks!(integrator, initialize_save)
        initialize!(integrator, integrator.cache)

        if _alg isa OrdinaryDiffEqCompositeAlgorithm
            # in case user mixes adaptive and non-adaptive algorithms
            ensure_behaving_adaptivity!(integrator, integrator.cache)
        end
        if alg_choice !== nothing && save_start
            # Loop to get all of the extra possible saves in callback initialization
            for i in 1:(integrator.saveiter)
                copyat_or_push!(alg_choice, i, integrator.cache.current)
            end
        end
    end

    if !(_tstops_cache isa AbstractArray || _tstops_cache isa Tuple || _tstops_cache isa Number)
        tstops = _tstops_cache(parameter_values(integrator), prob.tspan)
        for tstop in tstops
            add_tstop!(integrator, tstop)
        end
    end

    handle_dt!(integrator, dt)

    # Noise process initialization for SDE/RODE integrators
    if !isnothing(integrator.W)
        modify_dt_for_tstops!(integrator)
        integrator.sqdt = integrator.tdir * sqrt(abs(integrator.dt))
        integrator.W.dt = integrator.dt
        !isnothing(integrator.P) && (integrator.P.dt = integrator.dt)
    end

    # Starting-time discontinuity: if `d_discontinuities` names `t0`, advance
    # past it so the first step's RHS evaluations are on the post-discontinuity
    # side of `t`-dependent branches in `f` (e.g. `if t > t0`). The regular
    # `update_fsal!` mechanism only fires at iter > 0, so handle the iter==0
    # case explicitly here. Interior/end discontinuities are handled normally
    # via the tstops heap and `update_fsal!`. Gated on `initialize_integrator`
    # because otherwise the FSAL cache is not yet populated.
    if initialize_integrator
        handle_starting_time_discontinuity!(integrator)
    end

    return integrator
end

"""
    handle_starting_time_discontinuity!(integrator)

If `integrator.opts.d_discontinuities` has an entry at exactly `integrator.t`
(i.e. at `t0`), pop it, shift `t` by one ULP in the `tdir` direction, and
re-evaluate the FSAL cache on the post-discontinuity side. No-op otherwise.
"""
function handle_starting_time_discontinuity!(integrator)
    if has_discontinuity(integrator) &&
            first_discontinuity(integrator) == integrator.tdir * integrator.t
        handle_discontinuities!(integrator)
        shift_past_discontinuity!(integrator)
        get_current_isfsal(integrator.alg, integrator.cache) &&
            reset_fsal!(integrator)
    end
    return nothing
end

function SciMLBase.solve!(integrator::ODEIntegrator)
    @inbounds while !isempty(integrator.opts.tstops)
        first_tstop = first(integrator.opts.tstops)
        while integrator.tdir * integrator.t < first_tstop
            loopheader!(integrator)
            if integrator.do_error_check && check_error!(integrator) != ReturnCode.Success
                return integrator.sol
            end

            # Use special tstop handling if flag is set, otherwise normal stepping
            if integrator.next_step_tstop
                handle_tstop_step!(integrator)
            else
                perform_step!(integrator, integrator.cache)
            end

            should_exit = integrator.next_step_tstop

            loopfooter!(integrator)
            if isempty(integrator.opts.tstops) || should_exit
                break
            end
        end
        handle_tstop!(integrator)
    end
    postamble!(integrator)

    f = integrator.sol.prob.f
    # SDE split problems may store f as a Tuple; unwrap to check for analytic solution
    f = f isa Tuple ? f[1] : f

    if SciMLBase.has_analytic(f)
        SciMLBase.calculate_solution_errors!(
            integrator.sol;
            timeseries_errors = integrator.opts.timeseries_errors,
            dense_errors = integrator.opts.dense_errors
        )
    end
    if integrator.sol.retcode != ReturnCode.Default
        return integrator.sol
    end
    return integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, ReturnCode.Success)
end

# Helpers

function handle_dt!(integrator)
    return if iszero(integrator.dt) && integrator.opts.adaptive
        auto_dt_reset!(integrator)
        if sign(integrator.dt) != integrator.tdir && !iszero(integrator.dt) &&
                !isnan(integrator.dt)
            error("Automatic dt setting has the wrong sign. Exiting. Please report this error.")
        end
        if isnan(integrator.dt)
            @SciMLMessage(
                "Automatic dt set the starting dt as NaN, causing instability. Exiting.",
                integrator.opts.verbose, :dt_NaN
            )
        end
    elseif integrator.opts.adaptive && integrator.dt > zero(integrator.dt) &&
            integrator.tdir < 0
        integrator.dt *= integrator.tdir # Allow positive dt, but auto-convert
    end
end
function handle_dt!(integrator, dt)
    return if isnothing(dt) && iszero(integrator.dt) && integrator.opts.adaptive
        auto_dt_reset!(integrator)
        if sign(integrator.dt) != integrator.tdir && !iszero(integrator.dt) &&
                !isnan(integrator.dt)
            error("Automatic dt setting has the wrong sign. Exiting. Please report this error.")
        end
        if isnan(integrator.dt)
            @SciMLMessage(
                "Automatic dt set the starting dt as NaN, causing instability. Exiting.",
                integrator.opts.verbose, :dt_NaN
            )
        end
    elseif integrator.opts.adaptive && integrator.dt > zero(integrator.dt) &&
            integrator.tdir < 0
        integrator.dt *= integrator.tdir # Allow positive dt, but auto-convert
    end
end

"""
    initialize_tstops(::Type{T}, tstops, d_discontinuities, tspan) -> BinaryHeap{T}

Build the internal directional time-stop queue for a solver integrator.

# Arguments

- `T::Type`: Element type of the queue.
- `tstops`: Iterable of requested stopping times.
- `d_discontinuities`: Iterable of derivative-discontinuity times.
- `tspan::Tuple`: Integration start and end times.

# Returns

- `BinaryHeap{T}`: Directional times strictly inside `tspan`, followed by the final time.

# Rules

- Times are stored multiplied by the integration direction, so `pop!` visits the next
  physical time for both forward and reverse integrations.
- Entries at or outside the initial and final bounds are discarded; the final bound is
  inserted exactly once.
- Solver authors implementing a custom `init` path must use this helper so `tstops`
  and derivative discontinuities preserve the standard ordering semantics.

!!! warning "Developer API, not user API"
    Application code should pass `tstops` and `d_discontinuities` to `solve` or
    `init`; it must not construct this queue directly.

# Example
```julia
tstops_internal = initialize_tstops(Float64, (0.25, 0.75), (), (0.0, 1.0))
```
"""
@inline function initialize_tstops(::Type{T}, tstops, d_discontinuities, tspan) where {T}
    tstops_internal = BinaryHeap{T}(FasterForward())

    t0, tf = tspan
    tdir = sign(tf - t0)
    tdir_t0 = tdir * t0
    tdir_tf = tdir * tf

    for t in tstops
        tdir_t = tdir * t
        tdir_t0 < tdir_t < tdir_tf && push!(tstops_internal, tdir_t)
    end
    for t in d_discontinuities
        tdir_t = tdir * t
        tdir_t0 < tdir_t < tdir_tf && push!(tstops_internal, tdir_t)
    end
    push!(tstops_internal, tdir_tf)

    return tstops_internal
end

function reinit_tstops!(
        ::Type{T}, tstops_internal, tstops, d_discontinuities, tspan; p = nothing
    ) where {T}
    empty!(tstops_internal)

    # Evaluate callable tstops
    _tstops = if tstops isa AbstractArray || tstops isa Tuple || tstops isa Number
        tstops
    elseif p !== nothing
        tstops(p, tspan)
    else
        ()
    end

    t0, tf = tspan
    tdir = sign(tf - t0)
    tdir_t0 = tdir * t0
    tdir_tf = tdir * tf

    for t in _tstops
        tdir_t = tdir * t
        tdir_t0 < tdir_t < tdir_tf && push!(tstops_internal, tdir_t)
    end
    for t in d_discontinuities
        tdir_t = tdir * t
        tdir_t0 < tdir_t < tdir_tf && push!(tstops_internal, tdir_t)
    end
    return push!(tstops_internal, tdir_tf)
end

"""
    initialize_saveat(::Type{T}, saveat, tspan) -> BinaryHeap{T}

Build the internal directional queue of output times for a solver integrator.

# Arguments

- `T::Type`: Element type of the queue.
- `saveat`: A positive output interval or iterable of requested output times.
- `tspan::Tuple`: Integration start and end times.

# Returns

- `BinaryHeap{T}`: Directional output times accepted by the standard `saveat` rules.

# Rules

- A scalar `saveat` is treated as a positive interval in the integration direction.
- An iterable `saveat` contributes only times strictly after the initial bound and at or
  before the final bound.
- Solver authors implementing a custom `init` path must use this helper to preserve
  forward and reverse integration semantics.

!!! warning "Developer API, not user API"
    Application code should set the `saveat` keyword on `solve` or `init`.

# Example
```julia
saveat_internal = initialize_saveat(Float64, 0.1, (0.0, 1.0))
```
"""
function initialize_saveat(::Type{T}, saveat, tspan) where {T}
    saveat_internal = BinaryHeap{T}(FasterForward())

    t0, tf = tspan
    tdir = sign(tf - t0)
    tdir_t0 = tdir * t0
    tdir_tf = tdir * tf

    if saveat isa Number
        directional_saveat = tdir * abs(saveat)
        for t in (t0 + directional_saveat):directional_saveat:tf
            push!(saveat_internal, tdir * t)
        end
    elseif !isempty(saveat)
        for t in saveat
            tdir_t = tdir * t
            tdir_t0 < tdir_t ≤ tdir_tf && push!(saveat_internal, tdir_t)
        end
    end

    return saveat_internal
end

function reinit_saveat!(::Type{T}, saveat_internal, saveat, tspan) where {T}
    empty!(saveat_internal)

    t0, tf = tspan
    tdir = sign(tf - t0)
    tdir_t0 = tdir * t0
    tdir_tf = tdir * tf

    return if saveat isa Number
        directional_saveat = tdir * abs(saveat)
        for t in (t0 + directional_saveat):directional_saveat:tf
            push!(saveat_internal, tdir * t)
        end
    elseif !isempty(saveat)
        for t in saveat
            tdir_t = tdir * t
            tdir_t0 < tdir_t ≤ tdir_tf && push!(saveat_internal, tdir_t)
        end
    end
end

"""
    initialize_d_discontinuities(::Type{T}, d_discontinuities, tspan) -> BinaryHeap{T}

Build the internal directional queue of derivative-discontinuity times.

# Arguments

- `T::Type`: Element type of the queue.
- `d_discontinuities`: Iterable of derivative-discontinuity times.
- `tspan::Tuple`: Integration start and end times; its direction determines queue order.

# Returns

- `BinaryHeap{T}`: All requested discontinuities stored in integration direction.

# Rules

- Unlike `initialize_tstops`, this helper does not filter times against `tspan`.
- Solver authors use this queue only when their initialization path supports the
  `d_discontinuities` solve keyword.

!!! warning "Developer API, not user API"
    Application code should supply `d_discontinuities` to `solve` or `init`.

# Example
```julia
discontinuities = initialize_d_discontinuities(Float64, (0.5,), (0.0, 1.0))
```
"""
function initialize_d_discontinuities(::Type{T}, d_discontinuities, tspan) where {T}
    d_discontinuities_internal = BinaryHeap{T}(FasterForward())
    sizehint!(d_discontinuities_internal, length(d_discontinuities))

    t0, tf = tspan
    tdir = sign(tf - t0)

    for t in d_discontinuities
        push!(d_discontinuities_internal, tdir * t)
    end

    return d_discontinuities_internal
end

function reinit_d_discontinuities!(::Type{T}, d_discontinuities_internal, d_discontinuities, tspan) where {T}
    empty!(d_discontinuities_internal)

    t0, tf = tspan
    tdir = sign(tf - t0)

    for t in d_discontinuities
        push!(d_discontinuities_internal, tdir * t)
    end
    return
end

function initialize_callbacks!(integrator, initialize_save = true)
    t = integrator.t
    u = integrator.u
    callbacks = integrator.opts.callback
    integrator.derivative_discontinuity = true

    derivative_discontinuity = initialize!(callbacks, u, t, integrator)

    # if the user modifies u, we need to fix previous values before initializing
    # FSAL in order for the starting derivatives to be correct
    if derivative_discontinuity
        if isinplace(integrator.sol.prob)
            recursivecopy!(integrator.uprev, integrator.u)
        else
            integrator.uprev = integrator.u
        end

        if alg_extrapolates(integrator.alg)
            if isinplace(integrator.sol.prob)
                recursivecopy!(integrator.uprev2, integrator.uprev)
            else
                integrator.uprev2 = integrator.uprev
            end
        end

        if initialize_save &&
                (
                any((c) -> c.save_positions[2], callbacks.discrete_callbacks) ||
                    any((c) -> c.save_positions[2], callbacks.continuous_callbacks)
            )
            savevalues!(integrator, true)
        end
    end

    # reset this as it is now handled so the integrators should proceed as normal
    integrator.derivative_discontinuity = false

    return if initialize_save
        SciMLBase.save_discretes_if_enabled!(integrator, integrator.opts.callback; skip_duplicates = true)
    end
end
