diff --git a/src/ensemblegpukernel/integrators/integrator_utils.jl b/src/ensemblegpukernel/integrators/integrator_utils.jl
index fc3b31c8..bcd9f0fb 100644
--- a/src/ensemblegpukernel/integrators/integrator_utils.jl
+++ b/src/ensemblegpukernel/integrators/integrator_utils.jl
@@ -14,30 +14,49 @@ function build_adaptive_controller_cache(::GPUTsit5IController, ::Type{T}) where
     return T(1 / 5), zero(T), T(5), T(1 / 5), T(9 / 10), T(1.0e-4), T(1.0e-4)
 end
 
+# Adaptive steps land on pending tstops and on `tf` by construction: `_bounded_step`
+# makes a trial step that would reach (or nearly reach) the next landing target exactly
+# the remaining distance, and the accept branch marks `t = tf` only for an accepted step
+# of length exactly `tf - t`. A step whose state is not finite is never accepted. Only a
+# landing on a stop before `tf` retries with the controller's covering step (landing on
+# the stop by interpolation), so no accepted step reaches past `tf`.
+# There is no tolerance-based endpoint snap. (At large |t| a step shortened by error
+# control can still round onto `tf` in `t += dt`; see
+# https://github.com/SciML/DiffEqGPU.jl/issues/567.)
+
 # `qold^beta2` of the PI controller. The I controller has `beta2 = 0`, so the factor is
 # exactly one and the device `pow` call is skipped.
 @inline _qold_factor(alg, qold, beta2) = qold^beta2
 @inline _qold_factor(::GPUTsit5IController, qold, beta2) = one(qold)
 
-# Whether a pending tstop before `tf` falls inside the step just accepted from `integ.t`.
+# Whether a pending tstop before `tf` lies inside the step just accepted from `integ.t`.
 @inline function _tstop_in_step(integ, tf, ::Type{T}) where {T}
     tstops = integ.tstops
     (tstops === nothing || integ.tstops_idx > length(tstops)) && return false
-    stop = @inbounds tstops[integ.tstops_idx]
-    return stop < tf && stop - integ.t - integ.dt - T(100) * eps(T) < T(0)
+    stop = convert(typeof(tf), @inbounds tstops[integ.tstops_idx])
+    return stop < tf && stop - integ.t <= integ.dt
 end
 
 # The next time the adaptive integrator has to land on exactly: the next pending
 # tstop after `integ.t` and before `tf`, otherwise `tf`.
-@inline function _next_stop(integ, tf)
+@inline function _next_stop(integ, tf, ::Type{T}) where {T}
     tstops = integ.tstops
     if tstops !== nothing && integ.tstops_idx <= length(tstops)
-        stop = @inbounds tstops[integ.tstops_idx]
+        stop = convert(typeof(tf), @inbounds tstops[integ.tstops_idx])
         integ.t < stop < tf && return stop
     end
     return tf
 end
 
+# Trial step from `integ.t`: never past the next landing target, and exactly the
+# remaining distance when the proposal reaches it or would leave less than 1% of a
+# step before it.
+@inline function _bounded_step(integ, dt, tf, ::Type{T}) where {T}
+    target = _next_stop(integ, tf, T)
+    remaining = target - integ.t
+    return remaining <= dt + dt / 100 || integ.t + dt >= target ? remaining : dt
+end
+
 @inline function savevalues!(
         integrator::SciMLBase.AbstractODEIntegrator{
             AlgType,
@@ -62,21 +81,35 @@ end
     elseif saveat !== nothing
         saved = true
         savedexactly = true
-        while integrator.cur_t <= length(saveat) && saveat[integrator.cur_t] <= integrator.t
-            savet = saveat[integrator.cur_t]
-            Θ = (savet - integrator.tprev) / integrator.dt
-            @inbounds us[integrator.cur_t] = _ode_interpolant(
-                Θ, integrator.dt,
-                integrator.uprev, integrator
-            )
-            @inbounds ts[integrator.cur_t] = savet
-            integrator.cur_t += 1
-        end
+        _save_saveat!(integrator, ts, us, integrator.t)
     end
 
     return saved, savedexactly
 end
 
+# Save the pending `saveat` points up to `tsave` from the current step. The dense output
+# of several methods uses `integrator.u` as the step's end state, so this must run
+# before `integrator.u` is moved back to an event time.
+@inline function _save_saveat!(integrator, ts, us, tsave)
+    saveat = integrator.saveat
+    while integrator.cur_t <= length(saveat) && saveat[integrator.cur_t] <= tsave
+        savet = saveat[integrator.cur_t]
+        # The step's endpoints are known exactly; only interior points need the
+        # dense output, which is not exact at Θ = 0 or 1 for every method.
+        @inbounds us[integrator.cur_t] = if savet == integrator.t
+            integrator.u
+        elseif savet == integrator.tprev
+            integrator.uprev
+        else
+            Θ = (savet - integrator.tprev) / integrator.dt
+            _ode_interpolant(Θ, integrator.dt, integrator.uprev, integrator)
+        end
+        @inbounds ts[integrator.cur_t] = savet
+        integrator.cur_t += 1
+    end
+    return nothing
+end
+
 @inline function DiffEqBase.terminate!(
         integrator::SciMLBase.AbstractODEIntegrator{
             AlgType,
@@ -263,6 +296,7 @@ end
         cb_time, prev_sign, event_idx, ts,
         us
     ) where {AlgType <: GPUODEAlgorithm, IIP, S, T}
+    integrator.saveat === nothing || _save_saveat!(integrator, ts, us, cb_time)
     DiffEqBase.change_t_via_interpolation!(integrator, cb_time)
 
     # The new absolute-time callback handling can leave dtnew ≈ 0 when
diff --git a/src/ensemblegpukernel/linalg/linsolve.jl b/src/ensemblegpukernel/linalg/linsolve.jl
index 89b0bd64..e0abd041 100644
--- a/src/ensemblegpukernel/linalg/linsolve.jl
+++ b/src/ensemblegpukernel/linalg/linsolve.jl
@@ -14,9 +14,8 @@ end
     @inbounds return similar_type(b, typeof(a[1] \ b[1]))(a[1] \ b[1])
 end
 
-@inline function _linear_solve(
+@inline function _cramer_solve(
         ::Size{(2, 2)},
-        ::Size{(2,)},
         a::StaticMatrix{<:Any, <:Any, Ta},
         b::StaticVector{<:Any, Tb}
     ) where {Ta, Tb}
@@ -28,9 +27,8 @@ end
     )
 end
 
-@inline function _linear_solve(
+@inline function _cramer_solve(
         ::Size{(3, 3)},
-        ::Size{(3,)},
         a::StaticMatrix{<:Any, <:Any, Ta},
         b::StaticVector{<:Any, Tb}
     ) where {Ta, Tb}
@@ -56,6 +54,26 @@ end
     )
 end
 
+# Cramer's rule forms products of up to n matrix entries, which overflow (or underflow)
+# for badly scaled systems such as a Rosenbrock W = I/(γ dt) - J with a tiny dt. When
+# that happens, solve the system divided by its largest entry instead: the solution is
+# unchanged, and well-scaled systems keep their exact unscaled arithmetic.
+for Sa in [(2, 2), (3, 3)]
+    @eval @inline function _linear_solve(
+            ::Size{$Sa},
+            ::Size{($Sa[1],)},
+            a::StaticMatrix{<:Any, <:Any, Ta},
+            b::StaticVector{<:Any, Tb}
+        ) where {Ta, Tb}
+        d = det(a)
+        x = _cramer_solve(Size($Sa), a, b)
+        isfinite(d) && !iszero(d) && all(isfinite, x) && return x
+        s = maximum(abs, a)
+        (isfinite(s) && !iszero(s)) || return x
+        return _cramer_solve(Size($Sa), a / s, b / s)
+    end
+end
+
 for Sa in [(2, 2), (3, 3)]  # not needed for Sa = (1, 1);
     @eval begin
         @inline function _linear_solve(
diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl
index ba4eca59..33b30b60 100644
--- a/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl
+++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno3_perform_step.jl
@@ -99,6 +99,9 @@ end
     t = integ.t
     p = integ.p
     tf = integ.tf
+    dtprop = dt
+    dt = _bounded_step(integ, dt, tf, T)
+    shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf
 
     tmp = integ.tmp
     f = integ.f
@@ -122,7 +125,9 @@ end
     EEst = convert(T, Inf)
 
     while EEst > convert(T, 1.0)
-        dt < convert(T, 1.0f-14) && error("dt<dtmin")
+        (dt < convert(T, 1.0f-14) || t + dt == t) &&
+            dt != _next_stop(integ, tf, T) - t &&
+            error("dt<dtmin")
 
         ## Steps
 
@@ -181,6 +186,9 @@ end
         tmp = (err) ./
             (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
         EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)
+        # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs.
+        finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp)
+        finite || (EEst = T(Inf))
 
         q11 = EEst^beta1
         if iszero(EEst)
@@ -190,11 +198,15 @@ end
         end
 
         if EEst > 1
-            dt = dt / min(inv(qmin), q11 / gamma)
+            # A landing step on a stop that overflowed retries with the controller's covering
+            # step, which lands on the stop by interpolation; it never reaches past `tf`.
+            dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma)
+            shortened_to_land = false
         else # EEst <= 1
             q = max(inv(qmax), min(inv(qmin), q / gamma))
             qold = max(EEst, qoldinit)
             dtnew = dt / q #dtnew
+            shortened_to_land && (dtnew = max(dtnew, dtprop))
 
             @inbounds begin # Necessary for interpolation
                 integ.k1 = k1
@@ -208,18 +220,19 @@ end
 
             if _tstop_in_step(integ, tf, T)
                 integ.t = integ.tstops[integ.tstops_idx]
-                if abs(integ.t - (t + dt)) > eps(integ.t)
+                if integ.t - t != dt
                     integ.u = integ(integ.t)
+                    mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop")
                 end
                 dt = integ.t - integ.tprev
                 integ.tstops_idx += 1
-            elseif (tf - t - dt) < convert(T, 1.0f-14)
+            elseif dt == tf - t
                 integ.t = tf
             else
                 ##Advance the integrator
                 integ.t += dt
             end
-            integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t))
+            integ.dtnew = dtnew
         end
     end
 
diff --git a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl
index 56de80b6..f9d06599 100644
--- a/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl
+++ b/src/ensemblegpukernel/perform_step/gpu_kvaerno5_perform_step.jl
@@ -130,6 +130,9 @@ end
     t = integ.t
     p = integ.p
     tf = integ.tf
+    dtprop = dt
+    dt = _bounded_step(integ, dt, tf, T)
+    shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf
 
     tmp = integ.tmp
     f = integ.f
@@ -155,7 +158,9 @@ end
     EEst = convert(T, Inf)
 
     while EEst > convert(T, 1.0)
-        dt < convert(T, 1.0f-14) && error("dt<dtmin")
+        (dt < convert(T, 1.0f-14) || t + dt == t) &&
+            dt != _next_stop(integ, tf, T) - t &&
+            error("dt<dtmin")
 
         ## Steps
 
@@ -249,6 +254,9 @@ end
         tmp = (err) ./
             (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
         EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)
+        # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs.
+        finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp)
+        finite || (EEst = T(Inf))
 
         q11 = EEst^beta1
         if iszero(EEst)
@@ -258,11 +266,15 @@ end
         end
 
         if EEst > 1
-            dt = dt / min(inv(qmin), q11 / gamma)
+            # A landing step on a stop that overflowed retries with the controller's covering
+            # step, which lands on the stop by interpolation; it never reaches past `tf`.
+            dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma)
+            shortened_to_land = false
         else # EEst <= 1
             q = max(inv(qmax), min(inv(qmin), q / gamma))
             qold = max(EEst, qoldinit)
             dtnew = dt / q #dtnew
+            shortened_to_land && (dtnew = max(dtnew, dtprop))
 
             @inbounds begin # Necessary for interpolation
                 integ.k1 = k1
@@ -276,18 +288,19 @@ end
 
             if _tstop_in_step(integ, tf, T)
                 integ.t = integ.tstops[integ.tstops_idx]
-                if abs(integ.t - (t + dt)) > eps(integ.t)
+                if integ.t - t != dt
                     integ.u = integ(integ.t)
+                    mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop")
                 end
                 dt = integ.t - integ.tprev
                 integ.tstops_idx += 1
-            elseif (tf - t - dt) < convert(T, 1.0f-14)
+            elseif dt == tf - t
                 integ.t = tf
             else
                 ##Advance the integrator
                 integ.t += dt
             end
-            integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t))
+            integ.dtnew = dtnew
         end
     end
 
diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl
index adecbd81..53a012bc 100644
--- a/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl
+++ b/src/ensemblegpukernel/perform_step/gpu_rodas4_perform_step.jl
@@ -126,6 +126,9 @@ end
     t = integ.t
     p = integ.p
     tf = integ.tf
+    dtprop = dt
+    dt = _bounded_step(integ, dt, tf, T)
+    shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf
 
     tmp = integ.tmp
     f = integ.f
@@ -158,7 +161,9 @@ end
     EEst = convert(T, Inf)
 
     while EEst > convert(T, 1.0)
-        dt < convert(T, 1.0f-14) && error("dt<dtmin")
+        (dt < convert(T, 1.0f-14) || t + dt == t) &&
+            dt != _next_stop(integ, tf, T) - t &&
+            error("dt<dtmin")
 
         # Precalculations
         dtC21 = C21 / dt
@@ -225,6 +230,9 @@ end
 
         tmp = k6 ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
         EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)
+        # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs.
+        finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp)
+        finite || (EEst = T(Inf))
 
         q11 = EEst^beta1
         if iszero(EEst)
@@ -234,11 +242,15 @@ end
         end
 
         if EEst > 1
-            dt = dt / min(inv(qmin), q11 / gamma)
+            # A landing step on a stop that overflowed retries with the controller's covering
+            # step, which lands on the stop by interpolation; it never reaches past `tf`.
+            dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma)
+            shortened_to_land = false
         else # EEst <= 1
             q = max(inv(qmax), min(inv(qmin), q / gamma))
             qold = max(EEst, qoldinit)
             dtnew = dt / q #dtnew
+            shortened_to_land && (dtnew = max(dtnew, dtprop))
 
             @inbounds begin # Necessary for interpolation
                 @unpack h21, h22, h23, h24, h25, h31, h32, h33, h34, h35 = integ.tab
@@ -253,18 +265,19 @@ end
 
             if _tstop_in_step(integ, tf, T)
                 integ.t = integ.tstops[integ.tstops_idx]
-                if abs(integ.t - (t + dt)) > eps(integ.t)
+                if integ.t - t != dt
                     integ.u = integ(integ.t)
+                    mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop")
                 end
                 dt = integ.t - integ.tprev
                 integ.tstops_idx += 1
-            elseif (tf - t - dt) < convert(T, 1.0f-14)
+            elseif dt == tf - t
                 integ.t = tf
             else
                 ##Advance the integrator
                 integ.t += dt
             end
-            integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t))
+            integ.dtnew = dtnew
         end
     end
 
diff --git a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl
index 8638caaf..c8bae065 100644
--- a/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl
+++ b/src/ensemblegpukernel/perform_step/gpu_rodas5P_perform_step.jl
@@ -163,6 +163,9 @@ end
     t = integ.t
     p = integ.p
     tf = integ.tf
+    dtprop = dt
+    dt = _bounded_step(integ, dt, tf, T)
+    shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf
 
     tmp = integ.tmp
     f = integ.f
@@ -195,7 +198,9 @@ end
     EEst = convert(T, Inf)
 
     while EEst > convert(T, 1.0)
-        dt < convert(T, 1.0f-14) && error("dt<dtmin")
+        (dt < convert(T, 1.0f-14) || t + dt == t) &&
+            dt != _next_stop(integ, tf, T) - t &&
+            error("dt<dtmin")
 
         # Precalculations
         dtC21 = C21 / dt
@@ -294,6 +299,9 @@ end
 
         tmp = k8 ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
         EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)
+        # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs.
+        finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp)
+        finite || (EEst = T(Inf))
 
         q11 = EEst^beta1
         if iszero(EEst)
@@ -303,11 +311,15 @@ end
         end
 
         if EEst > 1
-            dt = dt / min(inv(qmin), q11 / gamma)
+            # A landing step on a stop that overflowed retries with the controller's covering
+            # step, which lands on the stop by interpolation; it never reaches past `tf`.
+            dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma)
+            shortened_to_land = false
         else # EEst <= 1
             q = max(inv(qmax), min(inv(qmin), q / gamma))
             qold = max(EEst, qoldinit)
             dtnew = dt / q #dtnew
+            shortened_to_land && (dtnew = max(dtnew, dtprop))
 
             @inbounds begin # Necessary for interpolation
                 @unpack h21, h22, h23, h24, h25, h26, h27, h28, h31, h32, h33, h34,
@@ -327,18 +339,19 @@ end
 
             if _tstop_in_step(integ, tf, T)
                 integ.t = integ.tstops[integ.tstops_idx]
-                if abs(integ.t - (t + dt)) > eps(integ.t)
+                if integ.t - t != dt
                     integ.u = integ(integ.t)
+                    mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop")
                 end
                 dt = integ.t - integ.tprev
                 integ.tstops_idx += 1
-            elseif (tf - t - dt) < convert(T, 1.0f-14)
+            elseif dt == tf - t
                 integ.t = tf
             else
                 ##Advance the integrator
                 integ.t += dt
             end
-            integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t))
+            integ.dtnew = dtnew
         end
     end
 
diff --git a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl
index 752f0247..31480bf8 100644
--- a/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl
+++ b/src/ensemblegpukernel/perform_step/gpu_rosenbrock23_perform_step.jl
@@ -81,6 +81,9 @@ end
     t = integ.t
     p = integ.p
     tf = integ.tf
+    dtprop = dt
+    dt = _bounded_step(integ, dt, tf, T)
+    shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf
 
     tmp = integ.tmp
     f = integ.f
@@ -104,7 +107,9 @@ end
     mass_matrix = integ.f.mass_matrix
 
     while EEst > convert(T, 1.0)
-        dt < convert(T, 1.0f-14) && error("dt<dtmin")
+        (dt < convert(T, 1.0f-14) || t + dt == t) &&
+            dt != _next_stop(integ, tf, T) - t &&
+            error("dt<dtmin")
 
         γ = dt * d
 
@@ -151,6 +156,9 @@ end
         tmp = dto6 * (k1 - 2 * k2 + k3)
         tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
         EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)
+        # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs.
+        finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp)
+        finite || (EEst = T(Inf))
 
         q11 = EEst^beta1
         if iszero(EEst)
@@ -160,11 +168,15 @@ end
         end
 
         if EEst > 1
-            dt = dt / min(inv(qmin), q11 / gamma)
+            # A landing step on a stop that overflowed retries with the controller's covering
+            # step, which lands on the stop by interpolation; it never reaches past `tf`.
+            dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma)
+            shortened_to_land = false
         else # EEst <= 1
             q = max(inv(qmax), min(inv(qmin), q / gamma))
             qold = max(EEst, qoldinit)
             dtnew = dt / q #dtnew
+            shortened_to_land && (dtnew = max(dtnew, dtprop))
 
             @inbounds begin # Necessary for interpolation
                 integ.k1 = k1
@@ -178,18 +190,19 @@ end
 
             if _tstop_in_step(integ, tf, T)
                 integ.t = integ.tstops[integ.tstops_idx]
-                if abs(integ.t - (t + dt)) > eps(integ.t)
+                if integ.t - t != dt
                     integ.u = integ(integ.t)
+                    mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop")
                 end
                 dt = integ.t - integ.tprev
                 integ.tstops_idx += 1
-            elseif (tf - t - dt) < convert(T, 1.0f-14)
+            elseif dt == tf - t
                 integ.t = tf
             else
                 ##Advance the integrator
                 integ.t += dt
             end
-            integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t))
+            integ.dtnew = dtnew
         end
     end
     _, saved_in_cb = handle_callbacks!(integ, ts, us)
diff --git a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl
index fd3aa416..4a7f05b9 100644
--- a/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl
+++ b/src/ensemblegpukernel/perform_step/gpu_tsit5_perform_step.jl
@@ -76,6 +76,9 @@ end
     t = integ.t
     p = integ.p
     tf = integ.tf
+    dtprop = dt
+    dt = _bounded_step(integ, dt, tf, T)
+    shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf
     a21, a31, a32, a41, a42, a43, a51, a52, a53, a54,
         a61, a62, a63, a64, a65, a71, a72, a73, a74, a75, a76 = integ.as
     btilde1, btilde2, btilde3, btilde4, btilde5, btilde6, btilde7 = integ.btildes
@@ -99,7 +102,9 @@ end
     EEst = convert(T, Inf)
 
     while EEst > T(1.0)
-        dt < T(1.0e-14) && error("dt<dtmin")
+        (dt < T(1.0e-14) || t + dt == t) &&
+            dt != _next_stop(integ, tf, T) - t &&
+            error("dt<dtmin")
 
         tmp = uprev + dt * a21 * k1
         k2 = f(tmp, p, t + c1 * dt)
@@ -120,6 +125,9 @@ end
         )
         tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
         EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)
+        # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs.
+        finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp)
+        finite || (EEst = T(Inf))
 
         q11 = EEst^beta1
         if iszero(EEst)
@@ -129,11 +137,15 @@ end
         end
 
         if EEst > 1
-            dt = dt / min(inv(qmin), q11 / gamma)
+            # A landing step on a stop that overflowed retries with the controller's covering
+            # step, which lands on the stop by interpolation; it never reaches past `tf`.
+            dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma)
+            shortened_to_land = false
         else # EEst <= 1
             q = max(inv(qmax), min(inv(qmin), q / gamma))
             qold = max(EEst, qoldinit)
             dtnew = dt / q #dtnew
+            shortened_to_land && (dtnew = max(dtnew, dtprop))
 
             @inbounds begin # Necessary for interpolation
                 integ.k1 = k1
@@ -152,18 +164,19 @@ end
 
             if _tstop_in_step(integ, tf, T)
                 integ.t = integ.tstops[integ.tstops_idx]
-                if abs(integ.t - (t + dt)) > eps(integ.t)
+                if integ.t - t != dt
                     integ.u = integ(integ.t)
+                    mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop")
                 end
                 dt = integ.t - integ.tprev
                 integ.tstops_idx += 1
-            elseif (tf - t - dt) < T(1.0e-14)
+            elseif dt == tf - t
                 integ.t = tf
             else
                 ##Advance the integrator
                 integ.t += dt
             end
-            integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t))
+            integ.dtnew = dtnew
         end
     end
     _, saved_in_cb = handle_callbacks!(integ, ts, us)
diff --git a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl
index 28ceeece..cbd1ce08 100644
--- a/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl
+++ b/src/ensemblegpukernel/perform_step/gpu_vern7_perform_step.jl
@@ -86,6 +86,9 @@ end
     t = integ.t
     p = integ.p
     tf = integ.tf
+    dtprop = dt
+    dt = _bounded_step(integ, dt, tf, T)
+    shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf
 
     @unpack c2, c3, c4, c5, c6, c7, c8, a021, a031, a032, a041, a043, a051, a053, a054,
         a061, a063, a064, a065, a071, a073, a074, a075, a076, a081, a083, a084,
@@ -112,7 +115,9 @@ end
     EEst = convert(T, Inf)
 
     while EEst > convert(T, 1.0)
-        dt < convert(T, 1.0f-14) && error("dt<dtmin")
+        (dt < convert(T, 1.0f-14) || t + dt == t) &&
+            dt != _next_stop(integ, tf, T) - t &&
+            error("dt<dtmin")
 
         k1 = f(uprev, p, t)
         a = dt * a021
@@ -153,6 +158,9 @@ end
 
         tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
         EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)
+        # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs.
+        finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp)
+        finite || (EEst = T(Inf))
 
         q11 = EEst^beta1
         if iszero(EEst)
@@ -162,11 +170,15 @@ end
         end
 
         if EEst > 1
-            dt = dt / min(inv(qmin), q11 / gamma)
+            # A landing step on a stop that overflowed retries with the controller's covering
+            # step, which lands on the stop by interpolation; it never reaches past `tf`.
+            dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma)
+            shortened_to_land = false
         else # EEst <= 1
             q = max(inv(qmax), min(inv(qmin), q / gamma))
             qold = max(EEst, qoldinit)
             dtnew = dt / q #dtnew
+            shortened_to_land && (dtnew = max(dtnew, dtprop))
 
             @inbounds begin # Necessary for interpolation
                 integ.k1 = k1
@@ -188,18 +200,19 @@ end
 
             if _tstop_in_step(integ, tf, T)
                 integ.t = integ.tstops[integ.tstops_idx]
-                if abs(integ.t - (t + dt)) > eps(integ.t)
+                if integ.t - t != dt
                     integ.u = integ(integ.t)
+                    mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop")
                 end
                 dt = integ.t - integ.tprev
                 integ.tstops_idx += 1
-            elseif (tf - t - dt) < convert(T, 1.0f-14)
+            elseif dt == tf - t
                 integ.t = tf
             else
                 ##Advance the integrator
                 integ.t += dt
             end
-            integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t))
+            integ.dtnew = dtnew
         end
     end
     _, saved_in_cb = handle_callbacks!(integ, ts, us)
diff --git a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl
index f1952034..62ea89ae 100644
--- a/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl
+++ b/src/ensemblegpukernel/perform_step/gpu_vern9_perform_step.jl
@@ -136,6 +136,9 @@ end
     t = integ.t
     p = integ.p
     tf = integ.tf
+    dtprop = dt
+    dt = _bounded_step(integ, dt, tf, T)
+    shortened_to_land = dt < dtprop && _next_stop(integ, tf, T) < tf
 
     @unpack c1, c2, c3, c4, c5, c6, c7, c8, c9, c10, c11, c12, c13, a0201, a0301, a0302,
         a0401, a0403, a0501, a0503, a0504, a0601, a0604, a0605, a0701, a0704, a0705, a0706,
@@ -166,7 +169,9 @@ end
     EEst = convert(T, Inf)
 
     while EEst > convert(T, 1.0)
-        dt < convert(T, 1.0f-14) && error("dt<dtmin")
+        (dt < convert(T, 1.0f-14) || t + dt == t) &&
+            dt != _next_stop(integ, tf, T) - t &&
+            error("dt<dtmin")
 
         k1 = f(uprev, p, t)
         a = dt * a0201
@@ -262,6 +267,9 @@ end
         )
         tmp = tmp ./ (abstol .+ max.(abs.(uprev), abs.(u)) * reltol)
         EEst = DiffEqBase.ODE_DEFAULT_NORM(tmp, t)
+        # `ODE_DEFAULT_NORM` is computed with fast-math, so NaN is checked on its inputs.
+        finite = mapreduce(isfinite, &, u) & mapreduce(isfinite, &, tmp)
+        finite || (EEst = T(Inf))
 
         q11 = EEst^beta1
         if iszero(EEst)
@@ -271,11 +279,15 @@ end
         end
 
         if EEst > 1
-            dt = dt / min(inv(qmin), q11 / gamma)
+            # A landing step on a stop that overflowed retries with the controller's covering
+            # step, which lands on the stop by interpolation; it never reaches past `tf`.
+            dt = !finite && shortened_to_land ? dtprop : dt / min(inv(qmin), q11 / gamma)
+            shortened_to_land = false
         else # EEst <= 1
             q = max(inv(qmax), min(inv(qmin), q / gamma))
             qold = max(EEst, qoldinit)
             dtnew = dt / q #dtnew
+            shortened_to_land && (dtnew = max(dtnew, dtprop))
 
             @inbounds begin # Necessary for interpolation
                 integ.k1 = k1
@@ -297,18 +309,19 @@ end
 
             if _tstop_in_step(integ, tf, T)
                 integ.t = integ.tstops[integ.tstops_idx]
-                if abs(integ.t - (t + dt)) > eps(integ.t)
+                if integ.t - t != dt
                     integ.u = integ(integ.t)
+                    mapreduce(isfinite, &, integ.u) || error("non-finite state at tstop")
                 end
                 dt = integ.t - integ.tprev
                 integ.tstops_idx += 1
-            elseif (tf - t - dt) < convert(T, 1.0e-14)
+            elseif dt == tf - t
                 integ.t = tf
             else
                 ##Advance the integrator
                 integ.t += dt
             end
-            integ.dtnew = min(abs(dtnew), abs(_next_stop(integ, tf) - integ.t))
+            integ.dtnew = dtnew
         end
     end
     _, saved_in_cb = handle_callbacks!(integ, ts, us)
diff --git a/test/gpu_kernel_de/adaptive_endpoint.jl b/test/gpu_kernel_de/adaptive_endpoint.jl
index d2139453..9a8666a9 100644
--- a/test/gpu_kernel_de/adaptive_endpoint.jl
+++ b/test/gpu_kernel_de/adaptive_endpoint.jl
@@ -1,4 +1,4 @@
-using DiffEqGPU, KernelAbstractions, SciMLBase, StaticArrays, Test
+using DiffEqGPU, KernelAbstractions, LinearAlgebra, SciMLBase, StaticArrays, Test
 
 const ADAPTIVE_ALGS = (
     GPUTsit5(), GPUTsit5IController(), GPUVern7(), GPUVern9(),
@@ -97,13 +97,7 @@ end
         tstops = Float32[0.875], abstol = 1.0f-9, reltol = 1.0f-6, save_everystep = false
     )
     @test all(s -> s.t[end] == tf, sol.u)
-    # Stopping at the tstop goes through the dense output, which for Float32 Vern7
-    # is ~60 ulps off even for u' = 1: https://github.com/SciML/DiffEqGPU.jl/issues/554
-    if alg isa GPUVern7
-        @test all(s -> all(isfinite, s.u[end]), sol.u)
-    else
-        @test all(s -> isapprox(s.u[end], SVector(1.25f0, 1.25f0); rtol = 1.0f-6), sol.u)
-    end
+    @test all(s -> isapprox(s.u[end], SVector(1.25f0, 1.25f0); rtol = 1.0f-6), sol.u)
 end
 
 # After landing on the first stop, the next proposed step reaches `tf`; the
@@ -144,11 +138,407 @@ end
         abstol = 1.0f-9, reltol = 1.0f-6, save_everystep = false
     )
     @test all(s -> s.t[end] == 1.0f0, sol.u)
-    # u0 + (tf - t0) + 10 from the callback; Vern7 is only checked for
-    # the callback having fired until issue 554 is fixed.
-    if alg isa GPUVern7
-        @test all(s -> all(>(11.0f0), s.u[end]), sol.u)
-    else
-        @test all(s -> isapprox(s.u[end], SVector(11.25f0, 11.25f0); rtol = 1.0f-6), sol.u)
+    # u0 + (tf - t0) + 10 from the callback
+    @test all(s -> isapprox(s.u[end], SVector(11.25f0, 11.25f0); rtol = 1.0f-6), sol.u)
+end
+
+# At t ≈ 2^26 a Float32 ulp is 8, so the initial step (32) would end one ulp past the
+# stop in time while having integrated 8 time units beyond it. `tf` keeps every step
+# a whole number of ulps (see https://github.com/SciML/DiffEqGPU.jl/issues/567).
+@testset "Adaptive tstop inside a coarse step, public solve ($(nameof(typeof(alg))))" for
+    alg in ADAPTIVE_ALGS
+    t0 = Float32(2)^26
+    tf = t0 + 32.0f0
+    prob = ODEProblem{false}(
+        (u, p, t) -> SVector(1.0f0, 1.0f0), SVector(0.0f0, 0.0f0), (t0, tf)
+    )
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = 32.0f0, tstops = Float32[t0 + 24.0f0],
+        abstol = 1.0f-6, reltol = 1.0f-3, save_everystep = false
+    )
+    @test all(s -> s.t[end] == tf, sol.u)
+    @test all(s -> s.u[end] ≈ SVector(32.0f0, 32.0f0), sol.u)
+end
+
+# Stops closer together than the minimum step size must both be visited without
+# clamping the next step below `dtmin`.
+@testset "Adaptive adjacent tstops with callbacks, public solve ($(nameof(typeof(alg))))" for
+    alg in ADAPTIVE_ALGS
+    stops = [0.8125, nextfloat(0.8125)]
+    cb = DiscreteCallback(
+        (u, t, integrator) -> t in stops,
+        integrator -> (integrator.u += SVector(10.0, 10.0));
+        save_positions = (false, false)
+    )
+    prob = ODEProblem{false}((u, p, t) -> SVector(1.0, 1.0), SVector(0.0, 0.0), (0.75, 2.0))
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = 0.0625, tstops = stops, callback = cb,
+        merge_callbacks = true, abstol = 1.0e-9, reltol = 1.0e-6, save_everystep = false
+    )
+    @test all(s -> s.t[end] == 2.0, sol.u)
+    @test all(s -> s.u[end] ≈ SVector(21.25, 21.25), sol.u)
+end
+
+# Stops at or within the minimum step size of `tf` must still be visited and their
+# callbacks run; only the forced final interval to `tf` may be shorter than `dtmin`.
+@testset "Adaptive tstops at the endpoint ($(nameof(typeof(alg))), $label)" for
+    alg in ADAPTIVE_ALGS,
+        (label, stops) in (
+            ("at tf", [2.0]), ("one ulp before tf", [prevfloat(2.0)]),
+            ("pair within dtmin of tf", [2.0 - 1.5e-14, 2.0 - 5.0e-15]),
+            ("one ulp before and at tf", [prevfloat(2.0), 2.0]),
+        )
+    cb = DiscreteCallback(
+        (u, t, integrator) -> t in stops,
+        integrator -> (integrator.u += SVector(10.0, 10.0));
+        save_positions = (false, false)
+    )
+    prob = ODEProblem{false}((u, p, t) -> SVector(1.0, 1.0), SVector(0.0, 0.0), (0.75, 2.0))
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = 0.0625, tstops = stops, callback = cb,
+        merge_callbacks = true, abstol = 1.0e-9, reltol = 1.0e-6, save_everystep = false
+    )
+    expected = 1.25 + 10 * length(stops)
+    @test all(s -> s.t[end] == 2.0, sol.u)
+    @test all(s -> isapprox(s.u[end], SVector(expected, expected); rtol = 1.0e-6), sol.u)
+end
+
+@testset "Adaptive terminating callback one ulp before tf ($(nameof(typeof(alg))))" for
+    alg in ADAPTIVE_ALGS
+    stop = prevfloat(2.0)
+    cb = DiscreteCallback(
+        (u, t, integrator) -> t == stop,
+        integrator -> (integrator.u += SVector(10.0, 10.0); terminate!(integrator));
+        save_positions = (false, false)
+    )
+    prob = ODEProblem{false}((u, p, t) -> SVector(1.0, 1.0), SVector(0.0, 0.0), (0.75, 2.0))
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = 1.25 - 2.0e-14, tstops = [stop],
+        callback = cb, merge_callbacks = true, abstol = 1.0e-9, reltol = 1.0e-6,
+        save_everystep = false
+    )
+    expected = stop - 0.75 + 10
+    @test all(s -> s.t[end] == stop, sol.u)
+    @test all(s -> isapprox(s.u[end], SVector(expected, expected); rtol = 1.0e-6), sol.u)
+end
+
+# The step that reaches a stop must keep its dense output intact for the saveat
+# points it covers: one ulp at these times is 8, so a step that overshoots the stop
+# and is then overwritten with the stop's state would corrupt the sample at t0 + 8.
+@testset "Adaptive saveat inside a step that reaches a tstop ($(nameof(typeof(alg))), $T)" for
+    alg in ADAPTIVE_ALGS, T in (Float32, Float64)
+    t0 = T === Float32 ? T(2)^26 : T(2)^55
+    @assert eps(t0) == T(8)
+    prob = ODEProblem{false}(
+        (u, p, t) -> SVector(one(T), one(T)), SVector(zero(T), zero(T)), (t0, t0 + T(32))
+    )
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = T(32), tstops = T[t0 + T(24)],
+        saveat = T[t0, t0 + T(8)], abstol = T(1.0e-6), reltol = T(1.0e-3),
+        save_everystep = false
+    )
+    @test all(s -> s.u[2] ≈ SVector(T(8), T(8)), sol.u)
+end
+
+# A final interval shorter than `dtmin` is only exempt from the `dtmin` check as a
+# whole. If error control rejects it, the solve must either fail with `dt<dtmin` or
+# integrate the full interval; it must never snap to `tf` after a partial step.
+function final_interval_outcome(solve_thunk, exact)
+    sol = try
+        solve_thunk()
+    catch err
+        return err isa ErrorException && occursin("dt<dtmin", err.msg)
+    end
+    return all(s -> isapprox(s.u[end][end], exact; rtol = 1.0e-5, atol = 1.0e-9), sol.u)
+end
+
+@testset "Adaptive rejected sub-dtmin final interval ($(nameof(typeof(alg))))" for
+    alg in ADAPTIVE_ALGS
+    prob = ODEProblem{false}((u, p, t) -> -1.0e15 * u, SVector(1.0, 1.0), (0.0, 4.0e-15))
+    @test final_interval_outcome(exp(-4.0)) do
+        solve(
+            EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+            trajectories = 2, adaptive = true, dt = 4.0e-15, abstol = 1.0e-9,
+            reltol = 1.0e-6, save_everystep = false
+        )
+    end
+end
+
+@testset "Adaptive stiff dynamics switched on just before tf ($(nameof(typeof(alg))))" for
+    alg in ADAPTIVE_ALGS
+    stop, tf = 1.0 - 4.0e-15, 1.0
+    cb = DiscreteCallback(
+        (u, t, integrator) -> t == stop,
+        integrator -> (integrator.u = SVector(1.0e15, 1.0));
+        save_positions = (false, false)
+    )
+    prob = ODEProblem{false}(
+        (u, p, t) -> SVector(0.0, -u[1] * u[2]), SVector(0.0, 1.0), (0.0, tf)
+    )
+    @test final_interval_outcome(exp(-1.0e15 * (tf - stop))) do
+        solve(
+            EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+            trajectories = 2, adaptive = true, dt = 0.5, tstops = [stop], callback = cb,
+            merge_callbacks = true, abstol = 1.0e-9, reltol = 1.0e-6, save_everystep = false
+        )
+    end
+end
+
+# Property: whatever the interval length, stiffness, initial step and stops near the
+# end, an adaptive solve either fails with `dt<dtmin` or reaches the analytic endpoint.
+# Integrating only part of the final interval and reporting success shows up as an
+# error far above the bound: across this grid every method's global error stays below
+# 2000 * reltol (Rosenbrock23, second order, is the largest at about 900 * reltol).
+@testset "Adaptive endpoint property grid ($(nameof(typeof(alg))))" for alg in ADAPTIVE_ALGS
+    decay(u, p, t) = -p[1] * u
+    failures = []
+    for reltol in (1.0e-6, 1.0e-9), L in (1.0e-13, 1.0e-3, 1.0), λL in (0.1, 1.0, 4.0),
+            dtfrac in (0.5, 0.999, 1.0),
+            stops in ([], [L / 4], [L * (1 - 1.0e-2)], [prevfloat(L)], [L - 5.0e-15])
+
+        stops = Float64[s for s in stops if 0 < s < L]
+        prob = ODEProblem{false}(decay, SVector(1.0, 1.0), (0.0, L), SVector(λL / L))
+        endpoint = try
+            sol = solve(
+                EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+                trajectories = 2, adaptive = true, dt = dtfrac * L, tstops = stops,
+                abstol = reltol / 1000, reltol, save_everystep = false
+            )
+            sol.u[1].u[end][1]
+        catch err
+            err isa ErrorException && occursin("dt<dtmin", err.msg) && continue
+            rethrow()
+        end
+        abs(endpoint - exp(-λL)) <= 2000 * reltol * exp(-λL) ||
+            push!(failures, (; reltol, L, λL, dtfrac, stops, endpoint))
+    end
+    @test isempty(failures)
+end
+
+# `tstops` given as a Float64 array for a Float32 problem: the step bound must stay in
+# the integrator's time type.
+@testset "Adaptive Float64 tstops on a Float32 problem ($(nameof(typeof(alg))))" for
+    alg in ADAPTIVE_ALGS
+    prob = ODEProblem{false}(
+        (u, p, t) -> SVector(1.0f0, 1.0f0), SVector(1.0f0, 1.0f0), (0.75f0, 1.0f0)
+    )
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = 0.125f0, tstops = [0.875],
+        abstol = 1.0f-9, reltol = 1.0f-6, save_everystep = false
+    )
+    @test all(s -> s.t[end] == 1.0f0 && s.u[end] ≈ SVector(1.25f0, 1.25f0), sol.u)
+end
+
+# A stop less than `dtmin` after `t` must be reached by an integrated step, not by
+# evaluating the dense output of a longer step near Θ = 1 (inaccurate for Float32
+# Verner interpolants). `u' = (1e14, 0)` is integrated exactly by every method.
+@testset "Adaptive tstop closer than dtmin, public solve ($(nameof(typeof(alg))))" for
+    alg in ADAPTIVE_ALGS
+    dt = 1.0f-14
+    stop = prevfloat(dt)
+    cb = DiscreteCallback(
+        (u, t, integrator) -> t == stop, terminate!; save_positions = (false, false)
+    )
+    prob = ODEProblem{false}(
+        (u, p, t) -> SVector(1.0f14, 0.0f0), SVector(0.0f0, 1.0f0), (0.0f0, 2dt)
+    )
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt, tstops = [stop], callback = cb,
+        merge_callbacks = true, abstol = 1.0f-9, reltol = 1.0f-6, save_everystep = false
+    )
+    exact = Float64(BigFloat(1.0f14) * BigFloat(stop))
+    @test all(s -> s.t[end] == stop && abs(s.u[end][1] - exact) <= 64eps(Float32), sol.u)
+end
+
+# Landing on a stop far below `dtmin` keeps the state accurate for every method family,
+# including Rosenbrock methods whose W = I/(γ dt) - J then has entries near 1/dt.
+@testset "Adaptive tiny landing step ($(nameof(typeof(alg))), rate $rate)" for
+    alg in ADAPTIVE_ALGS, rate in (1.0f14, 1.0f20)
+    stop = 1.0f-20
+    cb = DiscreteCallback(
+        (u, t, integrator) -> t == stop, terminate!; save_positions = (false, false)
+    )
+    prob = ODEProblem{false}(
+        (u, p, t) -> SVector(rate, rate), SVector(0.0f0, 0.0f0), (0.0f0, 2.0f-14)
+    )
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = 1.0f-14, tstops = [stop], callback = cb,
+        merge_callbacks = true, save_everystep = false, abstol = 1.0f-12, reltol = 1.0f-3
+    )
+    exact = Float64(BigFloat(rate) * BigFloat(stop))
+    @test all(
+        s -> s.t[end] == stop &&
+            all(x -> isapprox(Float64(x), exact; rtol = 64eps(Float32)), s.u[end]),
+        sol.u
+    )
+end
+
+# A stop so close that the landing step overflows (Rosenbrock `C/dt` terms) is reached by
+# interpolating the controller's covering step: the result is finite and accurate, never
+# NaN. Vern7's dense output is checked only for finiteness:
+# https://github.com/SciML/DiffEqGPU.jl/issues/554
+@testset "Adaptive stop near floatmin ($(nameof(typeof(alg))), $T)" for
+    alg in ADAPTIVE_ALGS, T in (Float32, Float64)
+    stop = T === Float32 ? 5.0f-38 : 1.0e-307
+    rate = T(1.0e14)
+    cb = DiscreteCallback(
+        (u, t, integrator) -> t == stop, terminate!; save_positions = (false, false)
+    )
+    prob = ODEProblem{false}(
+        (u, p, t) -> SVector(rate, rate), SVector(zero(T), zero(T)), (zero(T), T(2.0e-14))
+    )
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = T(1.0e-14), tstops = [stop], callback = cb,
+        merge_callbacks = true, save_everystep = false, abstol = T(1.0e-12), reltol = T(1.0e-3)
+    )
+    exact = Float64(BigFloat(rate) * BigFloat(stop))
+    @test all(s -> s.t[end] == stop && all(isfinite, s.u[end]), sol.u)
+    if !(alg isa GPUVern7)
+        @test all(
+            s -> all(x -> isapprox(Float64(x), exact; rtol = 128eps(T)), s.u[end]), sol.u
+        )
     end
 end
+
+# With a scaled mass matrix the landing step's W = J - M/(γ dt) overflows even when the
+# stop itself is representable. A step whose state or error estimate is not finite is
+# never accepted: the solver retries with the controller's covering step and lands on
+# the stop by interpolation.
+@testset "Adaptive overflowing landing step with mass matrix ($(nameof(typeof(alg))), $T)" for
+    alg in (GPURosenbrock23(), GPURodas4(), GPURodas5P()), T in (Float32, Float64)
+    stop, rate, mass = T(1024) / floatmax(T), T(1.0e14), T(1024)
+    f = ODEFunction{false}(
+        (u, p, t) -> SVector(rate, rate); mass_matrix = mass * SMatrix{2, 2, T}(I)
+    )
+    prob = ODEProblem{false}(f, zero(SVector{2, T}), (zero(T), T(2.0e-14)))
+    cb = DiscreteCallback(
+        (u, t, integrator) -> t == stop, terminate!; save_positions = (false, false)
+    )
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = T(1.0e-14), tstops = [stop], callback = cb,
+        merge_callbacks = true, save_everystep = false, abstol = T(1.0e-12), reltol = T(1.0e-3)
+    )
+    exact = BigFloat(rate) * BigFloat(stop) / BigFloat(mass)
+    @test all(
+        s -> s.t[end] == stop &&
+            all(x -> isapprox(BigFloat(x), exact; rtol = 128eps(T)), s.u[end]),
+        sol.u
+    )
+end
+
+# `saveat` points at the ends of an accepted step are the step's own states; the dense
+# output is only used strictly inside the step.
+@testset "Adaptive saveat at step endpoints ($(nameof(typeof(alg))), $T)" for
+    alg in ADAPTIVE_ALGS, T in (Float32, Float64)
+    stop = T(0.125)
+    cb = DiscreteCallback(
+        (u, t, integrator) -> t == stop, terminate!; save_positions = (false, false)
+    )
+    prob = ODEProblem{false}(
+        (u, p, t) -> SVector(one(T), one(T)), zero(SVector{2, T}), (zero(T), one(T))
+    )
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = T(0.5), tstops = [stop], callback = cb,
+        merge_callbacks = true, saveat = T[0, stop / 2, stop], abstol = T(1.0e-12),
+        reltol = T(1.0e-6)
+    )
+    @test all(s -> s.t[end] == stop && s.u[1] == zero(SVector{2, T}), sol.u)
+    @test all(s -> all(x -> isapprox(x, stop; rtol = 128eps(T)), s.u[end]), sol.u)
+end
+
+# A covering retry never reaches past `tf`: when the final interval cannot be integrated
+# the solve fails explicitly. The only acceptable outcomes are an accurate state at
+# exactly `tf`, or an explicit failure; success at another time, or with a non-finite
+# state, is never acceptable.
+function endpoint_or_failure(solve_thunk, tf, exact; rtol)
+    sol = try
+        solve_thunk()
+    catch err
+        return err isa ErrorException &&
+            (occursin("dt<dtmin", err.msg) || occursin("non-finite state", err.msg))
+    end
+    return all(sol.u) do s
+        !SciMLBase.successful_retcode(s.retcode) ||
+            s.t[end] == tf &&
+            all(x -> isfinite(x) && isapprox(BigFloat(x), exact; rtol), s.u[end])
+    end
+end
+
+@testset "Adaptive overflowing final landing step ($(nameof(typeof(alg))), $T, $mode)" for
+    alg in (GPURosenbrock23(), GPURodas4(), GPURodas5P()), T in (Float32, Float64),
+        mode in (:tiny, :scaled, :scaled_callback)
+    mass = mode == :tiny ? one(T) : T === Float32 ? T(1.0e30) : T(1.0e300)
+    rate = mode == :tiny ? T(1.0e14) : mass
+    tf = mode == :tiny ? (T === Float32 ? T(5.0e-38) : T(1.0e-307)) : T(1.0e-10)
+    f = ODEFunction{false}(
+        (u, p, t) -> SVector(rate, rate); mass_matrix = mass * SMatrix{2, 2, T}(I)
+    )
+    prob = ODEProblem{false}(f, zero(SVector{2, T}), (zero(T), tf))
+    cb = DiscreteCallback(
+        (u, t, integrator) -> t == tf, integrator -> (integrator.u += SVector(one(T), one(T)));
+        save_positions = (false, false)
+    )
+    kw = mode == :scaled_callback ? (; tstops = T[tf], callback = cb, merge_callbacks = true) :
+        (;)
+    solve_final() = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = mode == :tiny ? T(1.0e-14) : T(0.001),
+        save_everystep = false, abstol = T(1.0e-12), reltol = T(1.0e-3), kw...
+    )
+    exact = BigFloat(rate) * BigFloat(tf) / BigFloat(mass) + (mode == :scaled_callback)
+    @test endpoint_or_failure(tf, exact; rtol = 128eps(T)) do
+        solve_final()
+    end
+end
+
+# A right-hand side that is NaN exactly at `tf` ends either accurately at `tf` or in an
+# explicit failure.
+@testset "Adaptive NaN right-hand side at tf ($(nameof(typeof(alg))), $T)" for
+    alg in ADAPTIVE_ALGS, T in (Float32, Float64)
+    tf = T(0.125)
+    prob = ODEProblem{false}(
+        (u, p, t) -> t == tf ? SVector(T(NaN), T(NaN)) : SVector(one(T), one(T)),
+        zero(SVector{2, T}), (zero(T), tf)
+    )
+    @test endpoint_or_failure(tf, BigFloat(tf); rtol = 1.0e-3) do
+        solve(
+            EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+            trajectories = 2, adaptive = true, dt = T(0.5), save_everystep = false,
+            abstol = T(1.0e-12), reltol = T(1.0e-3)
+        )
+    end
+end
+
+# A `saveat` point inside a step that a ContinuousCallback cuts short at an event is
+# saved from that step's dense output, whose end state is the step's own, not the state
+# moved back to the event.
+@testset "Adaptive saveat before a continuous event ($(nameof(typeof(alg))), $T)" for
+    alg in ADAPTIVE_ALGS, T in (Float32, Float64)
+    cb = ContinuousCallback(
+        (u, t, integrator) -> u[1],
+        integrator -> (integrator.u = SVector(integrator.u[1], -integrator.u[2]));
+        save_positions = (false, false)
+    )
+    prob = ODEProblem{false}(
+        (u, p, t) -> SVector(u[2], zero(T)), SVector{2, T}(1, -1), (zero(T), T(2))
+    )
+    sol = solve(
+        EnsembleProblem(prob), alg, EnsembleGPUKernel(KernelAbstractions.CPU(), 0.0);
+        trajectories = 2, adaptive = true, dt = T(2), callback = cb,
+        merge_callbacks = true, saveat = T[0.5, 2]
+    )
+    @test all(s -> s.t == T[0.5, 2], sol.u)
+    @test all(s -> isapprox(s.u[1], SVector{2, T}(0.5, -1); atol = 1.0e-4), sol.u)
+    @test all(s -> isapprox(s.u[2], SVector{2, T}(1, 1); atol = 1.0e-4), sol.u)
+end