diff --git a/Project.toml b/Project.toml
index b84a1259e3..b6d6522f45 100644
--- a/Project.toml
+++ b/Project.toml
@@ -27,6 +27,7 @@ Logging = "56ddb016-857b-54e1-b83d-db4d58db5568"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 MuladdMacro = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+ODEProblemLibrary = "fdc4e326-1af4-4b90-96e7-779fcce2daa5"
 Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
@@ -43,6 +44,7 @@ SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 StaticArrayInterface = "0d7ed370-da01-4f52-bd93-41d350b8b718"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 TruncatedStacktraces = "781d530d-4396-4725-bb49-402e4bee1e77"
+UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
 
 [compat]
 ADTypes = "0.2, 1"
diff --git a/src/OrdinaryDiffEq.jl b/src/OrdinaryDiffEq.jl
index 3a0506f30d..d0d017dfc1 100644
--- a/src/OrdinaryDiffEq.jl
+++ b/src/OrdinaryDiffEq.jl
@@ -77,6 +77,8 @@ import DiffEqBase: resize!, deleteat!, addat!, full_cache, user_cache, u_cache,
                    add_saveat!, set_reltol!,
                    set_abstol!, postamble!, last_step_failed,
                    isautodifferentiable
+                   
+export change_u!, change_dt!, apriori_bounds_dt
 
 using DiffEqBase: check_error!, @def, _vec, _reshape
 
@@ -128,6 +130,10 @@ import Preferences
 
 DEFAULT_PRECS(W, du, u, p, t, newW, Plprev, Prprev, solverdata) = nothing, nothing
 
+export AbstractPerformStepCallback, NoPerformStepCallback, PerformStepCallback
+
+include("performstep_callback.jl")
+
 include("doc_utils.jl")
 include("misc_utils.jl")
 
@@ -218,6 +224,7 @@ include("perform_step/prk_perform_step.jl")
 include("perform_step/pdirk_perform_step.jl")
 include("perform_step/dae_perform_step.jl")
 include("perform_step/qprk_perform_step.jl")
+include("perform_step/test_for_relaxation_perform_step.jl")
 
 include("dense/generic_dense.jl")
 include("dense/interpolants.jl")
@@ -377,7 +384,7 @@ export FunctionMap, Euler, Heun, Ralston, Midpoint, RK4, ExplicitRK, OwrenZen3,
        OwrenZen5,
        BS3, BS5, RK46NL, DP5, Tsit5, DP8, Vern6, Vern7, Vern8, TanYam7, TsitPap8,
        Vern9, Feagin10, Feagin12, Feagin14, CompositeAlgorithm, Anas5, RKO65, FRK65, PFRK87,
-       RKM, MSRK5, MSRK6, Stepanov5, SIR54, QPRK98, PSRK4p7q6, PSRK3p6q5, PSRK3p5q4
+       RKM, MSRK5, MSRK6, Stepanov5, SIR54, QPRK98, PSRK4p7q6, PSRK3p6q5, PSRK3p5q4, Tsit5_for_relaxation
 
 export SSPRK22, SSPRK33, KYKSSPRK42, SSPRK53, SSPRK53_2N1, SSPRK53_2N2, SSPRK53_H, SSPRK63,
        SSPRK73, SSPRK83, SSPRK43, SSPRK432,
diff --git a/src/alg_utils.jl b/src/alg_utils.jl
index db2a8a414d..a53410d410 100644
--- a/src/alg_utils.jl
+++ b/src/alg_utils.jl
@@ -588,6 +588,7 @@ alg_order(alg::OwrenZen4) = 4
 alg_order(alg::OwrenZen5) = 5
 alg_order(alg::DP5) = 5
 alg_order(alg::Tsit5) = 5
+alg_order(alg::Tsit5_for_relaxation) = 5
 alg_order(alg::DP8) = 8
 alg_order(alg::Vern6) = 6
 alg_order(alg::Vern7) = 7
diff --git a/src/algorithms/explicit_rk.jl b/src/algorithms/explicit_rk.jl
index f18da8a86e..c2a5573332 100644
--- a/src/algorithms/explicit_rk.jl
+++ b/src/algorithms/explicit_rk.jl
@@ -361,6 +361,19 @@ function Tsit5(stage_limiter!, step_limiter! = trivial_limiter!)
     Tsit5(stage_limiter!, step_limiter!, False())
 end
 
+Base.@kwdef struct Tsit5_for_relaxation{StageLimiter, StepLimiter, Thread} <:
+    OrdinaryDiffEqAdaptiveAlgorithm
+stage_limiter!::StageLimiter = trivial_limiter!
+step_limiter!::StepLimiter = trivial_limiter!
+thread::Thread = False()
+end
+
+# for backwards compatibility
+function Tsit5_for_relaxation(stage_limiter!, step_limiter! = trivial_limiter!)
+Tsit5_for_relaxation(stage_limiter!, step_limiter!, False())
+end
+
+
 @doc explicit_rk_docstring(
     "Hairer's 8/5/3 adaption of the Dormand-Prince Runge-Kutta method. (7th order interpolant).",
     "DP8",
diff --git a/src/caches/low_order_rk_caches.jl b/src/caches/low_order_rk_caches.jl
index 54b67ea5b3..8e7685e2e6 100644
--- a/src/caches/low_order_rk_caches.jl
+++ b/src/caches/low_order_rk_caches.jl
@@ -448,6 +448,44 @@ if isdefined(Base, :Experimental) && isdefined(Base.Experimental, :silence!)
     Base.Experimental.silence!(Tsit5Cache)
 end
 
+@cache struct Tsit5Cache_for_relaxation{uType, rateType, uNoUnitsType, StageLimiter, StepLimiter,
+    Thread} <: OrdinaryDiffEqMutableCache
+    u::uType
+    uprev::uType
+    k1::rateType
+    k2::rateType
+    k3::rateType
+    k4::rateType
+    k5::rateType
+    k6::rateType
+    k7::rateType
+    utilde::uType
+    tmp::uType
+    atmp::uNoUnitsType
+    stage_limiter!::StageLimiter
+    step_limiter!::StepLimiter
+    thread::Thread
+end
+
+function alg_cache(alg::Tsit5_for_relaxation, u, rate_prototype, ::Type{uEltypeNoUnits},
+    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+    dt, reltol, p, calck,
+    ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    k1 = zero(rate_prototype)
+    k2 = zero(rate_prototype)
+    k3 = zero(rate_prototype)
+    k4 = zero(rate_prototype)
+    k5 = zero(rate_prototype)
+    k6 = zero(rate_prototype)
+    k7 = zero(rate_prototype)
+    utilde = zero(u)
+    atmp = similar(u, uEltypeNoUnits)
+    recursivefill!(atmp, false)
+    tmp = zero(u)
+    Tsit5Cache_for_relaxation(u, uprev, k1, k2, k3, k4, k5, k6, k7, utilde, tmp, atmp,
+        alg.stage_limiter!, alg.step_limiter!, alg.thread)
+end
+
 @cache struct RK46NLCache{uType, rateType, TabType, StageLimiter, StepLimiter, Thread} <:
               OrdinaryDiffEqMutableCache
     u::uType
@@ -547,6 +585,13 @@ function alg_cache(alg::Tsit5, u, rate_prototype, ::Type{uEltypeNoUnits},
     Tsit5ConstantCache()
 end
 
+function alg_cache(alg::Tsit5_for_relaxation, u, rate_prototype, ::Type{uEltypeNoUnits},
+    ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
+    dt, reltol, p, calck,
+    ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+Tsit5ConstantCache_for_relaxation()
+end
+
 @cache struct DP5Cache{uType, rateType, uNoUnitsType, StageLimiter, StepLimiter,
     Thread} <: OrdinaryDiffEqMutableCache
     u::uType
diff --git a/src/integrators/integrator_interface.jl b/src/integrators/integrator_interface.jl
index 4ba6935cb3..484dbe934c 100644
--- a/src/integrators/integrator_interface.jl
+++ b/src/integrators/integrator_interface.jl
@@ -506,3 +506,22 @@ function DiffEqBase.set_u!(integrator::ODEIntegrator, u)
 end
 
 DiffEqBase.has_stats(i::ODEIntegrator) = true
+
+
+function change_dt!(integrator::ODEIntegrator, dt)
+    integrator.dt_has_changed = true
+    integrator.dt_changed = dt
+end
+
+function change_u!(integrator::ODEIntegrator, u)
+    integrator.u = u
+end
+
+function change_fsallast!(integrator::ODEIntegrator, fsallast)
+    integrator.fsallast = fsallast
+end
+
+has_poststep_callback(integrator::ODEIntegrator) = has_poststep_callback(integrator.opts.performstepcallback)
+has_postfsal_callback(integrator::ODEIntegrator) = has_postfsal_callback(integrator.opts.performstepcallback)
+has_postEEst_callback(integrator::ODEIntegrator) = has_postEEst_callback(integrator.opts.performstepcallback)
+
diff --git a/src/integrators/type.jl b/src/integrators/type.jl
index a361302f9e..c958138842 100644
--- a/src/integrators/type.jl
+++ b/src/integrators/type.jl
@@ -1,6 +1,6 @@
 mutable struct DEOptions{absType, relType, QT, tType, Controller, F1, F2, F3, F4, F5, F6,
     F7, tstopsType, discType, ECType, SType, MI, tcache, savecache,
-    disccache}
+    disccache, MT}
     maxiters::MI
     save_everystep::Bool
     adaptive::Bool
@@ -46,6 +46,7 @@ mutable struct DEOptions{absType, relType, QT, tType, Controller, F1, F2, F3, F4
     force_dtmin::Bool
     advance_to_tstop::Bool
     stop_at_next_tstop::Bool
+    performstepcallback::AbstractPerformStepCallback
 end
 
 TruncatedStacktraces.@truncate_stacktrace DEOptions
@@ -134,6 +135,10 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori
     stats::SciMLBase.DEStats
     initializealg::IA
     differential_vars::DV
+
+    dt_has_changed::Bool
+    dt_changed::tType
+
     fsalfirst::FSALType
     fsallast::FSALType
 
@@ -155,7 +160,7 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori
             accept_step, isout, reeval_fsal, u_modified,
             reinitialize, isdae,
             opts, stats,
-            initializealg, differential_vars) where {algType, IIP, uType,
+            initializealg, differential_vars, dt_has_changed, dt_changed) where {algType, IIP, uType,
             duType, tType, pType,
             eigenType, EEstT,
             tTypeNoUnits, tdirType,
@@ -177,7 +182,7 @@ mutable struct ODEIntegrator{algType <: Union{OrdinaryDiffEqAlgorithm, DAEAlgori
             do_error_check,
             event_last_time, vector_event_last_time, last_event_error,
             accept_step, isout, reeval_fsal, u_modified, reinitialize, isdae,
-            opts, stats, initializealg, differential_vars) # Leave off fsalfirst and last
+            opts, stats, initializealg, differential_vars, dt_has_changed, dt_changed) # Leave off fsalfirst and last
     end
 end
 
diff --git a/src/perform_step/test_for_relaxation_perform_step.jl b/src/perform_step/test_for_relaxation_perform_step.jl
new file mode 100644
index 0000000000..42193c1c99
--- /dev/null
+++ b/src/perform_step/test_for_relaxation_perform_step.jl
@@ -0,0 +1,248 @@
+function initialize!(integrator, ::Tsit5ConstantCache_for_relaxation)
+    integrator.kshortsize = 7
+    integrator.k = typeof(integrator.k)(undef, integrator.kshortsize)
+    integrator.fsalfirst = integrator.f(integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
+    integrator.stats.nf += 1
+    # Avoid undefined entries if k is an array of arrays
+    integrator.fsallast = zero(integrator.fsalfirst)
+    integrator.k[1] = integrator.fsalfirst
+    @inbounds for i in 2:(integrator.kshortsize - 1)
+        integrator.k[i] = zero(integrator.fsalfirst)
+    end
+    integrator.k[integrator.kshortsize] = integrator.fsallast
+end
+
+function initialize!(integrator, cache::Tsit5Cache_for_relaxation)
+    integrator.kshortsize = 7
+    integrator.fsalfirst = cache.k1
+    integrator.fsallast = cache.k7 # setup pointers
+    resize!(integrator.k, integrator.kshortsize)
+    # Setup k pointers
+    integrator.k[1] = cache.k1
+    integrator.k[2] = cache.k2
+    integrator.k[3] = cache.k3
+    integrator.k[4] = cache.k4
+    integrator.k[5] = cache.k5
+    integrator.k[6] = cache.k6
+    integrator.k[7] = cache.k7
+    integrator.f(integrator.fsalfirst, integrator.uprev, integrator.p, integrator.t) # Pre-start fsal
+    integrator.stats.nf += 1
+    return nothing
+end
+
+function perform_step!(integrator, cache::Union{Tsit5Cache_for_relaxation,Tsit5ConstantCache_for_relaxation}, repeat_step = false)
+
+    # Caculate uₙ₊₁
+    calculate_nextstep!(integrator, cache, repeat_step)
+
+    # Perform customize modification right after the step, i.e modify uₙ₊₁
+    if has_poststep_callback(integrator)
+        apply_poststep_callback!(integrator)
+    end
+
+    # Calculate f(uₙ₊₁) if the aglortihm has the FSAL property
+    if isfsal(integrator.alg)
+        calculate_fsal!(integrator, cache)
+    end
+
+    # Perform customize modification right after fsal, i.e modify f(uₙ₊₁)
+    if has_postfsal_callback(integrator)
+        apply_postfsal_callback!(integrator)
+    end
+
+    # Calculate the error estimate needed for PID controller
+    if integrator.opts.adaptive
+        calculate_EEst!(integrator, cache)
+    end
+
+    #Perform customized modification right after error estimate
+    if has_postEEst_callback(integrator)
+        apply_postEEst_callback!(integrator)
+    end
+
+    # finalize_step (useless for the moment)
+    finalize_step!(integrator, cache)
+end
+
+function apply_poststep_callback!(integrator)
+
+    # Variable to know if dt has changed during perform_step
+    integrator.dt_has_changed = false
+
+    integrator.opts.performstepcallback.poststep_cb(integrator)
+
+    handle_dt_bound_and_tstop!(integrator)
+end
+
+function apply_postfsal_callback!(integrator)
+
+    # Variable to know if dt has changed during perform_step
+    integrator.dt_has_changed = false
+
+    integrator.opts.performstepcallback.postfsal_cb(integrator)
+
+    handle_dt_bound_and_tstop!(integrator)
+end
+
+function apply_postEEst_callback!(integrator)
+
+    # Variable to know if dt has changed during perform_step
+    integrator.dt_has_changed = false
+
+    integrator.opts.performstepcallback.postEEst_cb(integrator)
+
+    handle_dt_bound_and_tstop!(integrator)
+end
+
+function handle_dt_bound_and_tstop!(integrator)
+    
+    # Here we carry of the dt modified by the user in this step, if it has been changed.
+
+    if integrator.dt_has_changed
+
+        # Match dt in [dtmin, dtmax]
+        if integrator.tdir > 0
+            integrator.dt_changed = min(integrator.opts.dtmax, integrator.dt_changed)
+        else
+            integrator.dt_changed = max(integrator.opts.dtmax, integrator.dt_changed)
+        end
+        dtmin = timedepentdtmin(integrator)
+        if integrator.tdir > 0
+            integrator.dt_changed = max(integrator.dt_changed, dtmin)
+        else
+            integrator.dt_changed = min(integrator.dt_changed, dtmin)
+        end
+        # Match dt with tstops
+        if has_tstop(integrator)
+            tdir_t = integrator.tdir * integrator.t
+            tdir_tstop = first_tstop(integrator)
+            integrator.dt_changed = integrator.tdir * min(abs(integrator.dt_changed), abs(tdir_tstop - tdir_t)) 
+        end
+    
+        integrator.dt = integrator.dt_changed
+        integrator.dt_has_changed = false
+    end
+end
+
+@muladd function calculate_nextstep!(integrator, ::Tsit5ConstantCache_for_relaxation, repeat_step = false)
+    
+    @unpack t, dt, uprev, u, f, p = integrator
+    T = constvalue(recursive_unitless_bottom_eltype(u))
+    T2 = constvalue(typeof(one(t)))
+    @OnDemandTableauExtract Tsit5ConstantCacheActual T T2
+
+    k1 = integrator.fsalfirst
+    k2 = f(uprev + dt * a21 * k1, p, t + c1 * dt)
+    k3 = f(uprev + dt * (a31 * k1 + a32 * k2), p, t + c2 * dt)
+    k4 = f(uprev + dt * (a41 * k1 + a42 * k2 + a43 * k3), p, t + c3 * dt)
+    k5 = f(uprev + dt * (a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4), p, t + c4 * dt)
+    k6 = f(uprev + dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5), p, t + dt)
+    u = uprev + dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6)
+    
+    integrator.k[1] = k1
+    integrator.k[2] = k2
+    integrator.k[3] = k3
+    integrator.k[4] = k4
+    integrator.k[5] = k5
+    integrator.k[6] = k6
+    integrator.u = u
+
+    integrator.stats.nf += 5
+end
+
+@muladd function calculate_fsal!(integrator, ::Tsit5ConstantCache_for_relaxation)
+    @unpack t, dt, u, f, p = integrator
+    k7 = f(u, p, t + dt)
+    integrator.k[7] = k7
+    integrator.fsallast = k7
+    integrator.stats.nf += 1
+end
+
+@muladd function calculate_EEst!(integrator, ::Tsit5ConstantCache_for_relaxation)
+    T = constvalue(recursive_unitless_bottom_eltype(integrator.u))
+    T2 = constvalue(typeof(one(integrator.t)))
+    @OnDemandTableauExtract Tsit5ConstantCacheActual T T2
+    utilde = integrator.dt * (btilde1 * integrator.k[1] + btilde2 * integrator.k[2] 
+            + btilde3 * integrator.k[3] + btilde4 * integrator.k[4] 
+            + btilde5 * integrator.k[5] + btilde6 * integrator.k[6] 
+            + btilde7 * integrator.fsallast)
+    atmp = calculate_residuals(utilde, integrator.uprev, integrator.u, integrator.opts.abstol,
+    integrator.opts.reltol, integrator.opts.internalnorm, integrator.t)
+    integrator.EEst = integrator.opts.internalnorm(atmp, integrator.t)
+end
+
+function finalize_step!(integrator, ::Tsit5ConstantCache_for_relaxation)
+
+end
+
+
+## Non Constant cache
+
+@muladd function calculate_nextstep!(integrator, cache::Tsit5Cache_for_relaxation, repeat_step = false)
+    
+    @unpack t, dt, uprev, u, f, p = integrator
+    T = constvalue(recursive_unitless_bottom_eltype(u))
+    T2 = constvalue(typeof(one(t)))
+    @OnDemandTableauExtract Tsit5ConstantCacheActual T T2
+    @unpack k1, k2, k3, k4, k5, k6, k7, utilde, tmp, atmp, stage_limiter!, step_limiter!, thread = cache
+
+    @.. broadcast=false thread=thread tmp = uprev + dt * a21 * k1
+    stage_limiter!(tmp, f, p, t + c1 * dt)
+    f(k2, tmp, p, t + c1 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a31 * k1 + a32 * k2)
+    stage_limiter!(tmp, f, p, t + c2 * dt)
+    f(k3, tmp, p, t + c2 * dt)
+    @.. broadcast=false thread=thread tmp = uprev + dt * (a41 * k1 + a42 * k2 + a43 * k3)
+    stage_limiter!(tmp, f, p, t + c3 * dt)
+    f(k4, tmp, p, t + c3 * dt)
+    @.. broadcast=false thread=thread tmp = uprev +
+                                          dt * (a51 * k1 + a52 * k2 + a53 * k3 + a54 * k4)
+    stage_limiter!(tmp, f, p, t + c4 * dt)
+    f(k5, tmp, p, t + c4 * dt)
+    @.. broadcast=false thread=thread tmp = uprev +
+                            dt * (a61 * k1 + a62 * k2 + a63 * k3 + a64 * k4 + a65 * k5)
+    stage_limiter!(tmp, f, p, t + dt)
+    f(k6, tmp, p, t + dt)
+    @.. broadcast=false thread=thread u = uprev +
+                    dt * (a71 * k1 + a72 * k2 + a73 * k3 + a74 * k4 + a75 * k5 + a76 * k6)
+    stage_limiter!(u, integrator, p, t + dt)
+    step_limiter!(u, integrator, p, t + dt)
+
+    integrator.stats.nf += 5
+    return nothing
+end
+
+@muladd function calculate_fsal!(integrator, cache::Tsit5Cache_for_relaxation)
+    @unpack t, dt, u, f, p = integrator
+    f(cache.k7, u, p, t + dt)
+    integrator.stats.nf += 1
+end
+
+@muladd function calculate_EEst!(integrator, cache::Tsit5Cache_for_relaxation)
+    T = constvalue(recursive_unitless_bottom_eltype(integrator.u))
+    T2 = constvalue(typeof(one(integrator.t)))
+    @OnDemandTableauExtract Tsit5ConstantCacheActual T T2
+    @unpack k1, k2, k3, k4, k5, k6, k7, utilde, atmp, thread = cache
+    @.. broadcast = false thread = threa utilde = integrator.dt * 
+        (btilde1 * k1 + btilde2 * k2 + btilde3 * k3 + btilde4 * k4 +
+         btilde5 * k5 + btilde6 * k6 + btilde7 * k7)
+    calculate_residuals!(atmp, utilde, integrator.uprev, integrator.u,
+                         integrator.opts.abstol, integrator.opts.reltol, 
+                         integrator.opts.internalnorm, integrator.t, thread)
+    integrator.EEst = integrator.opts.internalnorm(atmp, integrator.t)
+end
+
+function finalize_step!(integrator, cache::Tsit5Cache_for_relaxation)
+
+end
+
+# Is this useful ?
+function apriori_bounds_dt(integrator)
+    dt_sup = if has_tstop(integrator)
+         integrator.tdir * min(abs(integrator.opts.dtmax) , abs(first_tstop(integrator) - integrator.t))
+    else
+        integrator.tdir * abs(integrator.opts.dtmax)
+    end
+    dt_inf = integrator.tdir * timedepentdtmin(integrator)
+    (dt_inf,dt_sup)
+end
\ No newline at end of file
diff --git a/src/performstep_callback.jl b/src/performstep_callback.jl
new file mode 100644
index 0000000000..85908f7c56
--- /dev/null
+++ b/src/performstep_callback.jl
@@ -0,0 +1,28 @@
+abstract type AbstractPerformStepCallback end
+
+struct NoPerformStepCallback <:AbstractPerformStepCallback end
+
+has_poststep_callback(::AbstractPerformStepCallback) = false
+has_postfsal_callback(::AbstractPerformStepCallback) = false
+has_postEEst_callback(::AbstractPerformStepCallback) = false
+
+struct PerformStepCallback <: AbstractPerformStepCallback
+    poststep_cb
+    has_poststep_cb::Bool
+    postfsal_cb
+    has_postfsal_cb::Bool
+    postEEst_cb
+    has_postEEst_cb::Bool
+
+    function PerformStepCallback(;poststep = nothing, postfsal = nothing, postEEst = nothing)
+        has_poststep_cb = poststep isa Nothing ? false : true
+        has_postfsal_cb = postfsal isa Nothing ? false : true
+        has_postEEst_cb = postEEst isa Nothing ? false : true
+        new(poststep, has_poststep_cb, postfsal, has_postfsal_cb, postEEst, has_postEEst_cb)
+    end
+end
+
+has_poststep_callback(pscb::PerformStepCallback) = pscb.has_poststep_cb
+has_postfsal_callback(pscb::PerformStepCallback) = pscb.has_postfsal_cb
+has_postEEst_callback(pscb::PerformStepCallback) = pscb.has_postEEst_cb
+
diff --git a/src/solve.jl b/src/solve.jl
index 6aee75aaf8..e7afc32133 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -72,6 +72,7 @@ function DiffEqBase.__init(
         alias_u0 = false,
         alias_du0 = false,
         initializealg = DefaultInit(),
+        modif = NoPerformStepCallback(),
         kwargs...) where {recompile_flag}
     if prob isa DiffEqBase.AbstractDAEProblem && alg isa OrdinaryDiffEqAlgorithm
         error("You cannot use an ODE Algorithm with a DAEProblem")
@@ -377,7 +378,7 @@ function DiffEqBase.__init(
         typeof(d_discontinuities_internal), typeof(userdata),
         typeof(save_idxs),
         typeof(maxiters), typeof(tstops),
-        typeof(saveat), typeof(d_discontinuities)}(maxiters, save_everystep,
+        typeof(saveat), typeof(d_discontinuities), typeof(modif)}(maxiters, save_everystep,
         adaptive, abstol_internal,
         reltol_internal,
         QT(gamma), QT(qmax),
@@ -409,7 +410,8 @@ function DiffEqBase.__init(
         unstable_check,
         verbose, calck, force_dtmin,
         advance_to_tstop,
-        stop_at_next_tstop)
+        stop_at_next_tstop,
+        modif)
 
     stats = SciMLBase.DEStats(0)
     differential_vars = prob isa DAEProblem ? prob.differential_vars :
@@ -442,6 +444,8 @@ function DiffEqBase.__init(
     tprev = t
     dtcache = tType(dt)
     dtpropose = tType(dt)
+    dt_has_changed = false
+    dt_changed = tType(dt)
     iter = 0
     kshortsize = 0
     reeval_fsal = false
@@ -489,7 +493,7 @@ function DiffEqBase.__init(
         last_event_error, accept_step,
         isout, reeval_fsal,
         u_modified, reinitiailize, isdae,
-        opts, stats, initializealg, differential_vars)
+        opts, stats, initializealg, differential_vars, dt_has_changed, dt_changed)
 
     if initialize_integrator
         if isdae || SciMLBase.has_initializeprob(prob.f)
diff --git a/src/tableaus/low_order_rk_tableaus.jl b/src/tableaus/low_order_rk_tableaus.jl
index 0dde8ad263..c62d6527d6 100644
--- a/src/tableaus/low_order_rk_tableaus.jl
+++ b/src/tableaus/low_order_rk_tableaus.jl
@@ -499,6 +499,7 @@ function OwrenZen5ConstantCache(T, T2)
 end
 
 struct Tsit5ConstantCache <: OrdinaryDiffEqConstantCache end
+struct Tsit5ConstantCache_for_relaxation <: OrdinaryDiffEqConstantCache end
 struct Tsit5ConstantCacheActual{T, T2}
     c1::T2
     c2::T2
diff --git a/test/relaxation/benchmark in time.ipynb b/test/relaxation/benchmark in time.ipynb
new file mode 100644
index 0000000000..5b45287870
--- /dev/null
+++ b/test/relaxation/benchmark in time.ipynb	
@@ -0,0 +1,255 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Study on time for solving"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 21,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using OrdinaryDiffEq, DiffEqDevTools,  Test, BenchmarkTools\n",
+    "\n",
+    "import ODEProblemLibrary: prob_ode_linear,\n",
+    "                          prob_ode_2Dlinear\n",
+    "\n",
+    "\n",
+    "include(\"relaxation.jl\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Check if the new structure is not slowing down the code\n",
+    "\n",
+    "For the previous structure :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 18,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "BenchmarkTools.Trial: 10000 samples with 10 evaluations.\n",
+       " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m):  \u001b[39m\u001b[36m\u001b[1m1.740 μs\u001b[22m\u001b[39m … \u001b[35m540.050 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m 0.00% … 98.64%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m):     \u001b[39m\u001b[34m\u001b[1m2.300 μs               \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m):    \u001b[39m 0.00%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m):   \u001b[39m\u001b[32m\u001b[1m2.801 μs\u001b[22m\u001b[39m ± \u001b[32m 10.732 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m):  \u001b[39m11.27% ±  3.09%\n",
+       "\n",
+       "  \u001b[39m \u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▇\u001b[39m█\u001b[34m█\u001b[39m\u001b[39m▇\u001b[39m▅\u001b[39m▄\u001b[39m▃\u001b[39m▃\u001b[39m▃\u001b[39m▃\u001b[32m▃\u001b[39m\u001b[39m▂\u001b[39m▃\u001b[39m▂\u001b[39m▂\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m \u001b[39m \u001b[39m▁\u001b[39m▁\u001b[39m \u001b[39m \u001b[39m▁\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▃\n",
+       "  \u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[34m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[32m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m█\u001b[39m▇\u001b[39m█\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▄\u001b[39m▄\u001b[39m▃\u001b[39m▅\u001b[39m▂\u001b[39m▆\u001b[39m \u001b[39m█\n",
+       "  1.74 μs\u001b[90m      \u001b[39m\u001b[90mHistogram: \u001b[39m\u001b[90m\u001b[1mlog(\u001b[22m\u001b[39m\u001b[90mfrequency\u001b[39m\u001b[90m\u001b[1m)\u001b[22m\u001b[39m\u001b[90m by time\u001b[39m      5.16 μs \u001b[0m\u001b[1m<\u001b[22m\n",
+       "\n",
+       " Memory estimate\u001b[90m: \u001b[39m\u001b[33m4.47 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m36\u001b[39m."
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "@benchmark solve(prob_ode_linear, Tsit5())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 15,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "BenchmarkTools.Trial: 10000 samples with 6 evaluations.\n",
+       " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m):  \u001b[39m\u001b[36m\u001b[1m5.933 μs\u001b[22m\u001b[39m … \u001b[35m397.350 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.00% … 94.93%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m):     \u001b[39m\u001b[34m\u001b[1m6.567 μs               \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m):    \u001b[39m0.00%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m):   \u001b[39m\u001b[32m\u001b[1m7.658 μs\u001b[22m\u001b[39m ± \u001b[32m 11.116 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m):  \u001b[39m5.87% ±  4.02%\n",
+       "\n",
+       "  \u001b[39m▃\u001b[39m▇\u001b[39m█\u001b[39m▇\u001b[39m▅\u001b[34m▅\u001b[39m\u001b[39m▅\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▅\u001b[32m▅\u001b[39m\u001b[39m▅\u001b[39m▅\u001b[39m▄\u001b[39m▂\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m \u001b[39m▁\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▂\n",
+       "  \u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[34m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[32m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▆\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m▇\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▄\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▆\u001b[39m▇\u001b[39m▇\u001b[39m▆\u001b[39m▆\u001b[39m▄\u001b[39m \u001b[39m█\n",
+       "  5.93 μs\u001b[90m      \u001b[39m\u001b[90mHistogram: \u001b[39m\u001b[90m\u001b[1mlog(\u001b[22m\u001b[39m\u001b[90mfrequency\u001b[39m\u001b[90m\u001b[1m)\u001b[22m\u001b[39m\u001b[90m by time\u001b[39m      13.6 μs \u001b[0m\u001b[1m<\u001b[22m\n",
+       "\n",
+       " Memory estimate\u001b[90m: \u001b[39m\u001b[33m14.70 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m110\u001b[39m."
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "@benchmark solve(prob_ode_2Dlinear, Tsit5())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For the new structure :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 17,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "BenchmarkTools.Trial: 10000 samples with 10 evaluations.\n",
+       " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m):  \u001b[39m\u001b[36m\u001b[1m1.770 μs\u001b[22m\u001b[39m … \u001b[35m505.320 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m 0.00% … 97.31%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m):     \u001b[39m\u001b[34m\u001b[1m2.000 μs               \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m):    \u001b[39m 0.00%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m):   \u001b[39m\u001b[32m\u001b[1m2.572 μs\u001b[22m\u001b[39m ± \u001b[32m 10.111 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m):  \u001b[39m11.65% ±  3.08%\n",
+       "\n",
+       "  \u001b[39m▃\u001b[39m▇\u001b[39m█\u001b[39m█\u001b[39m█\u001b[34m▆\u001b[39m\u001b[39m▅\u001b[39m▄\u001b[39m▄\u001b[39m▄\u001b[39m▅\u001b[39m▆\u001b[39m▅\u001b[39m▄\u001b[39m▂\u001b[39m▁\u001b[39m▂\u001b[32m▁\u001b[39m\u001b[39m▁\u001b[39m \u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▂\u001b[39m▂\u001b[39m▁\u001b[39m▂\u001b[39m▂\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m \u001b[39m▁\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▃\n",
+       "  \u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[34m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[32m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m▇\u001b[39m█\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▅\u001b[39m▅\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▅\u001b[39m▄\u001b[39m▄\u001b[39m \u001b[39m█\n",
+       "  1.77 μs\u001b[90m      \u001b[39m\u001b[90mHistogram: \u001b[39m\u001b[90m\u001b[1mlog(\u001b[22m\u001b[39m\u001b[90mfrequency\u001b[39m\u001b[90m\u001b[1m)\u001b[22m\u001b[39m\u001b[90m by time\u001b[39m      4.57 μs \u001b[0m\u001b[1m<\u001b[22m\n",
+       "\n",
+       " Memory estimate\u001b[90m: \u001b[39m\u001b[33m4.47 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m36\u001b[39m."
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "@benchmark solve(prob_ode_linear, Tsit5_for_relaxation())"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 16,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "BenchmarkTools.Trial: 10000 samples with 5 evaluations.\n",
+       " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m):  \u001b[39m\u001b[36m\u001b[1m5.920 μs\u001b[22m\u001b[39m … \u001b[35m366.900 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m0.00% … 95.12%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m):     \u001b[39m\u001b[34m\u001b[1m6.760 μs               \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m):    \u001b[39m0.00%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m):   \u001b[39m\u001b[32m\u001b[1m7.733 μs\u001b[22m\u001b[39m ± \u001b[32m 12.003 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m):  \u001b[39m5.79% ±  3.68%\n",
+       "\n",
+       "  \u001b[39m▁\u001b[39m▆\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[34m▆\u001b[39m\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[32m▅\u001b[39m\u001b[39m▆\u001b[39m▅\u001b[39m▄\u001b[39m▂\u001b[39m▁\u001b[39m▂\u001b[39m▁\u001b[39m \u001b[39m \u001b[39m \u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▂\u001b[39m▁\u001b[39m▂\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m▁\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▃\n",
+       "  \u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[34m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[32m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m▇\u001b[39m▆\u001b[39m▇\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▅\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▇\u001b[39m▆\u001b[39m▅\u001b[39m█\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▆\u001b[39m▆\u001b[39m \u001b[39m█\n",
+       "  5.92 μs\u001b[90m      \u001b[39m\u001b[90mHistogram: \u001b[39m\u001b[90m\u001b[1mlog(\u001b[22m\u001b[39m\u001b[90mfrequency\u001b[39m\u001b[90m\u001b[1m)\u001b[22m\u001b[39m\u001b[90m by time\u001b[39m      13.3 μs \u001b[0m\u001b[1m<\u001b[22m\n",
+       "\n",
+       " Memory estimate\u001b[90m: \u001b[39m\u001b[33m14.70 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m110\u001b[39m."
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "@benchmark solve(prob_ode_2Dlinear, Tsit5_for_relaxation())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "#### Compare time solving without and with relaxation steps\n",
+    "\n",
+    "We do the compareason on the non linear harmonic oscillator : "
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 22,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "Relaxation{DataType, var\"#21#22\"}(AlefeldPotraShi, var\"#21#22\"())"
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "f_nloscillator = (u, p, t) -> [-u[2]/(u[1]^2 + u[2]^2),u[1]/(u[1]^2 + u[2]^2)]\n",
+    "prob_nloscillator = ODEProblem(\n",
+    "    ODEFunction(f_nloscillator; analytic = (u0, p, t) -> [cos(t), sin(t)]),\n",
+    "    [1.0, 0.0],\n",
+    "    (0.0, 1.0))\n",
+    "r_nloscillator = Relaxation(AlefeldPotraShi, x-> norm(x))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "\n",
+    "\n",
+    "Without :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 26,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "BenchmarkTools.Trial: 10000 samples with 1 evaluation.\n",
+       " Range \u001b[90m(\u001b[39m\u001b[36m\u001b[1mmin\u001b[22m\u001b[39m … \u001b[35mmax\u001b[39m\u001b[90m):  \u001b[39m\u001b[36m\u001b[1m22.700 μs\u001b[22m\u001b[39m … \u001b[35m  4.386 ms\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmin … max\u001b[90m): \u001b[39m 0.00% … 97.69%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[34m\u001b[1mmedian\u001b[22m\u001b[39m\u001b[90m):     \u001b[39m\u001b[34m\u001b[1m26.700 μs               \u001b[22m\u001b[39m\u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmedian\u001b[90m):    \u001b[39m 0.00%\n",
+       " Time  \u001b[90m(\u001b[39m\u001b[32m\u001b[1mmean\u001b[22m\u001b[39m ± \u001b[32mσ\u001b[39m\u001b[90m):   \u001b[39m\u001b[32m\u001b[1m33.164 μs\u001b[22m\u001b[39m ± \u001b[32m116.550 μs\u001b[39m  \u001b[90m┊\u001b[39m GC \u001b[90m(\u001b[39mmean ± σ\u001b[90m):  \u001b[39m11.43% ±  3.25%\n",
+       "\n",
+       "  \u001b[39m▃\u001b[39m▇\u001b[39m█\u001b[39m▇\u001b[39m▆\u001b[39m▅\u001b[34m▅\u001b[39m\u001b[39m▅\u001b[39m▅\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▄\u001b[39m▂\u001b[39m▁\u001b[32m▁\u001b[39m\u001b[39m▂\u001b[39m▂\u001b[39m▁\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▂\u001b[39m▁\u001b[39m▁\u001b[39m \u001b[39m▁\u001b[39m \u001b[39m \u001b[39m▁\u001b[39m▁\u001b[39m▂\u001b[39m▁\u001b[39m▁\u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m \u001b[39m▂\n",
+       "  \u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[34m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[32m█\u001b[39m\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m█\u001b[39m▇\u001b[39m▇\u001b[39m▇\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▇\u001b[39m▅\u001b[39m▅\u001b[39m▅\u001b[39m▆\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▆\u001b[39m▆\u001b[39m▅\u001b[39m▆\u001b[39m▅\u001b[39m▅\u001b[39m▆\u001b[39m▅\u001b[39m▅\u001b[39m \u001b[39m█\n",
+       "  22.7 μs\u001b[90m       \u001b[39m\u001b[90mHistogram: \u001b[39m\u001b[90m\u001b[1mlog(\u001b[22m\u001b[39m\u001b[90mfrequency\u001b[39m\u001b[90m\u001b[1m)\u001b[22m\u001b[39m\u001b[90m by time\u001b[39m      63.9 μs \u001b[0m\u001b[1m<\u001b[22m\n",
+       "\n",
+       " Memory estimate\u001b[90m: \u001b[39m\u001b[33m50.98 KiB\u001b[39m, allocs estimate\u001b[90m: \u001b[39m\u001b[33m622\u001b[39m."
+      ]
+     },
+     "metadata": {},
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "@benchmark solve(prob_nloscillator, Tsit5())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "With :"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "@benchmark solve(prob_nloscillator, Tsit5_for_relaxation(); modif = r_nloscillator)"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.10.2",
+   "language": "julia",
+   "name": "julia-1.10"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.10.2"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}
diff --git a/test/relaxation/conserved_exponential_entropy.jl b/test/relaxation/conserved_exponential_entropy.jl
new file mode 100644
index 0000000000..6829b588af
--- /dev/null
+++ b/test/relaxation/conserved_exponential_entropy.jl
@@ -0,0 +1,45 @@
+using OrdinaryDiffEq, DiffEqDevTools
+
+include("relaxation.jl")
+
+printstyled("Conserved exponential entropy\n"; bold = true)
+
+dts = (1 / 2) .^ (6:-1:4)
+
+f = (u, p, t) -> [-exp(u[2]), exp(u[1])]
+prob  = ODEProblem(
+    ODEFunction(f ; 
+                analytic = (u0, p, t)->[log(exp(1) + exp(0.5)) - log(exp(0.5) + exp((exp(0.5)+exp(1))*t)), 
+                                        log(exp((exp(0.5)+exp(1))*t)*(exp(0.5)+exp(1)))/(exp(0.5) + exp((exp(0.5)+exp(1))*t))]),
+    [1.0, 0.5],
+    (0.0, 1.0))
+
+invariant(x) = exp(x[1]) + exp(x[2])
+
+# Convergence with the old method Tsit5()
+sim = test_convergence(dts, prob, Tsit5())
+println("order of convergence of older perform_step! : "*string(sim.𝒪est[:final]))
+
+# Convergence with the new method Tsit5_fors_relaxation() without relaxation
+sim = test_convergence(dts, prob, Tsit5_for_relaxation())
+println("order of convergence of new perform_step! without relaxation: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation without FSAL modification, i.e f(uₙ₊₁) ≈ f(uᵧ,ₙ₊₁), before EEst
+r = PerformStepCallback(;poststep = Relaxation(AlefeldPotraShi, invariant))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation without FSAL modification before EEst: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation without FSAL modification, i.e f(uᵧ,ₙ₊₁) ≈ f(uₙ₊₁) , after EEst
+r = PerformStepCallback(;postEEst = Relaxation(AlefeldPotraShi, invariant))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation without FSAL modification after EEst: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with FSAL-R, i.e  f(uᵧ,ₙ₊₁) ≈ f(uᵧ,ₙ) + γ ( f(uₙ₊₁) - f(uᵧ,ₙ)) 
+r = PerformStepCallback(;postEEst = Relaxation(AlefeldPotraShi, invariant, fsal_r))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation with FSAL-R modification: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with R-FSAL, i.e f(uₙ₊₁) ≈ f(uᵧ,ₙ) + 1/γ ( f(uᵧ,ₙ₊₁) - f(uᵧ,ₙ) )
+r = PerformStepCallback(;poststep = Relaxation(AlefeldPotraShi, invariant, r_fsal), postEEst = fsal_r_int)
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation with R-FSAL modification: "*string(sim.𝒪est[:final]))
\ No newline at end of file
diff --git a/test/relaxation/convergence.jl b/test/relaxation/convergence.jl
new file mode 100644
index 0000000000..192760eed6
--- /dev/null
+++ b/test/relaxation/convergence.jl
@@ -0,0 +1,59 @@
+using OrdinaryDiffEq, DiffEqDevTools,  Test
+
+import ODEProblemLibrary: prob_ode_linear,
+                          prob_ode_2Dlinear
+
+include("relaxation.jl")
+
+#########################################################################
+# Check that the code with the new structure without modification asked 
+# by the user gives the same result
+# Comparaison on Tsit5
+
+probnum = prob_ode_linear
+prob = prob_ode_2Dlinear
+
+
+testTol = 0.2
+
+### Tsit5()
+
+println("Tsit5")
+dts = (1 / 2) .^ (7:-1:3)
+sim = test_convergence(dts, probnum, Tsit5())
+@test abs.(sim.𝒪est[:l2] - 5) < testTol + 0.2
+sim = test_convergence(dts, prob, Tsit5())
+@test abs.(sim.𝒪est[:l2] - 5) < testTol + 0.2
+
+### Tsit5() with new structure
+
+println("Tsit5 with relaxation")
+dts = (1 / 2) .^ (7:-1:3)
+sim = test_convergence(dts, probnum, Tsit5_for_relaxation())
+@test abs.(sim.𝒪est[:l2] - 5) < testTol + 0.2
+sim = test_convergence(dts, prob, Tsit5_for_relaxation())  
+@test abs.(sim.𝒪est[:l2] - 5) < testTol + 0.2   
+
+
+#########################################################################
+##                            With Relaxation 
+
+########################################################
+# TEST  1 : Harmonic Oscillator
+include("harmonic_oscillator.jl")
+
+########################################################
+# TEST  2 : Non Linear Oscillator
+include("non_linear_oscillator.jl")
+
+########################################################
+# TEST  3 : Non Linear Pendulum
+include("non_linear_pendulum.jl")
+
+############################################################################
+# TEST  4 : Time dependent harmonic oscillator with bounded angular velocity
+include("time_dependant_harmonic_oscillator.jl")
+
+############################################################################
+# TEST  5 : Conserved exponential entropy
+include("conserved_exponential_entropy.jl")
\ No newline at end of file
diff --git a/test/relaxation/harmonic_oscillator.jl b/test/relaxation/harmonic_oscillator.jl
new file mode 100644
index 0000000000..f431019dcf
--- /dev/null
+++ b/test/relaxation/harmonic_oscillator.jl
@@ -0,0 +1,43 @@
+using OrdinaryDiffEq, DiffEqDevTools
+
+include("relaxation.jl")
+
+printstyled("Harmonic Oscillator\n"; bold = true)
+
+dts = (1 / 2) .^ (6:-1:4)
+
+f = (u, p, t) -> [-u[2],u[1]]
+prob = ODEProblem(
+    ODEFunction(f; analytic = (u0, p, t) -> [cos(t), sin(t)]),
+    [1.0, 0.0],
+    (0.0, 1.0))
+
+invariant(x) = norm(x)
+
+# Convergence with the old method Tsit5()
+sim = test_convergence(dts, prob, Tsit5())
+println("order of convergence of older perform_step! : "*string(sim.𝒪est[:final]))
+
+# Convergence with the new method Tsit5_fors_relaxation() without relaxation
+sim = test_convergence(dts, prob, Tsit5_for_relaxation())
+println("order of convergence of new perform_step! without relaxation: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation without FSAL modification, i.e f(uₙ₊₁) ≈ f(uᵧ,ₙ₊₁), before EEst
+r = PerformStepCallback(;poststep = Relaxation(AlefeldPotraShi, invariant))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation without FSAL modification before EEst: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation without FSAL modification, i.e f(uᵧ,ₙ₊₁) ≈ f(uₙ₊₁) , after EEst
+r = PerformStepCallback(;postEEst = Relaxation(AlefeldPotraShi, invariant))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation without FSAL modification after EEst: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with FSAL-R, i.e  f(uᵧ,ₙ₊₁) ≈ f(uᵧ,ₙ) + γ ( f(uₙ₊₁) - f(uᵧ,ₙ)) 
+r = PerformStepCallback(;postEEst = Relaxation(AlefeldPotraShi, invariant, fsal_r))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation with FSAL-R modification: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with R-FSAL, i.e f(uₙ₊₁) ≈ f(uᵧ,ₙ) + 1/γ ( f(uᵧ,ₙ₊₁) - f(uᵧ,ₙ) )
+r = PerformStepCallback(;poststep = Relaxation(AlefeldPotraShi, invariant, r_fsal), postEEst = fsal_r_int)
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation with R-FSAL modification: "*string(sim.𝒪est[:final]))
\ No newline at end of file
diff --git a/test/relaxation/non_linear_oscillator.jl b/test/relaxation/non_linear_oscillator.jl
new file mode 100644
index 0000000000..63b51f6128
--- /dev/null
+++ b/test/relaxation/non_linear_oscillator.jl
@@ -0,0 +1,43 @@
+using OrdinaryDiffEq, DiffEqDevTools
+
+include("relaxation.jl")
+
+printstyled("Non linear harmonic oscillator\n"; bold = true)
+
+dts = (1 / 2) .^ (6:-1:4)
+
+f = (u, p, t) -> [-u[2]/(u[1]^2 + u[2]^2),u[1]/(u[1]^2 + u[2]^2)]
+prob = ODEProblem(
+    ODEFunction(f; analytic = (u0, p, t) -> [cos(t), sin(t)]),
+    [1.0, 0.0],
+    (0.0, 1.0))
+
+invariant(x) = norm(x)
+
+# Convergence with the old method Tsit5()
+sim = test_convergence(dts, prob, Tsit5())
+println("order of convergence of older perform_step! : "*string(sim.𝒪est[:final]))
+
+# Convergence with the new method Tsit5_fors_relaxation() without relaxation
+sim = test_convergence(dts, prob, Tsit5_for_relaxation())
+println("order of convergence of new perform_step! without relaxation: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation without FSAL modification, i.e f(uₙ₊₁) ≈ f(uᵧ,ₙ₊₁), before EEst
+r = PerformStepCallback(;poststep = Relaxation(AlefeldPotraShi, invariant))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation without FSAL modification before EEst: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation without FSAL modification, i.e f(uᵧ,ₙ₊₁) ≈ f(uₙ₊₁) , after EEst
+r = PerformStepCallback(;postEEst = Relaxation(AlefeldPotraShi, invariant))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation without FSAL modification after EEst: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with FSAL-R, i.e  f(uᵧ,ₙ₊₁) ≈ f(uᵧ,ₙ) + γ ( f(uₙ₊₁) - f(uᵧ,ₙ)) 
+r = PerformStepCallback(;postEEst = Relaxation(AlefeldPotraShi, invariant, fsal_r))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation with FSAL-R modification: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with R-FSAL, i.e f(uₙ₊₁) ≈ f(uᵧ,ₙ) + 1/γ ( f(uᵧ,ₙ₊₁) - f(uᵧ,ₙ) )
+r = PerformStepCallback(;poststep = Relaxation(AlefeldPotraShi, invariant, r_fsal), postEEst = fsal_r_int)
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation with R-FSAL modification: "*string(sim.𝒪est[:final]))
\ No newline at end of file
diff --git a/test/relaxation/non_linear_pendulum.jl b/test/relaxation/non_linear_pendulum.jl
new file mode 100644
index 0000000000..15c61995dd
--- /dev/null
+++ b/test/relaxation/non_linear_pendulum.jl
@@ -0,0 +1,45 @@
+using OrdinaryDiffEq, DiffEqDevTools
+
+include("relaxation.jl")
+
+printstyled("Non Linear Pendulum\n"; bold = true)
+
+dts = (1 / 2) .^ (6:-1:4)
+
+f = (u, p, t) -> [-sin(u[2]), u[1]]
+prob = ODEProblem(
+    f,
+    [1.0, 0.0],
+    (0.0, 1.0))
+
+invariant(x) = x[1]^2/2 -  cos(x[2])
+
+test_setup = Dict(:alg => Vern9(), :reltol => 1e-14, :abstol => 1e-14)
+
+# Convergence with the old method Tsit5()
+sim = analyticless_test_convergence(dts, prob, Tsit5(), test_setup)
+println("order of convergence of older perform_step! : "*string(sim.𝒪est[:final]))
+
+# Convergence with the new method Tsit5_fors_relaxation() without relaxation
+sim = analyticless_test_convergence(dts, prob, Tsit5_for_relaxation(), test_setup)
+println("order of convergence of new perform_step! without relaxation: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation without FSAL modification, i.e f(uₙ₊₁) ≈ f(uᵧ,ₙ₊₁), before EEst
+r = PerformStepCallback(;poststep = Relaxation(AlefeldPotraShi, invariant))
+sim = analyticless_test_convergence(dts, prob, Tsit5_for_relaxation(), test_setup; modif = r)
+println("order with relaxation without FSAL modification before EEst: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation without FSAL modification, i.e f(uᵧ,ₙ₊₁) ≈ f(uₙ₊₁) , after EEst
+r = PerformStepCallback(;postEEst = Relaxation(AlefeldPotraShi, invariant))
+sim = analyticless_test_convergence(dts, prob, Tsit5_for_relaxation(), test_setup; modif = r)
+println("order with relaxation without FSAL modification after EEst: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with FSAL-R, i.e  f(uᵧ,ₙ₊₁) ≈ f(uᵧ,ₙ) + γ ( f(uₙ₊₁) - f(uᵧ,ₙ)) 
+r = PerformStepCallback(;postEEst = Relaxation(AlefeldPotraShi, invariant, fsal_r))
+sim = analyticless_test_convergence(dts, prob, Tsit5_for_relaxation(), test_setup; modif = r)
+println("order with relaxation with FSAL-R modification: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with R-FSAL, i.e f(uₙ₊₁) ≈ f(uᵧ,ₙ) + 1/γ ( f(uᵧ,ₙ₊₁) - f(uᵧ,ₙ) )
+r = PerformStepCallback(;poststep = Relaxation(AlefeldPotraShi, invariant, r_fsal), postEEst = fsal_r_int)
+sim = analyticless_test_convergence(dts, prob, Tsit5_for_relaxation(), test_setup; modif = r)
+println("order with relaxation with R-FSAL modification: "*string(sim.𝒪est[:final]))
\ No newline at end of file
diff --git a/test/relaxation/relaxation.jl b/test/relaxation/relaxation.jl
new file mode 100644
index 0000000000..41c1a3c725
--- /dev/null
+++ b/test/relaxation/relaxation.jl
@@ -0,0 +1,45 @@
+using Roots
+using LinearAlgebra
+using UnPack
+
+mutable struct Relaxation{OPT, INV, FSAL}
+    opt::OPT
+    invariant::INV
+    has_fsal::Bool
+    update_fsal::FSAL
+    γ
+    Relaxation(opt, inv, fsal = nothing) = new{typeof(opt), typeof(inv), typeof(fsal)}(opt, inv, fsal isa Nothing ? false : true, fsal, 1.0)
+end
+
+function (r::Relaxation)(integrator)
+
+    @unpack t, dt, uprev, u = integrator
+
+    # We fix here the bounds of interval where we are going to look for the relaxation
+    #(gamma_min, gamma_max) = apriori_bounds_dt(integrator) ./ dt
+    
+    S_u = u - uprev
+
+    # Find relaxation paramter gamma
+    gamma_min = 0.5
+    gamma_max = 1.5
+
+    if (r.invariant(gamma_min*S_u .+ uprev) .- r.invariant(uprev)) * (r.invariant(gamma_max*S_u .+ uprev) .- r.invariant(uprev)) ≤ 0
+        gamma = find_zero(gamma -> r.invariant(gamma*S_u .+ uprev) .- r.invariant(uprev),
+                            (gamma_min, gamma_max),
+                            r.opt())
+        r.γ = gamma
+        change_dt!(integrator, gamma*dt)
+        change_u!(integrator, uprev + gamma*S_u)
+
+        if r.has_fsal
+            integrator.fsallast = r.update_fsal(gamma, integrator.fsalfirst, integrator.fsallast)
+        end
+    end
+end
+
+
+fsal_r(gamma, fsalfirst, fsallast) = fsalfirst + gamma * (fsallast - fsalfirst)
+r_fsal(gamma, fsalfirst, fsallast) = fsalfirst + 1/gamma * (fsallast - fsalfirst)
+fsal_r_int(integrator) = integrator.fsallast = fsal_r(integrator.opts.performstepcallback.poststep_cb.γ, integrator.fsalfirst, integrator.fsallast)
+
diff --git a/test/relaxation/time_dependant_harmonic_oscillator.jl b/test/relaxation/time_dependant_harmonic_oscillator.jl
new file mode 100644
index 0000000000..afa2b7bab1
--- /dev/null
+++ b/test/relaxation/time_dependant_harmonic_oscillator.jl
@@ -0,0 +1,45 @@
+using OrdinaryDiffEq, DiffEqDevTools
+
+include("relaxation.jl")
+
+printstyled("Time-dependent harmonic oscillator\n"; bold = true)
+
+dts = (1 / 2) .^ (6:-1:4)
+
+f = (u, p, t) -> [-(1+0.5 * sin(t))*u[2], (1+0.5 * sin(t))*u[1]]
+prob = ODEProblem(
+    ODEFunction(f; 
+                analytic = (u0, p, t)->[cos(0.5)*cos(t-0.5*cos(t))-sin(0.5)*sin(t-0.5*cos(t)), 
+                                        sin(0.5)*cos(t-0.5*cos(t))+cos(0.5)*sin(t-0.5*cos(t))]),
+    [1.0, 0.0],
+    (0.0, 1.0))
+
+invariant(x) = norm(x)
+
+# Convergence with the old method Tsit5()
+sim = test_convergence(dts, prob, Tsit5())
+println("order of convergence of older perform_step! : "*string(sim.𝒪est[:final]))
+
+# Convergence with the new method Tsit5_fors_relaxation() without relaxation
+sim = test_convergence(dts, prob, Tsit5_for_relaxation())
+println("order of convergence of new perform_step! without relaxation: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation without FSAL modification, i.e f(uₙ₊₁) ≈ f(uᵧ,ₙ₊₁), before EEst
+r = PerformStepCallback(;poststep = Relaxation(AlefeldPotraShi, invariant))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation without FSAL modification before EEst: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation without FSAL modification, i.e f(uᵧ,ₙ₊₁) ≈ f(uₙ₊₁) , after EEst
+r = PerformStepCallback(;postEEst = Relaxation(AlefeldPotraShi, invariant))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation without FSAL modification after EEst: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with FSAL-R, i.e  f(uᵧ,ₙ₊₁) ≈ f(uᵧ,ₙ) + γ ( f(uₙ₊₁) - f(uᵧ,ₙ)) 
+r = PerformStepCallback(;postEEst = Relaxation(AlefeldPotraShi, invariant, fsal_r))
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation with FSAL-R modification: "*string(sim.𝒪est[:final]))
+
+# Convergence with relaxation with R-FSAL, i.e f(uₙ₊₁) ≈ f(uᵧ,ₙ) + 1/γ ( f(uᵧ,ₙ₊₁) - f(uᵧ,ₙ) )
+r = PerformStepCallback(;poststep = Relaxation(AlefeldPotraShi, invariant, r_fsal), postEEst = fsal_r_int)
+sim = test_convergence(dts, prob, Tsit5_for_relaxation(); modif = r)
+println("order with relaxation with R-FSAL modification: "*string(sim.𝒪est[:final]))
\ No newline at end of file