diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl
index 88dff7f11..32f3c559b 100644
--- a/src/NonlinearSolve.jl
+++ b/src/NonlinearSolve.jl
@@ -70,6 +70,7 @@ include("core/spectral_methods.jl")
 
 include("algorithms/raphson.jl")
 include("algorithms/pseudo_transient.jl")
+include("algorithms/multistep.jl")
 include("algorithms/broyden.jl")
 include("algorithms/klement.jl")
 include("algorithms/lbroyden.jl")
@@ -131,7 +132,8 @@ include("default.jl")
 end
 
 # Core Algorithms
-export NewtonRaphson, PseudoTransient, Klement, Broyden, LimitedMemoryBroyden, DFSane
+export NewtonRaphson, PseudoTransient, Klement, Broyden, LimitedMemoryBroyden, DFSane,
+    MultiStepNonlinearSolver
 export GaussNewton, LevenbergMarquardt, TrustRegion
 export NonlinearSolvePolyAlgorithm,
     RobustMultiNewton, FastShortcutNonlinearPolyalg, FastShortcutNLLSPolyalg
@@ -145,7 +147,9 @@ export GeneralizedFirstOrderAlgorithm, ApproximateJacobianSolveAlgorithm, Genera
 
 # Descent Algorithms
 export NewtonDescent, SteepestDescent, Dogleg, DampedNewtonDescent,
-    GeodesicAcceleration
+    GeodesicAcceleration, GenericMultiStepDescent
+## Multistep Algorithms
+export MultiStepSchemes
 
 # Globalization
 ## Line Search Algorithms
diff --git a/src/algorithms/multistep.jl b/src/algorithms/multistep.jl
new file mode 100644
index 000000000..35b204094
--- /dev/null
+++ b/src/algorithms/multistep.jl
@@ -0,0 +1,7 @@
+function MultiStepNonlinearSolver(; concrete_jac = nothing, linsolve = nothing,
+        scheme = MSS.PotraPtak3, precs = DEFAULT_PRECS, autodiff = nothing)
+    descent = GenericMultiStepDescent(; scheme, linsolve, precs)
+    # TODO: Use the scheme as the name
+    return GeneralizedFirstOrderAlgorithm(; concrete_jac, name = :MultiStepNonlinearSolver,
+        descent, jacobian_ad = autodiff)
+end
diff --git a/src/descent/common.jl b/src/descent/common.jl
index 10b14ad14..2a614d84a 100644
--- a/src/descent/common.jl
+++ b/src/descent/common.jl
@@ -5,11 +5,11 @@ Construct a `DescentResult` object.
 
 ### Keyword Arguments
 
-  * `δu`: The descent direction.
-  * `u`: The new iterate. This is provided only for multi-step methods currently.
-  * `success`: Certain Descent Algorithms can reject a descent direction for example
+  - `δu`: The descent direction.
+  - `u`: The new iterate. This is provided only for multi-step methods currently.
+  - `success`: Certain Descent Algorithms can reject a descent direction for example
     [`GeodesicAcceleration`](@ref).
-  * `extras`: A named tuple containing intermediates computed during the solve.
+  - `extras`: A named tuple containing intermediates computed during the solve.
     For example, [`GeodesicAcceleration`](@ref) returns `NamedTuple{(:v, :a)}` containing
     the "velocity" and "acceleration" terms.
 """
diff --git a/src/descent/multistep.jl b/src/descent/multistep.jl
index e69de29bb..2879a9bef 100644
--- a/src/descent/multistep.jl
+++ b/src/descent/multistep.jl
@@ -0,0 +1,134 @@
+"""
+    MultiStepSchemes
+
+This module defines the multistep schemes used in the multistep descent algorithms. The
+naming convention follows <name of method><order of convergence>. The name of method is
+typically the last names of the authors of the paper that introduced the method.
+"""
+module MultiStepSchemes
+
+abstract type AbstractMultiStepScheme end
+
+function Base.show(io::IO, mss::AbstractMultiStepScheme)
+    print(io, "MultiStepSchemes.$(string(nameof(typeof(mss)))[3:end])")
+end
+
+struct __PotraPtak3 <: AbstractMultiStepScheme end
+const PotraPtak3 = __PotraPtak3()
+
+alg_steps(::__PotraPtak3) = 1
+
+struct __SinghSharma4 <: AbstractMultiStepScheme end
+const SinghSharma4 = __SinghSharma4()
+
+alg_steps(::__SinghSharma4) = 3
+
+struct __SinghSharma5 <: AbstractMultiStepScheme end
+const SinghSharma5 = __SinghSharma5()
+
+alg_steps(::__SinghSharma5) = 3
+
+struct __SinghSharma7 <: AbstractMultiStepScheme end
+const SinghSharma7 = __SinghSharma7()
+
+alg_steps(::__SinghSharma7) = 4
+
+end
+
+const MSS = MultiStepSchemes
+
+@kwdef @concrete struct GenericMultiStepDescent <: AbstractDescentAlgorithm
+    scheme
+    linsolve = nothing
+    precs = DEFAULT_PRECS
+end
+
+supports_line_search(::GenericMultiStepDescent) = false
+supports_trust_region(::GenericMultiStepDescent) = false
+
+@concrete mutable struct GenericMultiStepDescentCache{S, INV} <: AbstractDescentCache
+    f
+    p
+    δu
+    δus
+    scheme::S
+    lincache
+    timer
+    nf::Int
+end
+
+@internal_caches GenericMultiStepDescentCache :lincache
+
+function __reinit_internal!(cache::GenericMultiStepDescentCache, args...; p = cache.p,
+        kwargs...)
+    cache.nf = 0
+    cache.p = p
+end
+
+function __δu_caches(scheme::MSS.__PotraPtak3, fu, u, ::Val{N}) where {N}
+    caches = ntuple(N) do i
+        @bb δu = similar(u)
+        @bb y = similar(u)
+        @bb fy = similar(fu)
+        @bb δy = similar(u)
+        @bb u_new = similar(u)
+        (δu, δy, fy, y, u_new)
+    end
+    return first(caches), (N ≤ 1 ? nothing : caches[2:end])
+end
+
+function __internal_init(prob::NonlinearProblem, alg::GenericMultiStepDescent, J, fu, u;
+        shared::Val{N} = Val(1), pre_inverted::Val{INV} = False, linsolve_kwargs = (;),
+        abstol = nothing, reltol = nothing, timer = get_timer_output(),
+        kwargs...) where {INV, N}
+    δu, δus = __δu_caches(alg.scheme, fu, u, shared)
+    INV && return GenericMultiStepDescentCache{true}(prob.f, prob.p, δu, δus,
+        alg.scheme, nothing, timer, 0)
+    lincache = LinearSolverCache(alg, alg.linsolve, J, _vec(fu), _vec(u); abstol, reltol,
+        linsolve_kwargs...)
+    return GenericMultiStepDescentCache{false}(prob.f, prob.p, δu, δus, alg.scheme,
+        lincache, timer, 0)
+end
+
+function __internal_init(prob::NonlinearLeastSquaresProblem, alg::GenericMultiStepDescent,
+        J, fu, u; kwargs...)
+    error("Multi-Step Descent Algorithms for NLLS are not implemented yet.")
+end
+
+function __internal_solve!(cache::GenericMultiStepDescentCache{MSS.__PotraPtak3, INV}, J,
+        fu, u, idx::Val = Val(1); skip_solve::Bool = false, new_jacobian::Bool = true,
+        kwargs...) where {INV}
+    (u_new, δy, fy, y, δu) = get_du(cache, idx)
+    skip_solve && return DescentResult(; u = u_new)
+
+    @static_timeit cache.timer "linear solve" begin
+        @static_timeit cache.timer "solve and step 1" begin
+            if INV
+                J !== nothing && @bb(δu=J × _vec(fu))
+            else
+                δu = cache.lincache(; A = J, b = _vec(fu), kwargs..., linu = _vec(δu),
+                    du = _vec(δu),
+                    reuse_A_if_factorization = !new_jacobian || (idx !== Val(1)))
+                δu = _restructure(u, δu)
+            end
+            @bb @. y = u - δu
+        end
+
+        fy = evaluate_f!!(cache.f, fy, y, cache.p)
+        cache.nf += 1
+
+        @static_timeit cache.timer "solve and step 2" begin
+            if INV
+                J !== nothing && @bb(δy=J × _vec(fy))
+            else
+                δy = cache.lincache(; A = J, b = _vec(fy), kwargs..., linu = _vec(δy),
+                    du = _vec(δy), reuse_A_if_factorization = true)
+                δy = _restructure(u, δy)
+            end
+            @bb @. u_new = y - δy
+        end
+    end
+
+    set_du!(cache, (u_new, δy, fy, y, δu), idx)
+    return DescentResult(; u = u_new)
+end
diff --git a/src/internal/tracing.jl b/src/internal/tracing.jl
index 667c6ce07..bfb93c6d7 100644
--- a/src/internal/tracing.jl
+++ b/src/internal/tracing.jl
@@ -187,6 +187,7 @@ function update_trace!(cache::AbstractNonlinearSolveCache, α = true)
     trace === nothing && return nothing
 
     J = __getproperty(cache, Val(:J))
+    # TODO: fix tracing for multi-step methods where du is not aliased properly
     if J === nothing
         update_trace!(trace, get_nsteps(cache) + 1, get_u(cache), get_fu(cache),
             nothing, cache.du, α)