diff --git a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl
index fb7ab3d305..3910a77c9f 100644
--- a/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl
+++ b/lib/OrdinaryDiffEqBDF/src/OrdinaryDiffEqBDF.jl
@@ -16,7 +16,8 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!,
                            trivial_limiter!,
                            issplit, qsteady_min_default, qsteady_max_default,
                            get_current_alg_order, get_current_adaptive_order,
-                           default_controller, stepsize_controller!,
+                           default_controller_v7,
+                           legacy_default_controller, stepsize_controller!,
                            step_accept_controller!,
                            step_reject_controller!, post_newton_controller!,
                            u_modified!, DAEAlgorithm, _unwrap_val, DummyController,
diff --git a/lib/OrdinaryDiffEqBDF/src/controllers.jl b/lib/OrdinaryDiffEqBDF/src/controllers.jl
index cd65919803..fd46c2b02e 100644
--- a/lib/OrdinaryDiffEqBDF/src/controllers.jl
+++ b/lib/OrdinaryDiffEqBDF/src/controllers.jl
@@ -1,4 +1,8 @@
-function default_controller(alg::Union{QNDF, FBDF}, args...)
+function legacy_default_controller(alg::Union{QNDF, FBDF}, args...)
+    DummyController()
+end
+
+function default_controller_v7(QT, alg::Union{QNDF, FBDF}, args...)
     DummyController()
 end
 
diff --git a/lib/OrdinaryDiffEqCore/src/alg_utils.jl b/lib/OrdinaryDiffEqCore/src/alg_utils.jl
index 4539082cbb..6c2282f609 100644
--- a/lib/OrdinaryDiffEqCore/src/alg_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/alg_utils.jl
@@ -306,12 +306,26 @@ alg_maximum_order(alg::CompositeAlgorithm) = maximum(alg_order(x) for x in alg.a
 
 alg_adaptive_order(alg::Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}) = alg_order(alg) - 1
 
-# this is actually incorrect and is purposefully decreased as this tends
-# to track the real error much better
 # this is actually incorrect and is purposefully decreased as this tends
 # to track the real error much better
 
-function default_controller(alg, cache, qoldinit, _beta1 = nothing, _beta2 = nothing)
+function default_controller_v7(QT, alg)
+    if ispredictive(alg)
+        return NewPredictiveController(QT, alg)
+    elseif isstandard(alg)
+        return NewIController(QT, alg)
+    else
+        return NewPIController(QT, alg)
+    end
+end
+
+function default_controller_v7(QT, alg::CompositeAlgorithm)
+    return CompositeController(
+        map(alg->default_controller_v7(QT, alg), alg.algs)
+    )
+end
+
+function legacy_default_controller(alg, cache, qoldinit, _beta1 = nothing, _beta2 = nothing)
     if ispredictive(alg)
         return PredictiveController()
     elseif isstandard(alg)
@@ -323,6 +337,9 @@ function default_controller(alg, cache, qoldinit, _beta1 = nothing, _beta2 = not
     end
 end
 
+# TODO remove this when done
+default_controller(args...) = legacy_default_controller(args...)
+
 function _digest_beta1_beta2(alg, cache, ::Val{QT}, _beta1, _beta2) where {QT}
     if alg isa OrdinaryDiffEqCompositeAlgorithm
         beta2 = _beta2 === nothing ?
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
index 9c58c48432..983b1942dd 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
@@ -1,16 +1,83 @@
 abstract type AbstractController end
-using OrdinaryDiffEqCore
+abstract type AbstractLegacyController <: AbstractController end
+
+abstract type AbstractControllerCache end
+
+"""
+    setup_controller_cache(alg, controller::AbstractController)::AbstractControllerCache
+
+This function takes a controller together with the time stepping algorithm to
+construct and initialize the respective cache for the controller.
+"""
+setup_controller_cache
+
+"""
+    accept_step_controller(integrator, alg)::Bool
+    accept_step_controller(integrator, controller_cache::AbstractControllerCache, alg)::Bool
+
+This function decides whether the current time step should be accepted or rejected.
+A return value of `false` corresponds to a rejection.
+"""
+accept_step_controller
+
+"""
+    stepsize_controller!(integrator, alg)
+    stepsize_controller!(integrator, controller_cache::AbstractControllerCache, alg)
+
+Update the cache to compute the new order of the time marching algorithm and prepare
+for an update of the time step. The update can be either due to a rejection or acceptance
+of the current time step.
+"""
+stepsize_controller!
+
+"""
+    step_accept_controller!(integrator, alg, q)
+    step_accept_controller!(integrator, controller_cache::AbstractControllerCache, alg, q)
+
+This function gets called in case of an accepted time step right after [`stepsize_controller!`](@ref).
+It returns the proposed new time step length. Please note that the time step length might not be
+applied as is and subject to further modification to e.g. match the next time stop.
+
+!!! warning
+    The parameter `q` will be removed the next release and will be passed through the
+    controller cache if needed.
+"""
+step_accept_controller!
+
+"""
+    step_reject_controller!(integrator, alg)
+    step_reject_controller!(integrator, controller_cache::AbstractControllerCache, alg)
+
+This function gets called in case of a rejected time step right after [`stepsize_controller!`](@ref).
+It directly sets the time step length (i.e. `integrator.dt`).
+"""
+step_reject_controller!
+
+
+# The legacy controllers do not have this concept.
+setup_controller_cache(alg, controller::AbstractLegacyController) = controller
+
+# checks whether the controller should accept a step based on the error estimate
+@inline function accept_step_controller(integrator, controller_or_cache::Union{<:AbstractLegacyController, <:AbstractControllerCache})
+    return integrator.EEst <= 1
+end
 
 @inline function stepsize_controller!(integrator, alg)
+    # TODO replace this when done - right now this holds the controller cache!
     stepsize_controller!(integrator, integrator.opts.controller, alg)
+    # stepsize_controller!(integrator, integrator.controller_cache, alg)
 end
 
 @inline function step_accept_controller!(integrator, alg, q)
+    # TODO replace this when done - right now this holds the controller cache!
     step_accept_controller!(integrator, integrator.opts.controller, alg, q)
+    # step_accept_controller!(integrator, integrator.controller_cache, alg, q)
 end
 
 @inline function step_reject_controller!(integrator, alg)
+     # TODO replace this when done - right now this holds the controller cache!
     step_reject_controller!(integrator, integrator.opts.controller, alg)
+    # step_reject_controller!(integrator, integrator.controller_cache, alg)
     cache = integrator.cache
     if hasfield(typeof(cache), :nlsolve)
         nlsolve = cache.nlsolve
@@ -23,6 +90,21 @@ reset_alg_dependent_opts!(controller::AbstractController, alg1, alg2) = nothing
 
 SciMLBase.reinit!(integrator::ODEIntegrator, controller::AbstractController) = nothing
 
+function post_newton_controller!(integrator, alg)
+    integrator.dt = integrator.dt / integrator.opts.failfactor
+    nothing
+end
+
+# This is a helper for Nordsieck and BDF methods, which come with an integrated controller.
+struct DummyController <: AbstractController
+end
+
+setup_controller_cache(alg, controller::DummyController) = controller
+
+@inline function accept_step_controller(integrator, controller::DummyController)
+    return integrator.EEst <= 1
+end
+
 # Standard integral (I) step size controller
 """
     IController()
@@ -54,10 +136,7 @@ the predicted step size.
     Solving Ordinary Differential Equations I Nonstiff Problems
     [DOI: 10.1007/978-3-540-78862-1](https://doi.org/10.1007/978-3-540-78862-1)
 """
-struct IController <: AbstractController
-end
-
-struct DummyController <: AbstractController
+struct IController <: AbstractLegacyController
 end
 
 @inline function stepsize_controller!(integrator, controller::IController, alg)
@@ -90,6 +169,77 @@ function step_reject_controller!(integrator, controller::IController, alg)
     integrator.dt = qold
 end
 
+struct NewIController{T} <: AbstractController
+    qmin::T
+    qmax::T
+    gamma::T
+    qsteady_min::T
+    qsteady_max::T
+end
+
+function NewIController(alg; kwargs...)
+    NewIController(Float64, alg; kwargs...)
+end
+
+function NewIController(QT, alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing)
+    NewIController{QT}(
+        qmin === nothing ? qmin_default(alg) : qmin,
+        qmax === nothing ? qmax_default(alg) : qmax,
+        gamma === nothing ? gamma_default(alg) : gamma,
+        qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min,
+        qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max,
+    )
+end
+
+mutable struct IControllerCache{C, T} <: AbstractControllerCache
+    controller::C
+    q::T
+    dtreject::T
+    # I believe this should go here or in the algorithm cache, but not in the integrator itself.
+    # EEst::T
+end
+
+function setup_controller_cache(alg, controller::NewIController{T}) where T
+    IControllerCache(
+        controller,
+        T(1),
+        T(1 // 10^4), # TODO which value?
+    )
+end
+
+@inline function stepsize_controller!(integrator, cache::IControllerCache, alg)
+    @unpack qmin, qmax, gamma = cache.controller
+    EEst = DiffEqBase.value(integrator.EEst)
+
+    if iszero(EEst)
+        cache.q = inv(qmax)
+    else
+        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
+        qtmp = fastpower(EEst, expo) / gamma
+        @fastmath cache.q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
+        # TODO: Shouldn't this be in `step_accept_controller!` as for the PI controller?
+        cache.dtreject = integrator.qold = DiffEqBase.value(integrator.dt) / cache.q
+    end
+    cache.q
+end
+
+# TODO change signature to remove the q input
+function step_accept_controller!(integrator, cache::IControllerCache, alg, q)
+    @unpack qsteady_min, qsteady_max = cache.controller
+    @assert q ≈ cache.q "Controller cache went out of sync with time stepping logic."
+
+    if qsteady_min <= q <= qsteady_max
+        cache.q = q = one(q)
+    end
+    integrator.dt / q # new dt
+end
+
+function step_reject_controller!(integrator, cache::IControllerCache, alg)
+    @assert cache.dtreject ≈ integrator.qold "Controller cache went out of sync with time stepping logic."
+    integrator.dt = cache.dtreject # TODO this does not look right.
+end
+
+
 # PI step size controller
 """
     PIController(beta1, beta2)
@@ -127,7 +277,7 @@ the predicted step size.
     Solving Ordinary Differential Equations I Nonstiff Problems
     [DOI: 10.1007/978-3-540-78862-1](https://doi.org/10.1007/978-3-540-78862-1)
 """
-mutable struct PIController{QT} <: AbstractController
+mutable struct PIController{QT} <: AbstractLegacyController
     beta1::QT
     beta2::QT
 end
@@ -175,6 +325,97 @@ function reset_alg_dependent_opts!(controller::PIController, alg1, alg2)
     end
 end
 
+struct NewPIController{T} <: AbstractController
+    beta1::T
+    beta2::T
+    qmin::T
+    qmax::T
+    gamma::T
+    qsteady_min::T
+    qsteady_max::T
+    qoldinit::T
+end
+
+function NewPIController(alg; kwargs...)
+    NewPIController(Float64, alg; kwargs...)
+end
+
+function NewPIController(QT, alg; beta1 = nothing, beta2 = nothing, qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing, qoldinit = nothing)
+    beta2 = beta2 === nothing ? beta2_default(alg) : beta2
+    beta1 = beta1 === nothing ? beta1_default(alg, beta2) : beta1
+    qoldinit = qoldinit === nothing ? 1 // 10^4 : qoldinit
+    NewPIController{QT}(
+        beta1,
+        beta2,
+        qmin === nothing ? qmin_default(alg) : qmin,
+        qmax === nothing ? qmax_default(alg) : qmax,
+        gamma === nothing ? gamma_default(alg) : gamma,
+        qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min,
+        qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max,
+        qoldinit,
+    )
+end
+
+mutable struct PIControllerCache{T} <: AbstractControllerCache
+    controller::NewPIController{T}
+    # Propsoed scaling factor for the time step length
+    q::T
+    # Cached εₙ₊₁^β₁
+    q11::T
+    # Previous EEst
+    errold::T
+end
+
+
+function setup_controller_cache(alg, controller::NewPIController{T}) where T
+    PIControllerCache(
+        controller,
+        T(1),
+        T(1),
+        T(1 // 10^4),
+    )
+end
+
+@inline function stepsize_controller!(integrator, cache::PIControllerCache, alg)
+    @unpack errold, controller = cache
+    @unpack qmin, qmax, gamma = controller
+    @unpack beta1, beta2 = controller
+    EEst = DiffEqBase.value(integrator.EEst)
+
+    if iszero(EEst)
+        q = qmax
+    else
+        # Legacy code
+        q11 = fastpower(EEst, beta1)
+        q = q11 / fastpower(errold, beta2)
+        cache.q11 = q11
+        integrator.q11 = q11 # TODO remove
+        @fastmath q = clamp(q / gamma, inv(qmax), inv(qmin))
+    end
+    cache.q = q
+    q
+end
+
+function step_accept_controller!(integrator, cache::PIControllerCache, alg, q)
+    @assert q ≈ cache.q "Controller cache went out of sync with time stepping logic (q=$q | cache.q=$(cache.q))."
+    @unpack controller = cache
+    @unpack qsteady_min, qsteady_max, qoldinit = controller
+    EEst = DiffEqBase.value(integrator.EEst)
+
+    if qsteady_min <= q <= qsteady_max
+        q = one(q)
+    end
+    cache.errold = max(EEst, qoldinit)
+    integrator.qold = cache.errold # TODO remove
+    return integrator.dt / q # new dt
+end
+
+function step_reject_controller!(integrator, cache::PIControllerCache, alg)
+    @unpack controller, q11 = cache
+    @unpack qmin, gamma = controller
+    integrator.dt /= min(inv(qmin), q11 / gamma)
+end
+
 # PID step size controller
 """
     PIDController(beta1, beta2, beta3=zero(beta1);
@@ -241,7 +482,7 @@ Some standard controller parameters suggested in the literature are
     Compressible Computational Fluid Dynamics    # is bigger than this parameter
     [arXiv:2104.06836](https://arxiv.org/abs/2104.06836)    # limiter of the dt factor (before clipping)
 """
-struct PIDController{QT, Limiter} <: AbstractController
+struct PIDController{QT, Limiter} <: AbstractLegacyController
     beta::MVector{3, QT} # controller coefficients
     err::MVector{3, QT} # history of the error estimates
     accept_safety::QT   # accept a step if the predicted change of the step size
@@ -289,7 +530,7 @@ end
     #   EEst = eps(typeof(EEst))
     # end
     # ```
-    EEst = ifelse(EEst > EEst_min, EEst, EEst_min)
+    EEst = max(EEst, EEst_min)
 
     controller.err[1] = inv(EEst)
     err1, err2, err3 = controller.err
@@ -311,11 +552,6 @@ end
     return dt_factor
 end
 
-# checks whether the controller should accept a step based on the error estimate
-@inline function accept_step_controller(integrator, controller::AbstractController)
-    return integrator.EEst <= 1
-end
-
 @inline function accept_step_controller(integrator, controller::PIDController)
     return integrator.qold >= controller.accept_safety
 end
@@ -337,6 +573,134 @@ function step_reject_controller!(integrator, controller::PIDController, alg)
     integrator.dt *= integrator.qold
 end
 
+
+
+struct NewPIDController{T, Limiter} <: AbstractController
+    beta::SVector{3, T} # controller coefficients
+    accept_safety::T   # accept a step if the predicted change of the step size
+    # is bigger than this parameter
+    limiter::Limiter    # limiter of the dt factor (before clipping)
+    qmin::T
+    qmax::T
+    qsteady_min::T
+    qsteady_max::T
+end
+
+function NewPIDController(alg; kwargs...)
+    NewPIDController(Float64, alg; kwargs...)
+end
+
+function NewPIDController(QT, alg; beta = nothing, accept_safety = 0.81, limiter = default_dt_factor_limiter, qmin = nothing, qmax = nothing, qsteady_min = nothing, qsteady_max = nothing)
+    if beta === nothing 
+        beta2 = QT(beta2_default(alg))
+        beta1 = QT(beta1_default(alg, beta2))
+        beta3 = QT(zero(beta1))
+    else
+        beta1, beta2, beta3 = beta
+    end
+    beta = SVector(map(float, promote(beta1, beta2, beta3))...)
+    NewPIDController{QT, typeof(limiter)}(
+        beta,
+        QT(accept_safety),
+        limiter,
+        QT(qmin === nothing ? qmin_default(alg) : qmin),
+        QT(qmax === nothing ? qmax_default(alg) : qmax),
+        QT(qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min),
+        QT(qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max),
+    )
+end
+
+function Base.show(io::IO, controller::NewPIDController)
+    print(io, "NewPIDController(beta=", controller.beta,
+        ", accept_safety=", controller.accept_safety,
+        ", limiter=", controller.limiter,
+        ")")
+end
+
+mutable struct PIDControllerCache{T, Limiter} <: AbstractControllerCache
+    controller::NewPIDController{T, Limiter}
+    err::MVector{3, T} # history of the error estimates
+    dt_factor::T
+end
+
+function setup_controller_cache(alg, controller::NewPIDController{QT}) where QT
+    err = MVector{3, QT}(true, true, true)
+    PIDControllerCache(
+        controller,
+        err,
+        QT(1 // 10^4),
+    )
+end
+
+@inline function stepsize_controller!(integrator, cache::PIDControllerCache, alg)
+    @unpack controller = cache
+    beta1, beta2, beta3 = controller.beta
+    @unpack qmax = controller
+
+    EEst = DiffEqBase.value(integrator.EEst)
+
+    # If the error estimate is zero, we can increase the step size as much as
+    # desired. This additional check fixes problems of the code below when the
+    # error estimates become zero
+    # -> err1, err2, err3 become Inf
+    # -> err1^positive_number * err2^negative_number becomes NaN
+    # -> dt becomes NaN
+    #
+    # `EEst_min` is smaller than PETSC_SMALL used in the equivalent logic in PETSc.
+    # For example, `eps(Float64) ≈ 2.2e-16` but `PETSC_SMALL ≈ 1.0e-10` for `double`.
+    EEst_min = eps(typeof(EEst))
+    # The code below is a bit more robust than
+    # ```
+    # if iszero(EEst)
+    #   EEst = eps(typeof(EEst))
+    # end
+    # ```
+    EEst = max(EEst, EEst_min)
+
+    cache.err[1] = inv(EEst)
+    err1, err2, err3 = cache.err
+
+    k = min(alg_order(alg), alg_adaptive_order(alg)) + 1
+    dt_factor = err1^(beta1 / k) * err2^(beta2 / k) * err3^(beta3 / k)
+    if isnan(dt_factor)
+        @warn "unlimited dt_factor" dt_factor err1 err2 err3 beta1 beta2 beta3 k
+    end
+    dt_factor = controller.limiter(dt_factor)
+
+    # Note: No additional limiting of the form
+    #   dt_factor = max(qmin, min(qmax, dt_factor))
+    # is necessary since the `limiter` should take care of that. The default limiter
+    # ensures
+    #   0.21 ≈ limiter(0) <= dt_factor <= limiter(Inf) ≈ 2.57
+    # See Söderlind, Wang (2006), Section 6.
+    cache.dt_factor = dt_factor
+    return dt_factor
+end
+
+@inline function accept_step_controller(integrator, cache::PIDControllerCache)
+    return cache.dt_factor >= cache.controller.accept_safety
+end
+
+function step_accept_controller!(integrator, cache::PIDControllerCache, alg, dt_factor)
+    @assert dt_factor ≈ cache.dt_factor "Controller cache went out of sync with time stepping logic."
+    @unpack controller = cache
+    @unpack qsteady_min, qsteady_max = controller
+
+    if qsteady_min <= inv(dt_factor) <= qsteady_max
+        dt_factor = one(dt_factor)
+        # cache.q = ...?
+    end
+    @inbounds begin
+        cache.err[3] = cache.err[2]
+        cache.err[2] = cache.err[1]
+    end
+    return integrator.dt * dt_factor # new dt
+end
+
+function step_reject_controller!(integrator, cache::PIDControllerCache, alg)
+    integrator.dt *= cache.qold
+end
+
 # Gustafsson predictive step size controller
 """
     PredictiveController()
@@ -395,11 +759,6 @@ end
 struct PredictiveController <: AbstractController
 end
 
-function post_newton_controller!(integrator, alg)
-    integrator.dt = integrator.dt / integrator.opts.failfactor
-    nothing
-end
-
 @inline function stepsize_controller!(integrator, controller::PredictiveController, alg)
     @unpack qmin, qmax, gamma = integrator.opts
     EEst = DiffEqBase.value(integrator.EEst)
@@ -452,3 +811,136 @@ function step_reject_controller!(integrator, controller::PredictiveController, a
     @unpack dt, success_iter, qold = integrator
     integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
 end
+
+
+struct NewPredictiveController{T} <: AbstractController
+    qmin::T
+    qmax::T
+    gamma::T
+    qsteady_min::T
+    qsteady_max::T
+end
+
+function NewPredictiveController(alg; kwargs...)
+    PIController(Float64, alg; kwargs...)
+end
+
+function NewPredictiveController(QT, alg; qmin = nothing, qmax = nothing, gamma = nothing, qsteady_min = nothing, qsteady_max = nothing)
+    NewPredictiveController{QT}(
+        qmin === nothing ? qmin_default(alg) : qmin,
+        qmax === nothing ? qmax_default(alg) : qmax,
+        gamma === nothing ? gamma_default(alg) : gamma,
+        qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min,
+        qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max,
+    )
+end
+
+mutable struct PredictiveControllerCache{T} <: AbstractControllerCache
+    controller::NewPredictiveController{T}
+    dtacc::T
+    erracc::T
+    qold::T
+    q::T
+end
+
+function setup_controller_cache(alg, controller::NewPredictiveController{T}) where T
+    PredictiveControllerCache{T}(
+        controller,
+        T(1),
+        T(1),
+        T(1),
+        T(1)
+    )
+end
+
+@inline function stepsize_controller!(integrator, cache::PredictiveControllerCache, alg)
+    @unpack qmin, qmax, gamma = integrator.opts
+    EEst = DiffEqBase.value(integrator.EEst)
+    if iszero(EEst)
+        q = inv(qmax)
+    else
+        if fac_default_gamma(alg)
+            fac = gamma
+        else
+            if isfirk(alg)
+                @unpack iter = integrator.cache
+                @unpack maxiters = alg
+            else
+                @unpack iter, maxiters = integrator.cache.nlsolver
+            end
+            fac = min(gamma, (1 + 2 * maxiters) * gamma / (iter + 2 * maxiters))
+        end
+        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
+        qtmp = fastpower(EEst, expo) / fac
+        @fastmath q = DiffEqBase.value(max(inv(qmax), min(inv(qmin), qtmp)))
+        cache.qold = q
+    end
+    cache.q = q
+    q
+end
+
+function step_accept_controller!(integrator, cache::PredictiveControllerCache, alg, q)
+    @assert q ≈ cache.q "Controller cache went out of sync with time stepping logic."
+    @unpack dtacc, erracc, controller = cache
+    @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = controller
+
+    EEst = DiffEqBase.value(integrator.EEst)
+
+    if integrator.success_iter > 0
+        expo = 1 / (get_current_adaptive_order(alg, integrator.cache) + 1)
+        qgus = (dtacc / integrator.dt) *
+               fastpower((EEst^2) / erracc, expo)
+        qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
+        qacc = max(q, qgus)
+    else
+        qacc = q
+    end
+    if qsteady_min <= qacc <= qsteady_max
+        qacc = one(qacc)
+    end
+    cache.dtacc = integrator.dt
+    cache.erracc = max(1e-2, EEst)
+
+    return integrator.dt / qacc
+end
+
+function step_reject_controller!(integrator, cache::PredictiveControllerCache, alg)
+    @unpack dt, success_iter = integrator
+    @unpack qold = cache
+    integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
+end
+
+# For a composite algorithm the default strategy is to switch forth and back between the controllers of the individual algorithms
+struct CompositeController{T} <: AbstractController
+    controllers::T
+end
+
+struct CompositeControllerCache{T} <: AbstractControllerCache
+    caches::T
+end
+
+function setup_controller_cache(alg::CompositeAlgorithm, cc::CompositeController)
+    CompositeControllerCache(
+        map((alg, controller)->setup_controller_cache(alg, controller), alg.algs, cc.controllers)
+    )
+end
+
+@inline function accept_step_controller(integrator, cache::CompositeControllerCache, alg::CompositeAlgorithm)
+    current_idx = integrator.cache.current
+    accept_step_controller(integrator, cache.caches[current_idx], alg.algs[current_idx])
+end
+
+@inline function stepsize_controller!(integrator, cache::CompositeControllerCache, alg::CompositeAlgorithm)
+    current_idx = integrator.cache.current
+    stepsize_controller!(integrator, cache.caches[current_idx], alg.algs[current_idx])
+end
+
+@inline function step_accept_controller!(integrator, cache::CompositeControllerCache, alg::CompositeAlgorithm, q)
+    current_idx = integrator.cache.current
+    step_accept_controller!(integrator, cache.caches[current_idx], alg.algs[current_idx], q)
+end
+
+@inline function step_reject_controller!(integrator, cache::CompositeControllerCache, alg::CompositeAlgorithm)
+    current_idx = integrator.cache.current
+    step_reject_controller!(integrator, cache.caches[current_idx], alg.algs[current_idx])
+end
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/type.jl b/lib/OrdinaryDiffEqCore/src/integrators/type.jl
index 7dc9e1a9c6..804fd275f3 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/type.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/type.jl
@@ -6,12 +6,14 @@ mutable struct DEOptions{absType, relType, QT, tType, Controller, F1, F2, F3, F4
     adaptive::Bool
     abstol::absType
     reltol::relType
+    # TODO vvv remove this block as these are controller and not integrator parameters vvv
     gamma::QT
     qmax::QT
     qmin::QT
     qsteady_max::QT
     qsteady_min::QT
     qoldinit::QT
+    # TODO ^^^ remove this block as these are controller and not integrator parameters ^^^
     failfactor::QT
     dtmax::tType
     dtmin::tType
@@ -84,7 +86,7 @@ For more info see the linked documentation page.
 mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgorithm}, IIP,
     uType, duType, tType, pType, eigenType, EEstT, QT, tdirType,
     ksEltype, SolType, F, CacheType, O, FSALType, EventErrorType,
-    CallbackCacheType, IA, DV} <:
+    CallbackCacheType, IA, DV, CC} <:
                SciMLBase.AbstractODEIntegrator{algType, IIP, uType, tType}
     sol::SolType
     u::uType
@@ -105,10 +107,13 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori
     tdir::tdirType
     eigen_est::eigenType
     EEst::EEstT
+    # TODO vvv remove these
     qold::QT
     q11::QT
     erracc::QT
     dtacc::tType
+    # TODO ^^^ remove these
+    controller_cache::CC
     success_iter::Int
     iter::Int
     saveiter::Int
diff --git a/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl b/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl
index c3c101a47f..a0159c843f 100644
--- a/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl
+++ b/lib/OrdinaryDiffEqCore/src/perform_step/composite_perform_step.jl
@@ -256,9 +256,6 @@ function choose_algorithm!(integrator, cache::DefaultCache)
             old_cache = cache.cache6
         end
 
-        integrator.opts.controller.beta2 = beta2 = beta2_default(algs[new_current])
-        integrator.opts.controller.beta1 = beta1_default(algs[new_current], beta2)
-
         reset_alg_dependent_opts!(integrator, algs[old_current], algs[new_current])
         transfer_cache!(integrator, old_cache, new_cache)
     end
@@ -275,12 +272,6 @@ function reset_alg_dependent_opts!(integrator, alg1, alg2)
     if integrator.opts.adaptive == isadaptive(alg1)
         integrator.opts.adaptive = isadaptive(alg2)
     end
-    if integrator.opts.qmin == qmin_default(alg1)
-        integrator.opts.qmin = qmin_default(alg2)
-    end
-    if integrator.opts.qmax == qmax_default(alg1)
-        integrator.opts.qmax == qmax_default(alg2)
-    end
     reset_alg_dependent_opts!(integrator.opts.controller, alg1, alg2)
     nothing
 end
@@ -289,3 +280,4 @@ end
 # Example: send the history variables from one multistep method to another
 
 transfer_cache!(integrator, alg1, alg2) = nothing
+reset_alg_dependent_opts!(::Union{AbstractController, AbstractControllerCache}, alg1, alg2) = nothing
diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl
index 6ddd866f87..57956b82f8 100644
--- a/lib/OrdinaryDiffEqCore/src/solve.jl
+++ b/lib/OrdinaryDiffEqCore/src/solve.jl
@@ -36,17 +36,17 @@ function SciMLBase.__init(
         dtmax = eltype(prob.tspan)((prob.tspan[end] - prob.tspan[1])),
         force_dtmin = false,
         adaptive = anyadaptive(alg),
-        gamma = gamma_default(alg),
         abstol = nothing,
         reltol = nothing,
-        qmin = qmin_default(alg),
-        qmax = qmax_default(alg),
-        qsteady_min = qsteady_min_default(alg),
-        qsteady_max = qsteady_max_default(alg),
+        gamma = nothing,
+        qmin = nothing,
+        qmax = nothing,
+        qsteady_min = nothing,
+        qsteady_max = nothing,
         beta1 = nothing,
         beta2 = nothing,
-        qoldinit = anyadaptive(alg) ? 1 // 10^4 : 0,
-        controller = nothing,
+        qoldinit = nothing,
+        controller = any((gamma, qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit) .!== nothing) ? nothing : default_controller_v7(typeof(one(eltype(prob.tspan))), alg), # We have to reconstruct the old controller before breaking release.,
         fullnormalize = true,
         failfactor = 2,
         maxiters = anyadaptive(alg) ? 1000000 : typemax(Int),
@@ -417,23 +417,60 @@ function SciMLBase.__init(
         throw(ArgumentError("Setting both the legacy PID parameters `beta1, beta2 = $((beta1, beta2))` and the `controller = $controller` is not allowed."))
     end
 
+    # Deprecation warnings for users to break down which parameters they accidentally set.
     if (beta1 !== nothing || beta2 !== nothing)
         message = "Providing the legacy PID parameters `beta1, beta2` is deprecated. Use the keyword argument `controller` instead."
         Base.depwarn(message, :init)
         Base.depwarn(message, :solve)
     end
+    if (qmin !== nothing || qmax !== nothing)
+        message = "Providing the legacy PID parameters `qmin, qmax` is deprecated. Use the keyword argument `controller` instead."
+        Base.depwarn(message, :init)
+        Base.depwarn(message, :solve)
+    end
+    if (qsteady_min !== nothing || qsteady_max !== nothing)
+        message = "Providing the legacy PID parameters `qsteady_min, qsteady_max` is deprecated. Use the keyword argument `controller` instead."
+        Base.depwarn(message, :init)
+        Base.depwarn(message, :solve)
+    end
+    if (qoldinit !== nothing)
+        message = "Providing the legacy PID parameters `qoldinit` is deprecated. Use the keyword argument `controller` instead."
+        Base.depwarn(message, :init)
+        Base.depwarn(message, :solve)
+    end
 
-    if controller === nothing
-        controller = default_controller(_alg, cache, qoldinit, beta1, beta2)
+    # The following code provides an upgrade path for users by preserving the old behavior.
+    legacy_controller_parameters = (gamma, qmin, qmax, qsteady_min, qsteady_max, beta1, beta2, qoldinit)
+    if controller === nothing # We have to reconstruct the old controller before breaking release.
+        if any(legacy_controller_parameters .== nothing)
+            gamma = gamma === nothing ? gamma_default(alg) : gamma
+            qmin = qmin === nothing ? qmin_default(alg) : qmin
+            qmax = qmax === nothing ? qmax_default(alg) : qmax
+            qsteady_min = qsteady_min === nothing ? qsteady_min_default(alg) : qsteady_min
+            qsteady_max = qsteady_max === nothing ? qsteady_max_default(alg) : qsteady_max
+            qoldinit = qoldinit === nothing ? (anyadaptive(alg) ? 1 // 10^4 : 0) : qoldinit
+            controller = legacy_default_controller(_alg, cache, qoldinit, beta1, beta2)
+        end
+    else # Controller has been passed
+        gamma = hasfield(typeof(controller), :gamma) ? controller.gamma : gamma_default(alg)
+        qmin = hasfield(typeof(controller), :qmin) ? controller.qmin : qmin_default(alg)
+        qmax = hasfield(typeof(controller), :qmax) ? controller.qmax : qmax_default(alg)
+        qsteady_min = hasfield(typeof(controller), :qsteady_min) ? controller.qsteady_min : qsteady_min_default(alg)
+        qsteady_max = hasfield(typeof(controller), :qsteady_max) ? controller.qsteady_max : qsteady_max_default(alg)
+        qoldinit = hasfield(typeof(controller), :qoldinit) ? controller.qoldinit : (anyadaptive(alg) ? 1 // 10^4 : 0)
     end
 
+    controller_cache = setup_controller_cache(_alg, controller)
+
     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(controller),
+        QT, tType, 
+        # typeof(controller),
+        typeof(controller_cache),
         typeof(internalnorm), typeof(internalopnorm),
         typeof(save_end_user),
         typeof(callbacks_internal),
@@ -446,14 +483,20 @@ function SciMLBase.__init(
         typeof(saveat), typeof(d_discontinuities)}(maxiters, save_everystep,
         adaptive, abstol_internal,
         reltol_internal,
-        QT(gamma), QT(qmax),
+        # TODO vvv remove this block as these are controller and not integrator parameters vvv
+        QT(gamma),
+        QT(qmax),
         QT(qmin),
         QT(qsteady_max),
         QT(qsteady_min),
         QT(qoldinit),
+        # TODO ^^^remove this block as these are controller and not integrator parameters ^^^
         QT(failfactor),
         tType(dtmax), tType(dtmin),
-        controller,
+        # TODO vvv remove this vvv
+        # controller,
+        controller_cache,
+        # TODO ^^^ remove this ^^^
         internalnorm,
         internalopnorm,
         save_idxs, tstops_internal,
@@ -530,13 +573,18 @@ function SciMLBase.__init(
         FType, cacheType,
         typeof(opts), typeof(fsalfirst),
         typeof(last_event_error), typeof(callback_cache),
-        typeof(initializealg), typeof(differential_vars)}(
+        typeof(initializealg), typeof(differential_vars),
+        typeof(controller_cache)}(
         sol, u, du, k, t, tType(dt), f, p,
         uprev, uprev2, duprev, tprev,
         _alg, dtcache, dtchangeable,
         dtpropose, tdir, eigen_est, EEst,
+        # TODO vvv remove these
         QT(qoldinit), q11,
-        erracc, dtacc, success_iter,
+        erracc, dtacc,
+        # TODO ^^^ remove these
+        controller_cache,
+        success_iter,
         iter, saveiter, saveiter_dense, cache,
         callback_cache,
         kshortsize, force_stepfail,
diff --git a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl
index 194b30e166..ba684b58a1 100644
--- a/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl
+++ b/lib/OrdinaryDiffEqExtrapolation/src/OrdinaryDiffEqExtrapolation.jl
@@ -2,23 +2,23 @@ module OrdinaryDiffEqExtrapolation
 
 import OrdinaryDiffEqCore: alg_order, alg_maximum_order, get_current_adaptive_order,
                            get_current_alg_order, calculate_residuals!,
-                           accept_step_controller,
-                           default_controller, beta2_default, beta1_default, gamma_default,
+                           accept_step_controller, default_controller_v7,
+                           legacy_default_controller, beta2_default, beta1_default, gamma_default,
                            initialize!, perform_step!, @unpack, @cache, unwrap_alg,
-                           isthreaded,
+                           isthreaded, isadaptive, PIController, AbstractControllerCache,
                            step_accept_controller!, calculate_residuals,
                            OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
-                           reset_alg_dependent_opts!, AbstractController,
+                           reset_alg_dependent_opts!, AbstractController, AbstractLegacyController,
                            step_accept_controller!, step_reject_controller!,
                            OrdinaryDiffEqAdaptiveAlgorithm,
                            OrdinaryDiffEqAdaptiveImplicitAlgorithm,
                            alg_cache, CompiledFloats, @threaded, stepsize_controller!,
-                           DEFAULT_PRECS, full_cache, qmin_default,
+                           DEFAULT_PRECS, full_cache, qmin_default, qmax_default,
                            constvalue, PolyesterThreads, Sequential, BaseThreads,
                            _digest_beta1_beta2, timedepentdtmin, _unwrap_val,
                            _reshape, _vec, get_fsalfirstlast, generic_solver_docstring,
                            differentiation_rk_docstring, _bool_to_ADType,
-                           _process_AD_choice, LinearAliasSpecifier
+                           _process_AD_choice, LinearAliasSpecifier, setup_controller_cache
 using FastBroadcast, Polyester, MuladdMacro, RecursiveArrayTools, LinearSolve
 import OrdinaryDiffEqCore
 import FastPower
diff --git a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl
index ec2215e511..425063990c 100644
--- a/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl
+++ b/lib/OrdinaryDiffEqExtrapolation/src/alg_utils.jl
@@ -32,7 +32,7 @@ alg_maximum_order(alg::ImplicitHairerWannerExtrapolation) = 2(alg.max_order + 1)
 alg_maximum_order(alg::ImplicitEulerExtrapolation) = 2(alg.max_order + 1)
 alg_maximum_order(alg::ImplicitEulerBarycentricExtrapolation) = alg.max_order
 
-function default_controller(
+function legacy_default_controller(
         alg::Union{ExtrapolationMidpointDeuflhard,
             ImplicitDeuflhardExtrapolation,
             ExtrapolationMidpointHairerWanner,
@@ -46,6 +46,31 @@ function default_controller(
     return ExtrapolationController(beta1)
 end
 
+function default_controller_v7(QT,
+    alg::Union{
+        ExtrapolationMidpointDeuflhard,
+        ImplicitDeuflhardExtrapolation,
+        ExtrapolationMidpointHairerWanner,
+        ImplicitHairerWannerExtrapolation,
+        ImplicitEulerExtrapolation,
+        ImplicitEulerBarycentricExtrapolation})
+    return NewExtrapolationController(QT, alg)
+end
+
+# FIXME AitkenNeville is missing integration with the extrapolation controller and picks up the PI controller instead.
+function legacy_default_controller(alg::AitkenNeville,
+    cache,
+    qoldinit, _beta1 = nothing, _beta2 = nothing)
+    QT = typeof(qoldinit)
+    beta1, beta2 = _digest_beta1_beta2(alg, cache, Val(QT), _beta1, _beta2)
+    return PIController(beta1, beta2)
+end
+function default_controller_v7(QT, alg::AitkenNeville)
+    beta2 = QT(beta2_default(alg))
+    beta1 = QT(beta1_default(alg, beta2))
+    return PIController(beta1, beta2)
+end
+
 beta2_default(alg::ExtrapolationMidpointDeuflhard) = 0 // 1
 beta2_default(alg::ImplicitDeuflhardExtrapolation) = 0 // 1
 beta2_default(alg::ExtrapolationMidpointHairerWanner) = 0 // 1
diff --git a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl
index a62d61cb77..2f2b25032b 100644
--- a/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl
+++ b/lib/OrdinaryDiffEqExtrapolation/src/controllers.jl
@@ -1,5 +1,5 @@
 # Extrapolation methods
-mutable struct ExtrapolationController{QT} <: AbstractController
+mutable struct ExtrapolationController{QT} <: AbstractLegacyController
     beta1::QT
 end
 
@@ -9,7 +9,40 @@ function reset_alg_dependent_opts!(controller::ExtrapolationController, alg1, al
     end
 end
 
+struct NewExtrapolationController{QT} <: AbstractController
+    qmin::QT
+    qmax::QT
+end
+
+function NewExtrapolationController(QT, alg; qmin = nothing, qmax = nothing, gamma = nothing)
+    NewExtrapolationController(
+        QT(qmin === nothing ? qmin_default(alg) : qmin),
+        QT(qmax === nothing ? qmax_default(alg) : qmax),
+    )
+end
+
+mutable struct ExtrapolationControllerCache{QT} <: AbstractControllerCache
+    controller::NewExtrapolationController{QT}
+    beta1::QT
+    gamma::QT
+end
+
+function setup_controller_cache(alg, controller::NewExtrapolationController{T}) where T
+    ExtrapolationControllerCache{T}(
+        controller,
+        T(1),
+        T(1)
+    )
+end
+
+# TODO replace this when done - right now this holds the controller cache!
+stepsize_controller_internal!(integrator, alg) = stepsize_controller_internal!(integrator, integrator.opts.controller, alg)
+stepsize_predictor!(integrator, alg, n_new) = stepsize_predictor!(integrator, integrator.opts.controller, alg, n_new)
+# stepsize_controller_internal!(integrator, alg) = stepsize_controller_internal!(integrator, integrator.controller_cache, alg)
+# stepsize_predictor!(integrator, alg, n_new) = stepsize_predictor!(integrator, integrator.controller_cache, alg, n_new)
+
 @inline function stepsize_controller!(integrator,
+        controller::ExtrapolationController,
         alg::Union{ExtrapolationMidpointDeuflhard,
             ImplicitDeuflhardExtrapolation})
     # Dummy function
@@ -19,7 +52,19 @@ end
     zero(typeof(integrator.opts.qmax))
 end
 
+@inline function stepsize_controller!(integrator,
+        cache::ExtrapolationControllerCache,
+        alg::Union{ExtrapolationMidpointDeuflhard,
+            ImplicitDeuflhardExtrapolation})
+    # Dummy function
+    # ExtrapolationMidpointDeuflhard's stepsize scaling is stored in the cache;
+    # it is computed by  stepsize_controller_internal! (in perform_step!) resp. stepsize_predictor!
+    # (in step_accept_controller! and step_reject_controller!)
+    zero(typeof(cache.controller.qmax))
+end
+
 function stepsize_controller_internal!(integrator,
+        controller::ExtrapolationController,
         alg::Union{ExtrapolationMidpointDeuflhard,
             ImplicitDeuflhardExtrapolation})
     # Standard step size controller
@@ -41,7 +86,32 @@ function stepsize_controller_internal!(integrator,
     integrator.cache.Q[integrator.cache.n_curr - alg.min_order + 1] = q
 end
 
+function stepsize_controller_internal!(integrator,
+        cache::ExtrapolationControllerCache,
+        alg::Union{ExtrapolationMidpointDeuflhard,
+            ImplicitDeuflhardExtrapolation})
+    # Standard step size controller
+    # Compute and save the stepsize scaling based on the latest error estimate of the current order
+    @unpack controller = cache
+
+    if iszero(integrator.EEst)
+        q = inv(controller.qmax)
+    else
+        # Update gamma and beta1
+        cache.beta1 = typeof(cache.beta1)(1 // (2integrator.cache.n_curr + 1))
+        cache.gamma = FastPower.fastpower(typeof(cache.gamma)(1 // 4),
+            cache.beta1)
+        # Compute new stepsize scaling
+        qtmp = FastPower.fastpower(integrator.EEst, cache.beta1) /
+               cache.gamma
+        @fastmath q = max(inv(controller.qmax), min(inv(controller.qmin), qtmp))
+    end
+    integrator.cache.Q[integrator.cache.n_curr - alg.min_order + 1] = q
+end
+
+
 function stepsize_predictor!(integrator,
+        controller::ExtrapolationController,
         alg::Union{ExtrapolationMidpointDeuflhard,
             ImplicitDeuflhardExtrapolation}, n_new::Int)
     # Compute and save the stepsize scaling for order n_new based on the latest error estimate of the current order.
@@ -69,6 +139,35 @@ function stepsize_predictor!(integrator,
     integrator.cache.Q[n_new - alg.min_order + 1] = q
 end
 
+function stepsize_predictor!(integrator,
+        cache::ExtrapolationControllerCache,
+        alg::Union{ExtrapolationMidpointDeuflhard,
+            ImplicitDeuflhardExtrapolation}, n_new::Int)
+    # Compute and save the stepsize scaling for order n_new based on the latest error estimate of the current order.
+    @unpack controller = cache
+
+    if iszero(integrator.EEst)
+        q = inv(controller.qmax)
+    else
+        # Initialize
+        @unpack t, EEst = integrator
+        @unpack stage_number = integrator.cache
+        tol = integrator.opts.internalnorm(integrator.opts.reltol, t) # Deuflhard's approach relies on EEstD ≈ ||relTol||
+        s_curr = stage_number[integrator.cache.n_curr - alg.min_order + 1]
+        s_new = stage_number[n_new - alg.min_order + 1]
+        # Update gamma and beta1
+        cache.beta1 = typeof(cache.beta1)(1 // (2integrator.cache.n_curr + 1))
+        cache.gamma = FastPower.fastpower(typeof(cache.gamma)(1 // 4),
+            cache.beta1)
+        # Compute new stepsize scaling
+        qtmp = EEst *
+               FastPower.fastpower(FastPower.fastpower(tol, (1.0 - s_curr / s_new)),
+            cache.beta1) / cache.gamma
+        @fastmath q = max(inv(controller.qmax), min(inv(controller.qmin), qtmp))
+    end
+    integrator.cache.Q[n_new - alg.min_order + 1] = q
+end
+
 function step_accept_controller!(integrator,
         alg::Union{ExtrapolationMidpointDeuflhard,
             ImplicitDeuflhardExtrapolation}, q)
@@ -136,6 +235,7 @@ end
 end
 
 function stepsize_controller_internal!(integrator,
+        controller::ExtrapolationController,
         alg::Union{ExtrapolationMidpointHairerWanner,
             ImplicitHairerWannerExtrapolation,
             ImplicitEulerExtrapolation,
@@ -191,6 +291,63 @@ function stepsize_controller_internal!(integrator,
     end
 end
 
+function stepsize_controller_internal!(integrator,
+        cache::ExtrapolationControllerCache,
+        alg::Union{ExtrapolationMidpointHairerWanner,
+            ImplicitHairerWannerExtrapolation,
+            ImplicitEulerExtrapolation,
+            ImplicitEulerBarycentricExtrapolation})
+    # Standard step size controller
+    # Compute and save the stepsize scaling based on the latest error estimate of the current order
+    @unpack controller = cache
+
+    if alg isa
+    Union{ImplicitEulerExtrapolation, ImplicitEulerBarycentricExtrapolation,
+        ImplicitHairerWannerExtrapolation}
+        if iszero(integrator.EEst)
+            q = inv(controller.qmax)
+        else
+            # Update gamma and beta1
+            if alg isa ImplicitHairerWannerExtrapolation
+                cache.beta1 = typeof(cache.beta1)(1 //
+                                                            (2integrator.cache.n_curr + 1))
+            elseif alg isa ImplicitEulerExtrapolation
+                cache.beta1 = typeof(cache.beta1)(1 // (integrator.cache.n_curr))
+            else
+                cache.beta1 = typeof(cache.beta1)(1 //
+                                                            (integrator.cache.n_curr - 1))
+            end
+            cache.gamma = FastPower.fastpower(
+                typeof(cache.gamma)(65 //
+                                            100),
+                cache.beta1)
+            # Compute new stepsize scaling
+            qtmp = FastPower.fastpower(integrator.EEst, cache.beta1) /
+                (cache.gamma)
+            @fastmath q = max(inv(controller.qmax),
+                min(inv(controller.qmin), qtmp))
+        end
+        integrator.cache.Q[integrator.cache.n_curr + 1] = q
+    else
+        if iszero(integrator.EEst)
+            q = inv(controller.qmax)
+        else
+            # Update gamma and beta1
+            cache.beta1 = typeof(cache.beta1)(1 // (2integrator.cache.n_curr + 1))
+            cache.gamma = FastPower.fastpower(
+                typeof(cache.gamma)(65 //
+                                            100),
+                cache.beta1)
+            # Compute new stepsize scaling
+            qtmp = FastPower.fastpower(integrator.EEst, cache.beta1) /
+                cache.gamma
+            @fastmath q = max(inv(controller.qmax),
+                min(inv(controller.qmin), qtmp))
+        end
+        integrator.cache.Q[integrator.cache.n_curr + 1] = q
+    end
+end
+
 function step_accept_controller!(integrator,
         alg::Union{ExtrapolationMidpointHairerWanner,
             ImplicitHairerWannerExtrapolation,
diff --git a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
index 7ab3cffccf..e4693dc3b8 100644
--- a/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/OrdinaryDiffEqFIRK.jl
@@ -14,7 +14,8 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!,
                            qmax_default, alg_adaptive_order, DEFAULT_PRECS,
                            stepsize_controller!, step_accept_controller!,
                            step_reject_controller!,
-                           PredictiveController, alg_can_repeat_jac, NewtonAlgorithm,
+                           PredictiveController, PredictiveControllerCache,
+                           NewPredictiveController, alg_can_repeat_jac, NewtonAlgorithm,
                            fac_default_gamma,
                            get_current_adaptive_order, get_fsalfirstlast,
                            isfirk, generic_solver_docstring, _bool_to_ADType,
diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
index f491a9e7d5..87fad2e5d6 100644
--- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
@@ -62,3 +62,70 @@ function step_reject_controller!(
         end
     end
 end
+
+function step_accept_controller!(
+        integrator, ccache::PredictiveControllerCache, alg::AdaptiveRadau, q)
+    @unpack controller = ccache
+    @unpack qmin, qmax, gamma, qsteady_min, qsteady_max = controller
+    @unpack cache = integrator
+    @unpack num_stages, step, iter, hist_iter, index = cache
+
+    EEst = DiffEqBase.value(integrator.EEst)
+
+    if integrator.success_iter > 0
+        expo = 1 / (get_current_adaptive_order(alg, cache) + 1)
+        qgus = (ccache.dtacc / integrator.dt) *
+               fastpow((EEst^2) / ccache.erracc, expo)
+        qgus = max(inv(qmax), min(inv(qmin), qgus / gamma))
+        qacc = max(q, qgus)
+    else
+        qacc = q
+    end
+    if qsteady_min <= qacc <= qsteady_max
+        qacc = one(qacc)
+    end
+    ccache.dtacc = integrator.dt
+    ccache.erracc = max(1e-2, EEst)
+    cache.step = step + 1
+    hist_iter = hist_iter * 0.8 + iter * 0.2
+    cache.hist_iter = hist_iter
+    max_stages = (alg.max_order - 1) ÷ 4 * 2 + 1
+    min_stages = (alg.min_order - 1) ÷ 4 * 2 + 1
+    if (step > 10)
+        if (hist_iter < 2.6 && num_stages < max_stages)
+            cache.num_stages += 2
+            cache.index += 1
+            cache.step = 1
+            cache.hist_iter = iter
+        elseif ((hist_iter > 8 || cache.status == VerySlowConvergence ||
+                 cache.status == Divergence) && num_stages > min_stages)
+            cache.num_stages -= 2
+            cache.index -= 1
+            cache.step = 1
+            cache.hist_iter = iter
+        end
+    end
+    return integrator.dt / qacc
+end
+
+function step_reject_controller!(
+        integrator, ccache::PredictiveControllerCache, alg::AdaptiveRadau)
+    @unpack controller = ccache
+    @unpack dt, success_iter = integrator
+    @unpack cache = integrator
+    @unpack num_stages, step, iter, hist_iter = cache
+    integrator.dt = success_iter == 0 ? 0.1 * dt : dt / ccache.qold
+    cache.step = step + 1
+    hist_iter = hist_iter * 0.8 + iter * 0.2
+    cache.hist_iter = hist_iter
+    min_stages = (alg.min_order - 1) ÷ 4 * 2 + 1
+    if (step > 10)
+        if ((hist_iter > 8 || cache.status == VerySlowConvergence ||
+             cache.status == Divergence) && num_stages > min_stages)
+            cache.num_stages -= 2
+            cache.index -= 1
+            cache.step = 1
+            cache.hist_iter = iter
+        end
+    end
+end
diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl
index 22596ff12b..0028ec16be 100644
--- a/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl
+++ b/lib/OrdinaryDiffEqLowStorageRK/src/OrdinaryDiffEqLowStorageRK.jl
@@ -7,7 +7,8 @@ import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, calculate_residuals!,
                            OrdinaryDiffEqAlgorithm, ispredictive,
                            OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
                            OrdinaryDiffEqAdaptiveAlgorithm, uses_uprev,
-                           default_controller, PIDController,
+                           default_controller_v7,
+                           legacy_default_controller, NewPIDController, PIDController,
                            alg_cache, _vec, _reshape, @cache, isfsal, full_cache,
                            constvalue, _unwrap_val,
                            trivial_limiter!, perform_step!, initialize!,
diff --git a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl
index 1dc1a9e81a..00d39c86be 100644
--- a/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl
+++ b/lib/OrdinaryDiffEqLowStorageRK/src/alg_utils.jl
@@ -84,32 +84,56 @@ uses_uprev(alg::CKLLSRK65_4M_4R, adaptive::Bool) = adaptive
 uses_uprev(alg::CKLLSRK85_4FM_4R, adaptive::Bool) = adaptive
 uses_uprev(alg::CKLLSRK75_4M_5R, adaptive::Bool) = adaptive
 
-function default_controller(alg::RDPK3Sp35, cache, qoldinit, args...)
+function legacy_default_controller(alg::RDPK3Sp35, cache, qoldinit, args...)
     QT = typeof(qoldinit)
     return PIDController(map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04))...)
 end
 
-function default_controller(alg::RDPK3SpFSAL35, cache, qoldinit, args...)
+function legacy_default_controller(alg::RDPK3SpFSAL35, cache, qoldinit, args...)
     QT = typeof(qoldinit)
     return PIDController(map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00))...)
 end
 
-function default_controller(alg::RDPK3Sp49, cache, qoldinit, args...)
+function legacy_default_controller(alg::RDPK3Sp49, cache, qoldinit, args...)
     QT = typeof(qoldinit)
     return PIDController(map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00))...)
 end
 
-function default_controller(alg::RDPK3SpFSAL49, cache, qoldinit, args...)
+function legacy_default_controller(alg::RDPK3SpFSAL49, cache, qoldinit, args...)
     QT = typeof(qoldinit)
     return PIDController(map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01))...)
 end
 
-function default_controller(alg::RDPK3Sp510, cache, qoldinit, args...)
+function legacy_default_controller(alg::RDPK3Sp510, cache, qoldinit, args...)
     QT = typeof(qoldinit)
     return PIDController(map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06))...)
 end
 
-function default_controller(alg::RDPK3SpFSAL510, cache, qoldinit, args...)
+function legacy_default_controller(alg::RDPK3SpFSAL510, cache, qoldinit, args...)
     QT = typeof(qoldinit)
     return PIDController(map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00))...)
 end
+
+function default_controller_v7(QT, alg::RDPK3Sp35)
+    return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.64, -0.31, 0.04)))
+end
+
+function default_controller_v7(QT, alg::RDPK3SpFSAL35)
+    return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.70, -0.23, 0.00)))
+end
+
+function default_controller_v7(QT, alg::RDPK3Sp49)
+    return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.25, -0.12, 0.00)))
+end
+
+function default_controller_v7(QT, alg::RDPK3SpFSAL49)
+    return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.38, -0.18, 0.01)))
+end
+
+function default_controller_v7(QT, alg::RDPK3Sp510)
+    return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.47, -0.20, 0.06)))
+end
+
+function default_controller_v7(QT, alg::RDPK3SpFSAL510)
+    return NewPIDController(QT, alg; beta=map(Base.Fix1(convert, QT), (0.45, -0.13, 0.00)))
+end
diff --git a/lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl b/lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl
index 5db6344b12..79c871d576 100644
--- a/lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl
+++ b/lib/OrdinaryDiffEqNordsieck/src/OrdinaryDiffEqNordsieck.jl
@@ -1,7 +1,7 @@
 module OrdinaryDiffEqNordsieck
 
 import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, qsteady_max_default,
-                           get_current_alg_order,
+                           get_current_alg_order, DummyController,
                            AbstractController, OrdinaryDiffEqAdaptiveAlgorithm,
                            OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm,
                            alg_cache, OrdinaryDiffEqMutableCache,
@@ -11,7 +11,8 @@ import OrdinaryDiffEqCore: alg_order, alg_adaptive_order, qsteady_max_default,
                            calculate_residuals, calculate_residuals!,
                            get_current_adaptive_order, get_fsalfirstlast,
                            ode_interpolant, ode_interpolant!, trivial_limiter!,
-                           generic_solver_docstring
+                           generic_solver_docstring,
+                           default_controller_v7, legacy_default_controller
 using MuladdMacro, FastBroadcast, RecursiveArrayTools
 import LinearAlgebra: rmul!
 import Static: False
diff --git a/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl b/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl
index 7083340c99..f09630bf0c 100644
--- a/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl
+++ b/lib/OrdinaryDiffEqNordsieck/src/alg_utils.jl
@@ -8,6 +8,10 @@ qsteady_max_default(alg::JVODE) = 3 // 2
 
 get_current_alg_order(alg::JVODE, cache) = get_current_adaptive_order(alg, cache)
 
-function default_controller(alg::Union{JVODE}, args...)
+function legacy_default_controller(alg::JVODE, args...)
+    DummyController()
+end
+
+function default_controller_v7(QT, alg::JVODE, args...)
     DummyController()
 end
diff --git a/lib/OrdinaryDiffEqNordsieck/src/algorithms.jl b/lib/OrdinaryDiffEqNordsieck/src/algorithms.jl
index c8e49d97a7..11fca50bfc 100644
--- a/lib/OrdinaryDiffEqNordsieck/src/algorithms.jl
+++ b/lib/OrdinaryDiffEqNordsieck/src/algorithms.jl
@@ -15,17 +15,20 @@ struct AN5 <: OrdinaryDiffEqAdaptiveAlgorithm end
 
     `JVODE` is experimental, the solver `VCABM` is generally preferred.
 """
-struct JVODE{bType, aType} <: OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm
+struct JVODE{bType, aType, qType} <: OrdinaryDiffEqAdamsVarOrderVarStepAlgorithm
     algorithm::Symbol
     bias1::bType
     bias2::bType
     bias3::bType
     addon::aType
+    qmax::qType
+    qsteady_min::qType
+    qsteady_max::qType
 end
 
 function JVODE(algorithm = :Adams; bias1 = 6, bias2 = 6, bias3 = 10,
-        addon = 1 // 10^6)
-    JVODE(algorithm, bias1, bias2, bias3, addon)
+        addon = 1 // 10^6, qmax = float(10), qsteady_min = float(1), qsteady_max = float(3 // 2))
+    JVODE(algorithm, bias1, bias2, bias3, addon, qmax, qsteady_min, qsteady_max)
 end
 """
 !!! warning "Experimental"
diff --git a/lib/OrdinaryDiffEqNordsieck/src/controllers.jl b/lib/OrdinaryDiffEqNordsieck/src/controllers.jl
index 216eb36168..cd6745b67e 100644
--- a/lib/OrdinaryDiffEqNordsieck/src/controllers.jl
+++ b/lib/OrdinaryDiffEqNordsieck/src/controllers.jl
@@ -1,25 +1,22 @@
-struct DummyController <: AbstractController
-end
-
 # JVODE
 function stepsize_controller!(integrator, alg::JVODE)
     if iszero(integrator.EEst)
-        η = integrator.opts.qmax
+        η = alg.qmax
     else
         η = integrator.cache.η
-        integrator.qold = η
+        integrator.cache.ηold = η
     end
     η
 end
 
 function step_accept_controller!(integrator, alg::JVODE, η)
     q = inv(η)
-    if q <= integrator.opts.qsteady_max && q >= integrator.opts.qsteady_min
+    if q <= alg.qsteady_max && q >= alg.qsteady_min
         q = one(q)
     end
     return integrator.dt / q  # dtnew
 end
 
 function step_reject_controller!(integrator, alg::JVODE)
-    integrator.dt *= integrator.qold
+    integrator.dt *= integrator.cache.η
 end
diff --git a/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl b/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl
index a9ecb45215..4c72b0f126 100644
--- a/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl
+++ b/lib/OrdinaryDiffEqNordsieck/src/nordsieck_caches.jl
@@ -133,6 +133,7 @@ mutable struct JVODEConstantCache{zType, lType, dtsType, dType, tsit5Type, etaTy
     n_wait::Int
     # `η` is `dtₙ₊₁/dtₙ`
     η::etaType
+    ηold::etaType
     ηq::etaType
     η₊₁::etaType
     η₋₁::etaType
@@ -154,7 +155,7 @@ function alg_cache(alg::JVODE, u, rate_prototype, ::Type{uEltypeNoUnits},
     η = zero(dt / dt)
     JVODEConstantCache(z, l, m,
         c_LTE₊₁, c_LTE, c_LTE₋₁, c_conv, c_𝒟, prev_𝒟,
-        dts, Δ, tsit5tab, 2, 1, 1, 2, η, η, η, η, η)
+        dts, Δ, tsit5tab, 2, 1, 1, 2, η, η, η, η, η, η)
 end
 
 mutable struct JVODECache{
@@ -205,6 +206,7 @@ mutable struct JVODECache{
     n_wait::Int
     # `η` is `dtₙ₊₁/dtₙ`
     η::etaType
+    ηold::etaType
     ηq::etaType
     η₊₁::etaType
     η₋₁::etaType
@@ -262,7 +264,7 @@ function alg_cache(alg::JVODE, u, rate_prototype, ::Type{uEltypeNoUnits},
     JVODECache(u, uprev, tmp, fsalfirst, ratetmp,
         z, l, m,
         c_LTE₊₁, c_LTE, c_LTE₋₁, c_conv, c_𝒟, prev_𝒟,
-        dts, Δ, atmp, tsit5cache, 2, 1, 1, 2, η, η, η, η, η)
+        dts, Δ, atmp, tsit5cache, 2, 1, 1, 2, η, η, η, η, η, η)
 end
 
 function get_fsalfirstlast(cache::Union{JVODECache, AN5Cache}, u)
diff --git a/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl b/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl
index 4fd7a60f84..f9b7223400 100644
--- a/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl
+++ b/lib/OrdinaryDiffEqVerner/src/OrdinaryDiffEqVerner.jl
@@ -4,7 +4,7 @@ import OrdinaryDiffEqCore: alg_order, calculate_residuals!,
                            initialize!, perform_step!, @unpack, unwrap_alg,
                            calculate_residuals, alg_stability_size,
                            OrdinaryDiffEqAlgorithm,
-                           CompositeAlgorithm, AbstractController, PIDController,
+                           CompositeAlgorithm, accept_step_controller,
                            OrdinaryDiffEqMutableCache, OrdinaryDiffEqConstantCache,
                            OrdinaryDiffEqAdaptiveAlgorithm, CompiledFloats, uses_uprev,
                            alg_cache, _vec, _reshape, @cache, isfsal, full_cache,
@@ -30,7 +30,6 @@ include("verner_caches.jl")
 include("verner_addsteps.jl")
 include("interp_func.jl")
 include("interpolants.jl")
-include("controllers.jl")
 include("verner_rk_perform_step.jl")
 
 import PrecompileTools
diff --git a/lib/OrdinaryDiffEqVerner/src/controllers.jl b/lib/OrdinaryDiffEqVerner/src/controllers.jl
deleted file mode 100644
index e61ba34b6c..0000000000
--- a/lib/OrdinaryDiffEqVerner/src/controllers.jl
+++ /dev/null
@@ -1,8 +0,0 @@
-# checks whether the controller should accept a step based on the error estimate
-@inline function accept_step_controller(integrator, controller::AbstractController)
-    return integrator.EEst <= 1
-end
-
-@inline function accept_step_controller(integrator, controller::PIDController)
-    return integrator.qold >= controller.accept_safety
-end