diff --git a/lib/OrdinaryDiffEqBDF/src/controllers.jl b/lib/OrdinaryDiffEqBDF/src/controllers.jl
index 171078d996b..bfc3dfe9cbd 100644
--- a/lib/OrdinaryDiffEqBDF/src/controllers.jl
+++ b/lib/OrdinaryDiffEqBDF/src/controllers.jl
@@ -83,9 +83,14 @@ end
 function step_accept_controller!(integrator, cache::Union{QNDFCache, QNDFConstantCache}, alg::QNDF{max_order}, q) where {max_order}
     #step is accepted, reset count of consecutive failed steps
     cache.consfailcnt = 0
+    is_disco = integrator.is_disco_step
+    if is_disco
+        integrator.is_disco_step = false
+        cache.nconsteps = 0
+    end
     cache.nconsteps += 1
-    if iszero(OrdinaryDiffEqCore.get_EEst(integrator))
-        return integrator.dt * get_current_qmax(integrator, get_qmax(integrator))
+    new_dt = if iszero(OrdinaryDiffEqCore.get_EEst(integrator))
+        integrator.dt * get_current_qmax(integrator, get_qmax(integrator))
     else
         est = OrdinaryDiffEqCore.get_EEst(integrator)
         estₖ₋₁ = cache.EEst1
@@ -147,16 +152,16 @@ function step_accept_controller!(integrator, cache::Union{QNDFCache, QNDFConstan
         end
         cache.order = kₙ
         q = integrator.dt / hₙ
-    end
-    if prefer_const_step
-        if q < 1.2 && q > 0.6
-            return integrator.dt
+
+        if prefer_const_step && 0.6 < q < 1.2
+            h
+        elseif q <= get_qsteady_max(integrator) && q >= get_qsteady_min(integrator)
+            h
+        else
+            h / q
         end
     end
-    if q <= get_qsteady_max(integrator) && q >= get_qsteady_min(integrator)
-        return integrator.dt
-    end
-    return integrator.dt / q
+    return is_disco ? min((integrator.disco_checkpoint - integrator.t) / 4, new_dt) : new_dt
 end
 
 function bdf_step_reject_controller!(integrator, cache, EEst1)
@@ -174,8 +179,8 @@ function bdf_step_reject_controller!(integrator, cache, EEst1)
     end
 
     if discontinuity_detection
-        disco_dt = set_discontinuity(integrator.u, integrator.uprev, integrator)
-        if disco_dt != -1
+        disco_dt = set_discontinuity(integrator)
+        if disco_dt > zero(disco_dt)
             integrator.dt = disco_dt
             return integrator.dt
         end
@@ -416,6 +421,11 @@ function step_accept_controller!(
         q
     ) where {max_order}
     cache.consfailcnt = 0
+    is_disco = integrator.is_disco_step
+    if is_disco
+        integrator.is_disco_step = false
+        cache.nconsteps = 0
+    end
     if q <= get_qsteady_max(integrator) && q >= get_qsteady_min(integrator)
         q = one(q)
     end
@@ -427,7 +437,8 @@ function step_accept_controller!(
     elseif cache.qwait > 0
         cache.qwait -= 1 # countdown
     end
-    return integrator.dt / q
+    new_dt = integrator.dt / q
+    return is_disco ? min((integrator.disco_checkpoint - integrator.t) / 4, new_dt) : new_dt
 end
 
 function step_reject_controller!(integrator, alg::DFBDF)
@@ -579,6 +590,11 @@ function step_accept_controller!(
         q
     ) where {max_order}
     cache.consfailcnt = 0
+    is_disco = integrator.is_disco_step
+    if is_disco
+        integrator.is_disco_step = false
+        cache.nconsteps = 0
+    end
     if q <= get_qsteady_max(integrator) && q >= get_qsteady_min(integrator)
         q = one(q)
     end
@@ -590,5 +606,6 @@ function step_accept_controller!(
     elseif cache.qwait > 0
         cache.qwait -= 1 # countdown
     end
-    return integrator.dt / q
+    new_dt = integrator.dt / q
+    return is_disco ? min((integrator.disco_checkpoint - integrator.t) / 4, new_dt) : new_dt
 end
diff --git a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
index d1b7cd834cf..8f818665c71 100644
--- a/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
+++ b/lib/OrdinaryDiffEqCore/src/OrdinaryDiffEqCore.jl
@@ -122,6 +122,7 @@ import EnzymeCore
     Tableau        # tableau-derived predictor (α / const_stage_guess)
 end
 export Predictor
+export set_discontinuity
 
 const CompiledFloats = Union{Float32, Float64}
 import Preferences
@@ -261,7 +262,7 @@ include("precompilation_setup.jl")
         # Integrator step / cache / initialization hooks.
         :_ode_init, :_determine_initdt, :ode_determine_initdt, :_initialize_dae!,
         :postamble!, :apply_step!, :last_step_failed, :reset_alg_dependent_opts!,
-        :handle_callback_modifiers!, :set_discontinuity, :resolve_basic,
+        :handle_callback_modifiers!, :resolve_basic,
         # Noise hooks used by the SDE/RODE solver sublibs.
         :accept_noise!, :reinit_noise!, :reject_noise!, :save_noise!, :noise_curt,
         :is_noise_saveable,
diff --git a/lib/OrdinaryDiffEqCore/src/disco.jl b/lib/OrdinaryDiffEqCore/src/disco.jl
index 5c18b9f8952..593ee83be1d 100644
--- a/lib/OrdinaryDiffEqCore/src/disco.jl
+++ b/lib/OrdinaryDiffEqCore/src/disco.jl
@@ -1,8 +1,9 @@
-function set_discontinuity(u, uprev, integrator)
-    breakpointθ = find_discontinuity(u, uprev, integrator)
+function set_discontinuity(integrator)
+    breakpointθ = find_discontinuity(integrator)
     dt = integrator.dt
-    if 1.0e-10 < breakpointθ < 1.0
-        return breakpointθ * dt
+    if breakpointθ < one(breakpointθ)
+        return (0.5 + 0.4 * sin(π * (breakpointθ - 0.5))) * dt
+        return (0.5 + 1.2372 * (breakpointθ - 0.5) - 1.7487 * (breakpointθ - 0.5)^3) * dt
     end
     return -one(dt)
 end
@@ -11,17 +12,20 @@ get_disco_probs(cache::AbstractControllerCache) = cache.controller.basic.disco_p
 get_disco_probs(cache::DummyControllerCache) = cache.disco_probs
 get_disco_probs(cache::CompositeControllerCache) = get_disco_probs(first(cache.caches))
 
-function find_discontinuity(u, uprev, integrator)
+function find_discontinuity(integrator)
     cb = integrator.opts.callback
     dt = integrator.dt
-    cb === nothing && return -one(dt)
-    isempty(cb.continuous_callbacks) && return -one(dt)
+    u = integrator.u
+    uprev = integrator.uprev
+    cb === nothing && return one(dt)
+    isempty(cb.continuous_callbacks) && return one(dt)
     p = integrator.p
     t = integrator.t
     k = integrator.k
-    breakpointθ = -one(dt)
+    breakpointθ = one(dt)
     disco_probs = get_disco_probs(integrator.controller_cache)
     idx = 1
+    addsteps_called = false
     for i in cb.continuous_callbacks
         if (!(i.maybe_discontinuity))
             continue
@@ -40,22 +44,38 @@ function find_discontinuity(u, uprev, integrator)
             i.condition(disco_zero.out_high, u, t + dt, integrator)
             for j in 1:len_cb
                 if (disco_zero.out_low[j] * disco_zero.out_high[j] < zero(disco_zero.out_low[j]))
+                    if (!addsteps_called)
+                        addsteps_called = true
+                        _ode_addsteps!(disco_zero.k, disco_zero.tprev, disco_zero.uprev, disco_zero.u,
+                                    disco_zero.dt, disco_zero.f, disco_zero.p, disco_zero.cache, false, true, false)
+                    end
                     disco_zero.ind = j
+                    disco_prob.tspan[2] = breakpointθ
                     sol = solve(disco_prob)
                     tmp = sol[]
-                    if (!isnan(tmp) && (breakpointθ < zero(breakpointθ) || tmp < breakpointθ))
+                    if (!isnan(tmp) && tmp < breakpointθ)
                         breakpointθ = tmp
+                        integrator.is_disco_step = true
+                        integrator.disco_checkpoint = integrator.t + dt #our prev rejected step, we shouldn't step too far past this
                     end
                 end
             end
         else
-            out_prev = i.condition(uprev, t, integrator)
-            out_curr = i.condition(u, t + dt, integrator)
-            if (out_prev * out_curr < zero(out_prev))
+            disco_zero.out_low[1] = i.condition(uprev, t, integrator)
+            disco_zero.out_high[1] = i.condition(u, t + dt, integrator)
+            if (disco_zero.out_low[1] * disco_zero.out_high[1] < zero(disco_zero.out_low[1]))
+                if (!addsteps_called)
+                    addsteps_called = true
+                    _ode_addsteps!(disco_zero.k, disco_zero.tprev, disco_zero.uprev, disco_zero.u,
+                                disco_zero.dt, disco_zero.f, disco_zero.p, disco_zero.cache, false, true, false)
+                end
+                disco_prob.tspan[2] = breakpointθ
                 sol = solve(disco_prob)
                 tmp = sol[]
-                if (!isnan(tmp) && (breakpointθ < zero(breakpointθ) || tmp < breakpointθ))
+                if (!isnan(tmp) && tmp < breakpointθ)
                     breakpointθ = tmp
+                    integrator.is_disco_step = true
+                    integrator.disco_checkpoint = integrator.t + dt #our prev rejected step, we shouldn't step past this
                 end
             end
         end
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
index c653adfbf30..d511244c99b 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/controllers.jl
@@ -43,6 +43,13 @@ The fields are:
   inside this interval, dt is held constant.
 - `failfactor`: post-Newton-failure shrink factor used by
   [`post_newton_controller!`](@ref).
+- 'discontinuity_detection': If `discontinuity_detection` is set to true, the algorithm will run the autonomous
+    discontinuity detection to predict the best next timestep after step rejection. 
+    Otherwise, it follows the default step rejection algorithm. This feature is currently
+    defaulted off. 
+- `disco_probs`: If `discontinuity_detection` is set to true, this field holds the vector of 
+`IntervalNonlinearProblem`s used for discontinuity detection. Otherwise, it can be left empty.
+
 
 User-supplied overrides flow through controllers as a `NamedTuple` of
 keyword arguments (whatever subset the user passed). At
@@ -485,6 +492,22 @@ for (accessor, default) in (
     end
 end
 
+# Discontinuity-detection helpers 
+# Shared logic in every step_accept/reject_controller! below.
+
+function handle_disco_accept!(integrator, controller_basic, t, nominal_new_dt)
+    controller_basic.discontinuity_detection && integrator.is_disco_step || return nominal_new_dt
+    integrator.is_disco_step = false
+    return min((integrator.disco_checkpoint - t) / 4, nominal_new_dt)
+end
+
+function handle_disco_reject!(integrator, controller_basic)
+    controller_basic.discontinuity_detection || return false
+    disco_dt = set_discontinuity(integrator)
+    disco_dt > zero(disco_dt) || return false
+    integrator.dt = disco_dt
+    return true
+end
 
 # Standard integral (I) step size controller
 """
@@ -510,10 +533,6 @@ the interval `[qmin, qmax]`.
 A step will be accepted whenever the estimated error `get_EEst(integrator)` is
 less than or equal to unity. Otherwise, the step is rejected and re-tried with
 the predicted step size.
-If `discontinuity_detection` is set to true, the algorithm will run the autonomous
-discontinuity detection to predict the best next timestep after step rejection. 
-Otherwise, it follows the default step rejection algorithm. This feature is currently
-defaulted off. 
 
 ## References
 
@@ -569,21 +588,16 @@ end
 function step_accept_controller!(integrator, cache::IControllerCache, alg, q)
     (; qsteady_min, qsteady_max) = cache.controller.basic
 
+    t = integrator.t
+    dt = integrator.dt
     if qsteady_min <= q <= qsteady_max
         q = one(q)
     end
-    return integrator.dt / q # new dt
+    return handle_disco_accept!(integrator, cache.controller.basic, t, dt / q)
 end
 
 function step_reject_controller!(integrator, cache::IControllerCache, alg)
-    discontinuity_detection = cache.controller.basic.discontinuity_detection
-    if discontinuity_detection
-        disco_dt = set_discontinuity(integrator.u, integrator.uprev, integrator)
-        if disco_dt > zero(disco_dt)
-            integrator.dt = disco_dt
-            return integrator.dt
-        end
-    end
+    handle_disco_reject!(integrator, cache.controller.basic) && return integrator.dt
     return integrator.dt = cache.dtreject # TODO this does not look right.
 end
 
@@ -615,10 +629,7 @@ the interval `[qmin, qmax]`.
 A step will be accepted whenever the estimated error `get_EEst(integrator)` is
 less than or equal to unity. Otherwise, the step is rejected and re-tried with
 the predicted step size.
-If `discontinuity_detection` is set to true, the algorithm will run the autonomous
-discontinuity detection to predict the best next timestep after step rejection. 
-Otherwise, it follows the default step rejection algorithm. This feature is currently
-defaulted off. 
+
 !!! note
 
     The coefficients `beta1, beta2` are not scaled by the order of the method,
@@ -712,23 +723,19 @@ function step_accept_controller!(integrator, cache::PIControllerCache, alg, q)
     qoldinit = controller.qoldinit
     EEst = SciMLBase.value(get_EEst(integrator))
 
+    t = integrator.t
+    dt = integrator.dt
     if qsteady_min <= q <= qsteady_max
         q = one(q)
     end
     cache.errold = max(EEst, qoldinit)
-    return integrator.dt / q # new dt
+    return handle_disco_accept!(integrator, controller.basic, t, dt / q)
 end
 
 function step_reject_controller!(integrator, cache::PIControllerCache, alg)
     (; controller, q11) = cache
-    (; qmin, gamma, discontinuity_detection) = controller.basic
-    if discontinuity_detection
-        disco_dt = set_discontinuity(integrator.u, integrator.uprev, integrator)
-        if disco_dt > zero(disco_dt)
-            integrator.dt = disco_dt
-            return integrator.dt
-        end
-    end
+    (; qmin, gamma) = controller.basic
+    handle_disco_reject!(integrator, controller.basic) && return integrator.dt
     return integrator.dt /= min(inv(qmin), q11 / gamma)
 end
 
@@ -781,11 +788,6 @@ Some standard controller parameters suggested in the literature are
 | H211PI     | `1//6`  | `1//6`  | `0`     |
 | H312PID    | `1//18` | `1//9`  | `1//18` |
 
-If `discontinuity_detection` is set to true, the algorithm will run the autonomous
-discontinuity detection to predict the best next timestep after step rejection. 
-Otherwise, it follows the default step rejection algorithm. This feature is currently
-defaulted off. 
-
 !!! note
 
     In contrast to the [`PIController`](@ref), the coefficients `beta1, beta2, beta3`
@@ -947,6 +949,8 @@ function step_accept_controller!(integrator, cache::PIDControllerCache, alg, dt_
     (; controller) = cache
     (; qsteady_min, qsteady_max) = controller.basic
 
+    t = integrator.t
+    dt = integrator.dt
     if qsteady_min <= inv(dt_factor) <= qsteady_max
         dt_factor = one(dt_factor)
     end
@@ -954,18 +958,11 @@ function step_accept_controller!(integrator, cache::PIDControllerCache, alg, dt_
         cache.err[3] = cache.err[2]
         cache.err[2] = cache.err[1]
     end
-    return integrator.dt * dt_factor # new dt
+    return handle_disco_accept!(integrator, controller.basic, t, dt * dt_factor)
 end
 
 function step_reject_controller!(integrator, cache::PIDControllerCache, alg)
-    discontinuity_detection = cache.controller.basic.discontinuity_detection
-    if discontinuity_detection
-        disco_dt = set_discontinuity(integrator.u, integrator.uprev, integrator)
-        if disco_dt > zero(disco_dt)
-            integrator.dt = disco_dt
-            return integrator.dt
-        end
-    end
+    handle_disco_reject!(integrator, cache.controller.basic) && return integrator.dt
     return integrator.dt *= cache.dt_factor
 end
 
@@ -1025,10 +1022,6 @@ integrator.dt / qacc
 ```
 
 When it rejects, it's the same as the [`IController`](@ref):
-If `discontinuity_detection` is set to true, the algorithm will run the autonomous
-discontinuity detection to predict the best next timestep after step rejection. 
-Otherwise, it follows the default step rejection algorithm. This feature is currently
-defaulted off. 
 ```julia
 if integrator.success_iter == 0
     integrator.dt *= 0.1
@@ -1127,20 +1120,13 @@ function step_accept_controller!(integrator, cache::PredictiveControllerCache, a
     cache.dtacc = SciMLBase.value(integrator.dt)
     cache.erracc = max(1.0e-2, EEst)
 
-    return integrator.dt / qacc
+    return handle_disco_accept!(integrator, cache.controller.basic, integrator.t, integrator.dt / qacc)
 end
 
 function step_reject_controller!(integrator, cache::PredictiveControllerCache, alg)
     (; dt, success_iter) = integrator
     (; qold) = cache
-    discontinuity_detection = cache.controller.basic.discontinuity_detection
-    if discontinuity_detection
-        disco_dt = set_discontinuity(integrator.u, integrator.uprev, integrator)
-        if disco_dt > zero(disco_dt)
-            integrator.dt = disco_dt
-            return integrator.dt
-        end
-    end
+    handle_disco_reject!(integrator, cache.controller.basic) && return integrator.dt
     return integrator.dt = success_iter == 0 ? 0.1 * dt : dt / qold
 end
 
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
index b405cdd5ac3..e965920f364 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/integrator_utils.jl
@@ -620,6 +620,18 @@ function fixed_t_for_tstop_error!(integrator, ttmp)
     end
 end
 
+# Type-stable check: did the callback that fired have maybe_discontinuity = true?
+@generated function _fired_cb_maybe_discontinuity(
+        cb_idx,
+        callbacks::NTuple{N, Union{ContinuousCallback, VectorContinuousCallback}}
+    ) where {N}
+    ex = :(false)
+    for i in N:-1:1
+        ex = :(cb_idx == $i ? callbacks[$i].maybe_discontinuity : $ex)
+    end
+    return ex
+end
+
 # Use a generated function to call apply_callback! in a type-stable way
 @generated function apply_ith_callback!(
         integrator,
@@ -681,6 +693,9 @@ function handle_callbacks!(integrator)
                 idx,
                 continuous_callbacks
             )
+            if _fired_cb_maybe_discontinuity(idx, continuous_callbacks)
+                reinit_controller!(integrator, integrator.controller_cache)
+            end
         else
             integrator.event_last_time = 0
             integrator.vector_event_last_time = 1
diff --git a/lib/OrdinaryDiffEqCore/src/integrators/type.jl b/lib/OrdinaryDiffEqCore/src/integrators/type.jl
index 180f1af5ccf..7fc0729e8a0 100644
--- a/lib/OrdinaryDiffEqCore/src/integrators/type.jl
+++ b/lib/OrdinaryDiffEqCore/src/integrators/type.jl
@@ -175,6 +175,8 @@ mutable struct ODEIntegrator{
     fsalfirst::FSALType
     fsallast::FSALType
     rng::RNGType
+    is_disco_step::Bool
+    disco_checkpoint::tType
     W::WType
     P::PType
     sqdt::SqdtType
diff --git a/lib/OrdinaryDiffEqCore/src/solve.jl b/lib/OrdinaryDiffEqCore/src/solve.jl
index e1240f747e6..66c1699fb31 100644
--- a/lib/OrdinaryDiffEqCore/src/solve.jl
+++ b/lib/OrdinaryDiffEqCore/src/solve.jl
@@ -16,7 +16,7 @@ determine_controller_datatype(u, internalnorm, ts::Tuple{<:Number, <:Number}) =
 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, FunctionType, tType2, ParameterType}
+mutable struct zero_func_struct{u1Type, uType, tType, kType, CacheType, idxsType, varsType, callbackType, outType, outCacheType, FunctionType, tType2, ParameterType}
     u₁::u1Type
     callback::callbackType
     dt::tType
@@ -28,8 +28,8 @@ mutable struct zero_func_struct{u1Type, uType, tType, kType, CacheType, idxsType
     differential_vars::varsType
     ind::Int
     out::outType
-    out_low::outType
-    out_high::outType
+    out_low::outCacheType
+    out_high::outCacheType
     f::FunctionType
     tprev::tType2
     p::ParameterType
@@ -38,7 +38,8 @@ end
 parameter_values(z::zero_func_struct) = z.p
 
 function (z::zero_func_struct)(θ, p)
-    _ode_addsteps!(z.k, z.tprev, z.uprev, z.u, z.dt, z.f, z.p, z.cache, false, true, false)
+    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
@@ -672,7 +673,7 @@ Base.@constprop :aggressive function _ode_init(
                 arr = (u isa AbstractArray) ? similar(u, i.len) : zeros(typeof(u), i.len)
                 arr, similar(arr), similar(arr)
             else
-                nothing, nothing, nothing
+                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)
@@ -688,6 +689,8 @@ Base.@constprop :aggressive function _ode_init(
 
     controller_cache = setup_controller_cache(_alg, cache, controller, EEstT, disco_probs)
 
+    is_disco_step = false
+    disco_checkpoint = zero(tType)
     # Seed the initial EEst on the controller cache (was previously
     # `integrator.EEst = oneunit(EEstT)`).
     set_EEst!(controller_cache, EEst)
@@ -721,7 +724,7 @@ Base.@constprop :aggressive function _ode_init(
         isout, reeval_fsal,
         derivative_discontinuity, reinitialize, isdae,
         opts, stats, initializealg, differential_vars,
-        fsalfirst, fsallast, _rng,
+        fsalfirst, fsallast, _rng, is_disco_step, disco_checkpoint,
         W, P, sqdt,
         noise, c, rate_constants
     )
diff --git a/lib/OrdinaryDiffEqCore/test/runtests.jl b/lib/OrdinaryDiffEqCore/test/runtests.jl
index 46342deb64b..1f79736c7a8 100644
--- a/lib/OrdinaryDiffEqCore/test/runtests.jl
+++ b/lib/OrdinaryDiffEqCore/test/runtests.jl
@@ -33,5 +33,4 @@ end
 # Functional tests
 if TEST_GROUP == "Core" || TEST_GROUP == "ALL"
     @time @safetestset "Sparse isdiag Performance" include("sparse_isdiag_tests.jl")
-    @time @safetestset "Discontinuity Detection" include("disco_tests.jl")
 end
diff --git a/lib/OrdinaryDiffEqFIRK/src/controllers.jl b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
index 140707e95de..dc7c3399019 100644
--- a/lib/OrdinaryDiffEqFIRK/src/controllers.jl
+++ b/lib/OrdinaryDiffEqFIRK/src/controllers.jl
@@ -21,8 +21,6 @@ function step_accept_controller!(
     if qsteady_min <= qacc <= qsteady_max
         qacc = one(qacc)
     end
-    ccache.dtacc = integrator.dt
-    ccache.erracc = max(1.0e-2, EEst)
     cache.step = step + 1
     hist_iter = hist_iter * 0.8 + iter * 0.2
     cache.hist_iter = hist_iter
@@ -46,6 +44,12 @@ function step_accept_controller!(
             cache.hist_iter = iter
         end
     end
+    ccache.dtacc = integrator.dt
+    ccache.erracc = max(1.0e-2, EEst)
+    if controller.basic.discontinuity_detection && integrator.is_disco_step
+        integrator.is_disco_step = false
+        return min((integrator.disco_checkpoint - integrator.t) / 4, integrator.dt / qacc)
+    end
     return integrator.dt / qacc
 end
 
@@ -63,7 +67,7 @@ function step_reject_controller!(
     min_stages = (alg.min_order - 1) ÷ 4 * 2 + 1
 
     if (discontinuity_detection)
-        disco_dt = set_discontinuity(integrator.u, integrator.uprev, integrator)
+        disco_dt = set_discontinuity(integrator)
         if disco_dt > zero(dt)
             integrator.dt = disco_dt
         end
diff --git a/lib/OrdinaryDiffEqCore/test/disco_tests.jl b/test/Integrators_I/disco_tests.jl
similarity index 52%
rename from lib/OrdinaryDiffEqCore/test/disco_tests.jl
rename to test/Integrators_I/disco_tests.jl
index 02e9108e724..0310c1df5ed 100644
--- a/lib/OrdinaryDiffEqCore/test/disco_tests.jl
+++ b/test/Integrators_I/disco_tests.jl
@@ -18,26 +18,25 @@ prob = ODEProblem(f, u0, tspan)
 #define callback
 condition(u, t, integrator) = u[1] - 1
 cb = ContinuousCallback(condition, default_affect!; maybe_discontinuity = true)
-cb2 = ContinuousCallback(condition, default_affect!; maybe_discontinuity = false)
 
-sol_disco_radau = solve(prob, RadauIIA5(); callback = cb, reltol = 1.0e-6, controller = predictive_disco_controller(RadauIIA5()))
-#  253.208 μs (7657 allocations: 239.47 KiB)
-sol_no_disco_radau = solve(prob, RadauIIA5(); callback = cb2, reltol = 1.0e-6)
-#  347.917 μs (10308 allocations: 320.86 KiB)
+sol_disco_radau = solve(prob, RadauIIA5(); callback = cb, reltol = 1e-6, controller = predictive_disco_controller(RadauIIA5()))
+#  314.500 μs (8909 allocations: 278.73 KiB)
+sol_no_disco_radau = solve(prob, RadauIIA5(); callback = cb, reltol = 1e-6)
+#  376.125 μs (10628 allocations: 331.61 KiB)
 @test sol_disco_radau.retcode == ReturnCode.Success
 @test sol_disco_radau.stats.nreject <= sol_no_disco_radau.stats.nreject
 
-sol_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb, reltol = 1.0e-6, controller = PI_disco_controller(Rodas5P()))
-#  456.292 μs (16771 allocations: 587.12 KiB)
-sol_no_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb2, reltol = 1.0e-6)
-#  589.666 μs (21804 allocations: 763.61 KiB)
+sol_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb, reltol = 1e-6, controller = PI_disco_controller(Rodas5P()))
+#  444.500 μs (15596 allocations: 547.62 KiB)
+sol_no_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb, reltol = 1e-6)
+#  620.709 μs (21830 allocations: 765.77 KiB)
 @test sol_disco_rosenbrock.retcode == ReturnCode.Success
 @test sol_disco_rosenbrock.stats.nreject <= sol_no_disco_rosenbrock.stats.nreject
 
-sol_disco_tsit5 = solve(prob, Tsit5(); callback = cb, reltol = 1.0e-6, controller = PI_disco_controller(Tsit5()))
-#  48.917 μs (7173 allocations: 227.61 KiB)
-sol_no_disco_tsit5 = solve(prob, Tsit5(); callback = cb2, reltol = 1.0e-6)
-#  48.750 μs (7158 allocations: 227.05 KiB)
+sol_disco_tsit5 = solve(prob, Tsit5(); callback = cb, reltol = 1e-6, controller = PI_disco_controller(Tsit5()))
+#  55.625 μs (7576 allocations: 240.67 KiB)
+sol_no_disco_tsit5 = solve(prob, Tsit5(); callback = cb, reltol = 1e-6)
+#  55.083 μs (7570 allocations: 240.58 KiB)
 @test sol_disco_tsit5.retcode == ReturnCode.Success
 @test sol_disco_tsit5.stats.nreject <= sol_no_disco_tsit5.stats.nreject
 
@@ -60,25 +59,29 @@ prob = ODEProblem(f, u0, tspan)
 #define callbacks
 condition1(u, t, integrator) = u[1] - 1
 cb1 = ContinuousCallback(condition1, default_affect!; maybe_discontinuity = true)
-cb1f = ContinuousCallback(condition1, default_affect!; maybe_discontinuity = false)
 
 condition2(u, t, integrator) = u[1] - 2
 cb2 = ContinuousCallback(condition2, default_affect!; maybe_discontinuity = true)
-cb2f = ContinuousCallback(condition2, default_affect!; maybe_discontinuity = false)
 cb = CallbackSet(cb1, cb2)
-cb2 = CallbackSet(cb1f, cb2f)
 
-sol_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb, reltol = 1.0e-6, controller = PI_disco_controller(Rodas5P()))
-#  1.306 ms (45129 allocations: 1.54 MiB)
-sol_no_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb2, reltol = 1.0e-6)
-#  1.500 ms (52871 allocations: 1.80 MiB)
+sol_disco_radau = solve(prob, RadauIIA5(); callback = cb, reltol = 1e-6, controller = predictive_disco_controller(RadauIIA5()))
+#  930.334 μs (25771 allocations: 802.09 KiB)
+sol_no_disco_radau = solve(prob, RadauIIA5(); callback = cb, reltol = 1e-6)
+#  1.256 ms (34747 allocations: 1.05 MiB)
+@test sol_disco_radau.retcode == ReturnCode.Success
+@test sol_disco_radau.stats.nreject <= sol_no_disco_radau.stats.nreject
+
+sol_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb, reltol = 1e-6, controller = PI_disco_controller(Rodas5P()))
+#  1.224 ms (42698 allocations: 1.46 MiB)
+sol_no_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb, reltol = 1e-6)
+#  1.551 ms (52905 allocations: 1.80 MiB)
 @test sol_disco_rosenbrock.retcode == ReturnCode.Success
 @test sol_disco_rosenbrock.stats.nreject <= sol_no_disco_rosenbrock.stats.nreject
 
-sol_disco_tsit5 = solve(prob, Tsit5(); callback = cb, reltol = 1.0e-6, controller = PI_disco_controller(Tsit5()))
-#  261.375 μs (34373 allocations: 1.06 MiB)
-sol_no_disco_tsit5 = solve(prob, Tsit5(); callback = cb2, reltol = 1.0e-6)
-#  266.125 μs (38745 allocations: 1.19 MiB)
+sol_disco_tsit5 = solve(prob, Tsit5(); callback = cb, reltol = 1e-6, controller = PI_disco_controller(Tsit5()))
+#  267.167 μs (34737 allocations: 1.07 MiB)
+sol_no_disco_tsit5 = solve(prob, Tsit5(); callback = cb, reltol = 1e-6)
+#  279.416 μs (39645 allocations: 1.22 MiB)
 @test sol_disco_tsit5.retcode == ReturnCode.Success
 @test sol_disco_tsit5.stats.nreject <= sol_no_disco_tsit5.stats.nreject
 
@@ -101,71 +104,67 @@ prob_multi = ODEProblem(f_multi_exp!, u0_multi, tspan_multi)
 #define callbacks
 cond_multi_1(u, t, integrator) = u[1] - 0.3
 cb_multi_1 = ContinuousCallback(cond_multi_1, default_affect!; maybe_discontinuity = true)
-cb_multi_1f = ContinuousCallback(cond_multi_1, default_affect!; maybe_discontinuity = false)
 
 cond_multi_2(u, t, integrator) = u[1] - 0.8
 cb_multi_2 = ContinuousCallback(cond_multi_2, default_affect!; maybe_discontinuity = true)
-cb_multi_2f = ContinuousCallback(cond_multi_2, default_affect!; maybe_discontinuity = false)
 
 cb_multi = CallbackSet(cb_multi_1, cb_multi_2)
-cb_multi2 = CallbackSet(cb_multi_1f, cb_multi_2f)
 
-#disco solve
-sol_disco_radau = solve(prob_multi, RadauIIA5(); callback = cb_multi, reltol = 1.0e-7, abstol = 1.0e-9, controller = predictive_disco_controller(RadauIIA5()))
-#  145.542 μs (1641 allocations: 64.64 KiB)
-sol_no_disco_radau = solve(prob_multi, RadauIIA5(); callback = cb_multi2, reltol = 1.0e-7, abstol = 1.0e-9)
-#  125.667 μs (1044 allocations: 46.12 KiB)
+sol_disco_radau = solve(prob_multi, RadauIIA5(); callback=cb_multi, reltol=1e-7, abstol=1e-9, controller = PI_disco_controller(RadauIIA5()))
+#  148.875 μs (1564 allocations: 64.56 KiB)
+sol_no_disco_radau = solve(prob_multi, RadauIIA5(); callback=cb_multi, reltol=1e-7, abstol=1e-9)
+#  129.875 μs (1073 allocations: 48.44 KiB)
 @test sol_disco_radau.retcode == ReturnCode.Success
 @test sol_disco_radau.stats.nreject <= sol_no_disco_radau.stats.nreject
 
-sol_disco_rosenbrock = solve(prob_multi, Rodas5P(); callback = cb_multi, reltol = 1.0e-7, abstol = 1.0e-9, controller = PI_disco_controller(Rodas5P()))
-#  253.291 μs (1758 allocations: 75.62 KiB)
-sol_no_disco_rosenbrock = solve(prob_multi, Rodas5P(); callback = cb_multi2, reltol = 1.0e-7, abstol = 1.0e-9)
-#  239.292 μs (947 allocations: 43.89 KiB)
+sol_disco_rosenbrock = solve(prob_multi, Rodas5P(); callback=cb_multi, reltol=1e-7, abstol=1e-9, controller = PI_disco_controller(Rodas5P()))
+#  209.000 μs (1474 allocations: 60.94 KiB)
+sol_no_disco_rosenbrock = solve(prob_multi, Rodas5P(); callback=cb_multi, reltol=1e-7, abstol=1e-9)
+#  237.250 μs (969 allocations: 45.44 KiB)
 @test sol_disco_rosenbrock.retcode == ReturnCode.Success
 @test sol_disco_rosenbrock.stats.nreject <= sol_no_disco_rosenbrock.stats.nreject
 
-sol_disco_tsit5 = solve(prob_multi, Tsit5(); callback = cb_multi, reltol = 1.0e-7, abstol = 1.0e-9, controller = PI_disco_controller(Tsit5()))
-#  99.250 μs (1603 allocations: 59.85 KiB)
-sol_no_disco_tsit5 = solve(prob_multi, Tsit5(); callback = cb_multi2, reltol = 1.0e-7, abstol = 1.0e-9)
-#  88.000 μs (1105 allocations: 45.27 KiB)
+sol_disco_tsit5 = solve(prob_multi, Tsit5(); callback=cb_multi, reltol=1e-7, abstol=1e-9, controller = PI_disco_controller(Tsit5()))
+#  105.250 μs (1592 allocations: 60.34 KiB)
+sol_no_disco_tsit5 = solve(prob_multi, Tsit5(); callback=cb_multi, reltol = 1e-7, abstol = 1e-9)
+#  88.000 μs (1133 allocations: 46.84 KiB)
 @test sol_disco_tsit5.retcode == ReturnCode.Success
 @test sol_disco_tsit5.stats.nreject <= sol_no_disco_tsit5.stats.nreject
 
-#TEST 4: STIFF DISCONTINUITY
-# very stiff discontinuous system
+#TEST 4: STIFF MULTI-COMPONENT DISCONTINUITY
 function f_stiff_disc!(du, u, p, t)
-    λ = p[1] # stiffness parameter
-    return if u[1] < 0.5
-        du[1] = -λ * u[1] + λ * exp(-t) # stiff decay with forcing
+    λ = p[1]
+    if u[1] < 0.5
+        du[1] = -λ * u[1] + λ * exp(-t)
+        du[2] = -λ * (u[2] - u[1])
+        du[3] = -(u[3] - u[1]^2)
     else
-        du[1] = u[1]
+        du[1] = u[1] - u[2]
+        du[2] = -λ * (u[2] - u[1])
+        du[3] = -(u[3] - u[1])
     end
 end
 
-u0_stiff = [0.1]
-tspan_stiff = (0.0, 3.0)
-prob_stiff = ODEProblem(f_stiff_disc!, u0_stiff, tspan_stiff, [100.0])
+u0_stiff = [0.1, 0.1, 0.01]
+tspan_stiff = (0.0, 2.0)
+prob_stiff = ODEProblem(f_stiff_disc!, u0_stiff, tspan_stiff, [500.0])
 
-#define callback
 cond_stiff(u, t, integrator) = u[1] - 0.5
 cb_stiff = ContinuousCallback(cond_stiff, default_affect!; maybe_discontinuity = true)
-cb_stiff_f = ContinuousCallback(cond_stiff, default_affect!; maybe_discontinuity = false)
 
-#disco solve
-sol_disco_radau = solve(prob_stiff, RadauIIA5(); callback = cb_stiff, reltol = 1.0e-9, abstol = 1.0e-11, controller = predictive_disco_controller(RadauIIA5()))
-#  119.250 μs (1666 allocations: 66.52 KiB)
-sol_no_disco_radau = solve(prob_stiff, RadauIIA5(); callback = cb_stiff_f, reltol = 1.0e-9, abstol = 1.0e-11)
-#  124.167 μs (1476 allocations: 61.16 KiB)
+sol_disco_radau = solve(prob_stiff, RadauIIA5(); callback=cb_stiff, reltol=1e-9, abstol=1e-11, controller = predictive_disco_controller(RadauIIA5()))
+#  207.500 μs (2285 allocations: 101.72 KiB)
+sol_no_disco_radau = solve(prob_stiff, RadauIIA5(); callback=cb_stiff, reltol=1e-9, abstol=1e-11)
+#  196.458 μs (2054 allocations: 95.58 KiB)
 @test sol_disco_radau.retcode == ReturnCode.Success
 @test sol_disco_radau.stats.nreject <= sol_no_disco_radau.stats.nreject
 
-sol_disco_tsit5 = solve(prob_stiff, Tsit5(); callback = cb_stiff, reltol = 1.0e-9, abstol = 1.0e-11, controller = PI_disco_controller(Tsit5()))
-#  80.375 μs (1870 allocations: 70.51 KiB)
-sol_no_disco_tsit5 = solve(prob_stiff, Tsit5(); callback = cb_stiff_f, reltol = 1.0e-9, abstol = 1.0e-11)
-#  79.583 μs (1783 allocations: 68.51 KiB)
-@test sol_disco_tsit5.retcode == ReturnCode.Success
-@test sol_disco_tsit5.stats.nreject <= sol_no_disco_tsit5.stats.nreject
+sol_disco_rosenbrock = solve(prob_stiff, Rodas5P(); callback=cb_stiff, reltol=1e-9, abstol=1e-11, controller = PI_disco_controller(Rodas5P()))
+#  364.833 μs (2030 allocations: 85.34 KiB)
+sol_no_disco_rosenbrock = solve(prob_stiff, Rodas5P(); callback=cb_stiff, reltol=1e-9, abstol=1e-11)
+#  332.750 μs (1765 allocations: 79.78 KiB)
+@test sol_disco_rosenbrock.retcode == ReturnCode.Success
+@test sol_disco_rosenbrock.stats.nreject <= sol_no_disco_rosenbrock.stats.nreject
 
 #TEST 5: DISCONTINUOUS DAE
 # discontinuous DAE with mass matrix
@@ -192,14 +191,13 @@ prob_dae = ODEProblem(f_dae_func, u0_dae, tspan_dae)
 
 cond_dae(u, t, integrator) = u[1] - 0.5
 cb_dae = ContinuousCallback(cond_dae, default_affect!; maybe_discontinuity = true)
-cb_daef = ContinuousCallback(cond_dae, default_affect!; maybe_discontinuity = false)
 
-radau_disco = solve(prob_dae, RadauIIA5(); callback = cb_dae, reltol = 1.0e-8, abstol = 1.0e-10, controller = predictive_disco_controller(RadauIIA5()))
-#  62.875 μs (754 allocations: 33.98 KiB)
-radau_no_disco = solve(prob_dae, RadauIIA5(); callback = cb_daef, reltol = 1.0e-8, abstol = 1.0e-10)
-#  58.417 μs (639 allocations: 30.61 KiB)
-@test radau_disco.retcode == ReturnCode.Success
-@test radau_disco.stats.nreject <= radau_no_disco.stats.nreject
+sol_disco_rosenbrock = solve(prob_dae, Rodas5P(); callback=cb_dae, reltol=1e-8, abstol=1e-10, controller = PI_disco_controller(Rodas5P()))
+#  127.083 μs (831 allocations: 35.53 KiB)
+sol_no_disco_rosenbrock = solve(prob_dae, Rodas5P(); callback=cb_dae, reltol=1e-8, abstol=1e-10)
+#  121.625 μs (625 allocations: 31.69 KiB)
+@test sol_disco_rosenbrock.retcode == ReturnCode.Success
+@test sol_disco_rosenbrock.stats.nreject <= sol_no_disco_rosenbrock.stats.nreject
 
 #TEST 6: VECTOR CALLBACK
 function f_vec_disc!(du, u, p, t)
@@ -218,19 +216,18 @@ function condition_vec!(out, u, t, integrator)
 end
 default_affect_vec!(integrator, idx) = nothing
 cb_vec = VectorContinuousCallback(condition_vec!, default_affect_vec!, 2; maybe_discontinuity = true)
-cb_vec2 = VectorContinuousCallback(condition_vec!, default_affect_vec!, 2; maybe_discontinuity = false)
 
-sol_disco_rosenbrock = solve(prob_vec, Rodas5P(); callback = cb_vec, reltol = 1.0e-7, abstol = 1.0e-9, controller = PI_disco_controller(Rodas5P()))
-#  175.292 μs (1213 allocations: 54.08 KiB)
-sol_no_disco_rosenbrock = solve(prob_vec, Rodas5P(); callback = cb_vec2, reltol = 1.0e-7, abstol = 1.0e-9)
-#  194.709 μs (888 allocations: 46.53 KiB)
+sol_disco_rosenbrock = solve(prob_vec, Rodas5P(); callback = cb_vec, reltol=1e-7, abstol=1e-9, controller = PI_disco_controller(Rodas5P()))
+#  191.250 μs (1212 allocations: 53.62 KiB)
+sol_no_disco_rosenbrock = solve(prob_vec, Rodas5P(); callback = cb_vec, reltol=1e-7, abstol=1e-9)
+#  202.000 μs (939 allocations: 49.47 KiB)
 @test sol_disco_rosenbrock.retcode == ReturnCode.Success
 @test sol_disco_rosenbrock.stats.nreject <= sol_no_disco_rosenbrock.stats.nreject
 
-sol_disco_tsit5 = solve(prob_vec, Tsit5(); callback = cb_vec, reltol = 1.0e-7, abstol = 1.0e-9, controller = PI_disco_controller(Tsit5()))
-#  73.333 μs (1349 allocations: 57.71 KiB)
-sol_no_disco_tsit5 = solve(prob_vec, Tsit5(); callback = cb_vec2, reltol = 1.0e-7, abstol = 1.0e-9)
-#  74.708 μs (1256 allocations: 58.48 KiB)
+sol_disco_tsit5 = solve(prob_vec, Tsit5(); callback = cb_vec, reltol=1e-7, abstol=1e-9, controller = PI_disco_controller(Tsit5()))
+#  94.459 μs (1626 allocations: 67.73 KiB)
+sol_no_disco_tsit5 = solve(prob_vec, Tsit5(); callback = cb_vec, reltol=1e-7, abstol=1e-9)
+#  93.084 μs (1356 allocations: 63.04 KiB)
 @test sol_disco_tsit5.retcode == ReturnCode.Success
 @test sol_disco_tsit5.stats.nreject <= sol_no_disco_tsit5.stats.nreject
 
@@ -251,18 +248,17 @@ prob = ODEProblem(f!, u, tspan)
 
 cond(u, t, integrator) = u[2]
 cb = ContinuousCallback(cond, default_affect!; maybe_discontinuity = true)
-cb2 = ContinuousCallback(cond, default_affect!; maybe_discontinuity = false)
 
-sol_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb, reltol = 1.0e-8, abstol = 1.0e-10, controller = PI_disco_controller(Rodas5P()))
-#  215.792 μs (1621 allocations: 61.61 KiB)
-sol_no_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb2, reltol = 1.0e-8, abstol = 1.0e-10)
-#  172.333 μs (845 allocations: 43.62 KiB)
+sol_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb, reltol = 1e-8, abstol = 1e-10, controller = PI_disco_controller(Rodas5P()))
+#  169.208 μs (1190 allocations: 51.78 KiB)
+sol_no_disco_rosenbrock = solve(prob, Rodas5P(); callback = cb, reltol = 1e-8, abstol = 1e-10)
+#  184.792 μs (877 allocations: 45.72 KiB)
 @test sol_disco_rosenbrock.retcode == ReturnCode.Success
 @test sol_disco_rosenbrock.stats.nreject <= sol_no_disco_rosenbrock.stats.nreject
 
-sol_disco_tsit5 = solve(prob, Tsit5(); callback = cb, reltol = 1.0e-8, abstol = 1.0e-10, controller = PI_disco_controller(Tsit5()))
-#  66.292 μs (1450 allocations: 59.10 KiB)
-sol_no_disco_tsit5 = solve(prob, Tsit5(); callback = cb2, reltol = 1.0e-8, abstol = 1.0e-10)
-#  52.500 μs (1117 allocations: 52.55 KiB)
+sol_disco_tsit5 = solve(prob, Tsit5(); callback = cb, reltol = 1e-8, abstol = 1e-10, controller = PI_disco_controller(Tsit5()))
+#  59.625 μs (1248 allocations: 53.79 KiB)
+sol_no_disco_tsit5 = solve(prob, Tsit5(); callback = cb, reltol = 1e-8, abstol = 1e-10)
+#  54.041 μs (1155 allocations: 54.55 KiB)
 @test sol_disco_tsit5.retcode == ReturnCode.Success
 @test sol_disco_tsit5.stats.nreject <= sol_no_disco_tsit5.stats.nreject