diff --git a/lib/ImplicitDiscreteSolve/Project.toml b/lib/ImplicitDiscreteSolve/Project.toml
index d25e9f75d1..f609fb5bac 100644
--- a/lib/ImplicitDiscreteSolve/Project.toml
+++ b/lib/ImplicitDiscreteSolve/Project.toml
@@ -4,38 +4,40 @@ authors = ["vyudu <vincent.duyuan@gmail.com>"]
 version = "1.2.0"
 
 [deps]
-SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
-SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-UnPack = "3a884ed6-31ef-47d7-9d2a-63182c4928ed"
+ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
-Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
+NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
+NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
 OrdinaryDiffEqCore = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
+Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
+SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
 
-[extras]
-JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
-OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
-Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
-Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
-AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a"
+[sources]
+OrdinaryDiffEqCore = {path = "../OrdinaryDiffEqCore"}
 
 [compat]
-Test = "1.10.0"
+AllocCheck = "0.2"
+Aqua = "0.8.11"
+ConcreteStructs = "0.2.3"
+DiffEqBase = "6.176"
+JET = "0.9.18, 0.10.4"
+NonlinearSolveBase = "2.0.0"
+NonlinearSolveFirstOrder = "1.9.0"
+OrdinaryDiffEqCore = "1.29.0"
 OrdinaryDiffEqSDIRK = "1.6.0"
+Reexport = "1.2"
 SciMLBase = "2.99"
-SimpleNonlinearSolve = "2.7"
-OrdinaryDiffEqCore = "1.29.0"
-Aqua = "0.8.11"
 SymbolicIndexingInterface = "0.3.38"
+Test = "1.10.0"
 julia = "1.10"
-JET = "0.9.18, 0.10.4"
-UnPack = "1.0.2"
-AllocCheck = "0.2"
-DiffEqBase = "6.176"
-Reexport = "1.2"
+
+[extras]
+AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a"
+Aqua = "4c88cf16-eb10-579e-8560-4a9242c79595"
+JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
+OrdinaryDiffEqSDIRK = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
 test = ["OrdinaryDiffEqSDIRK", "Test", "JET", "Aqua", "AllocCheck"]
-
-[sources.OrdinaryDiffEqCore]
-path = "../OrdinaryDiffEqCore"
diff --git a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
index b52829f5ed..60d05501aa 100644
--- a/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
+++ b/lib/ImplicitDiscreteSolve/src/ImplicitDiscreteSolve.jl
@@ -1,28 +1,25 @@
 module ImplicitDiscreteSolve
 
 using SciMLBase
-using SimpleNonlinearSolve
-using UnPack
+using NonlinearSolveBase
+using NonlinearSolveFirstOrder
+using ConcreteStructs
 using SymbolicIndexingInterface: parameter_symbols
 import OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, alg_cache, OrdinaryDiffEqMutableCache,
                            OrdinaryDiffEqConstantCache, get_fsalfirstlast, isfsal,
                            initialize!, perform_step!, isdiscretecache, isdiscretealg,
                            alg_order, beta2_default, beta1_default, dt_required,
-                           _initialize_dae!, DefaultInit, BrownFullBasicInit, OverrideInit
+                           _initialize_dae!, DefaultInit, BrownFullBasicInit, OverrideInit,
+                           @muladd, @.., _unwrap_val, OrdinaryDiffEqCore, isadaptive
 
 using Reexport
 @reexport using SciMLBase
 
-"""
-    IDSolve()
-
-Simple solver for `ImplicitDiscreteSystems`. Uses `SimpleNewtonRaphson` to solve for the next state at every timestep.
-"""
-struct IDSolve <: OrdinaryDiffEqAlgorithm end
-
+include("algorithms.jl")
 include("cache.jl")
 include("solve.jl")
 include("alg_utils.jl")
+include("controller.jl")
 
 export IDSolve
 
diff --git a/lib/ImplicitDiscreteSolve/src/alg_utils.jl b/lib/ImplicitDiscreteSolve/src/alg_utils.jl
index 5c25c2c217..a2d6e6ecb5 100644
--- a/lib/ImplicitDiscreteSolve/src/alg_utils.jl
+++ b/lib/ImplicitDiscreteSolve/src/alg_utils.jl
@@ -11,9 +11,33 @@ end
 SciMLBase.isdiscrete(alg::IDSolve) = true
 
 isfsal(alg::IDSolve) = false
-alg_order(alg::IDSolve) = 0
+alg_order(alg::IDSolve) = 1
 beta2_default(alg::IDSolve) = 0
 beta1_default(alg::IDSolve, beta2) = 0
 
 dt_required(alg::IDSolve) = false
 isdiscretealg(alg::IDSolve) = true
+
+isadaptive(alg::IDSolve) = true
+
+# @concrete struct ConvergenceRateTracing
+#     inner_tracing
+# end
+
+# @concrete struct ConvergenceRateTraceTrick
+#     incrementL2norms
+#     residualL2norms
+#     trace_wrapper
+# end
+
+# function NonlinearSolveBase.init_nonlinearsolve_trace(
+#         prob, alg::IDSolve, u, fu, J, δu;
+#         trace_level::ConvergenceRateTracing, kwargs... # This kind of dispatch does not work. Need to figure out a different way.
+# )
+#     inner_trace = NonlinearSolveBase.init_nonlinearsolve_trace(
+#         prob, alg, u, fu, J, δu;
+#         trace_level.inner_tracing, kwargs...
+#     )
+
+#     return ConvergenceRateTraceTrick(eltype(δu)[], eltype(fu)[], inner_trace)
+# end
diff --git a/lib/ImplicitDiscreteSolve/src/algorithms.jl b/lib/ImplicitDiscreteSolve/src/algorithms.jl
new file mode 100644
index 0000000000..aa2c567598
--- /dev/null
+++ b/lib/ImplicitDiscreteSolve/src/algorithms.jl
@@ -0,0 +1,19 @@
+"""
+    IDSolve()
+
+First order solver for `ImplicitDiscreteSystems`.
+"""
+struct IDSolve{NLS} <:
+       OrdinaryDiffEqAlgorithm
+    nlsolve::NLS
+    extrapolant::Symbol
+    controller::Symbol
+end
+
+function IDSolve(;
+        nlsolve = NewtonRaphson(),
+        extrapolant = :constant,
+        controller = :PI
+)
+    IDSolve{typeof(nlsolve)}(nlsolve, extrapolant, controller)
+end
diff --git a/lib/ImplicitDiscreteSolve/src/cache.jl b/lib/ImplicitDiscreteSolve/src/cache.jl
index f9c7835d82..f3200762d4 100644
--- a/lib/ImplicitDiscreteSolve/src/cache.jl
+++ b/lib/ImplicitDiscreteSolve/src/cache.jl
@@ -1,14 +1,15 @@
-mutable struct ImplicitDiscreteState{uType, pType, tType}
+struct ImplicitDiscreteState{uType, pType, tType}
     u::uType
     p::pType
-    t_next::tType
+    t::tType
 end
 
-mutable struct IDSolveCache{uType} <: OrdinaryDiffEqMutableCache
+mutable struct IDSolveCache{uType, cType, thetaType} <: OrdinaryDiffEqMutableCache
     u::uType
     uprev::uType
-    state::ImplicitDiscreteState
-    prob::Union{Nothing, SciMLBase.AbstractNonlinearProblem}
+    z::uType
+    nlcache::cType
+    Θks::thetaType
 end
 
 function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits},
@@ -16,21 +17,50 @@ function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits},
         dt, reltol, p, calck,
         ::Val{true}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
     state = ImplicitDiscreteState(isnothing(u) ? nothing : zero(u), p, t)
-    IDSolveCache(u, uprev, state, nothing)
-end
+    f_nl = (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t)
 
-isdiscretecache(cache::IDSolveCache) = true
+    u_len = isnothing(u) ? 0 : length(u)
+    nlls = !isnothing(f.resid_prototype) && (length(f.resid_prototype) != u_len)
+    unl = isnothing(u) ? Float64[] : u # FIXME nonlinear solve cannot handle nothing for u
+    prob = if nlls
+        NonlinearLeastSquaresProblem{isinplace(f)}(
+            NonlinearFunction(f_nl; resid_prototype = f.resid_prototype),
+            unl, state)
+    else
+        NonlinearProblem{isinplace(f)}(f_nl, unl, state)
+    end
+
+    nlcache = init(prob, alg.nlsolve)
 
-struct IDSolveConstantCache <: OrdinaryDiffEqConstantCache
-    prob::Union{Nothing, SciMLBase.AbstractNonlinearProblem}
+    IDSolveCache(u, uprev, state.u, nlcache, uBottomEltypeNoUnits[])
 end
 
+isdiscretecache(cache::IDSolveCache) = true
+
 function alg_cache(alg::IDSolve, u, rate_prototype, ::Type{uEltypeNoUnits},
         ::Type{uBottomEltypeNoUnits}, ::Type{tTypeNoUnits}, uprev, uprev2, f, t,
         dt, reltol, p, calck,
         ::Val{false}) where {uEltypeNoUnits, uBottomEltypeNoUnits, tTypeNoUnits}
+    @assert !isnothing(u) "Empty u not supported with out of place functions yet."
+
     state = ImplicitDiscreteState(isnothing(u) ? nothing : zero(u), p, t)
-    IDSolveCache(u, uprev, state, nothing)
+    f_nl = (u_next, p) -> f(u_next, p.u, p.p, p.t)
+
+    u_len = isnothing(u) ? 0 : length(u)
+    nlls = !isnothing(f.resid_prototype) && (length(f.resid_prototype) != u_len)
+    unl = isnothing(u) ? Float64[] : u # FIXME nonlinear solve cannot handle nothing for u
+    prob = if nlls
+        NonlinearLeastSquaresProblem{isinplace(f)}(
+            NonlinearFunction(f_nl; resid_prototype = f.resid_prototype),
+            unl, state)
+    else
+        NonlinearProblem{isinplace(f)}(f_nl, unl, state)
+    end
+
+    nlcache = init(prob, alg.nlsolve)
+
+    # FIXME Use IDSolveConstantCache?
+    IDSolveCache(u, uprev, state.u, nlcache, uBottomEltypeNoUnits[])
 end
 
 get_fsalfirstlast(cache::IDSolveCache, rate_prototype) = (nothing, nothing)
diff --git a/lib/ImplicitDiscreteSolve/src/controller.jl b/lib/ImplicitDiscreteSolve/src/controller.jl
new file mode 100644
index 0000000000..63a27c9f84
--- /dev/null
+++ b/lib/ImplicitDiscreteSolve/src/controller.jl
@@ -0,0 +1,52 @@
+Base.@kwdef struct KantorovichTypeController <: OrdinaryDiffEqCore.AbstractController
+    Θmin::Float64
+    p::Int64
+    Θreject::Float64 = 0.95
+    Θbar::Float64 = 0.5
+    γ::Float64 = 0.95
+    qmin::Float64 = 1 / 5
+    qmax::Float64 = 5.0
+end
+
+function OrdinaryDiffEqCore.default_controller(
+        alg::IDSolve, cache::IDSolveCache, _1, _2, _3)
+    return KantorovichTypeController(; Θmin = 1 // 8, p = 1)
+end
+
+function OrdinaryDiffEqCore.stepsize_controller!(
+        integrator, controller::KantorovichTypeController, alg::IDSolve
+)
+    @inline g(x) = √(1 + 4x) - 1
+
+    # Adapt dt with a priori estimate (Eq. 5.24)
+    (; Θks) = integrator.cache
+    (; Θbar, γ, Θmin, qmin, qmax, p) = controller
+
+    Θ₀ = length(Θks) > 0 ? max(first(Θks), Θmin) : Θmin
+    q = clamp(γ * (g(Θbar) / (g(Θ₀)))^(1 / p), qmin, qmax)
+
+    return q
+end
+
+function OrdinaryDiffEqCore.step_accept_controller!(
+        integrator, controller::KantorovichTypeController, alg::IDSolve, q
+)
+    return q * integrator.dt
+end
+
+function OrdinaryDiffEqCore.step_reject_controller!(
+        integrator, controller::KantorovichTypeController, alg::IDSolve
+)
+    @inline g(x) = √(1 + 4x) - 1
+
+    # Shorten dt according to (Eq. 5.24)
+    (; Θks) = integrator.cache.nlcache
+    (; Θbar, Θreject, γ, Θmin, qmin, qmax, p) = controller
+    for Θk in Θks
+        if Θk > Θreject
+            q = clamp(γ * (g(Θbar) / g(Θk))^(1 / p), qmin, qmax)
+            integrator.dt = q * integrator.dt
+            return
+        end
+    end
+end
diff --git a/lib/ImplicitDiscreteSolve/src/solve.jl b/lib/ImplicitDiscreteSolve/src/solve.jl
index e0580799d8..d5a511fbab 100644
--- a/lib/ImplicitDiscreteSolve/src/solve.jl
+++ b/lib/ImplicitDiscreteSolve/src/solve.jl
@@ -1,38 +1,72 @@
-# Remake the nonlinear problem, then update
 function perform_step!(integrator, cache::IDSolveCache, repeat_step = false)
-    @unpack alg, u, uprev, dt, t, f, p = integrator
-    @unpack state, prob = cache
-    state.u .= uprev
-    state.t_next = t
-    prob = remake(prob, p = state)
+    (; alg, u, uprev, dt, t, tprev, f, p) = integrator
+    (; nlcache, Θks) = cache
 
-    u = solve(prob, SimpleNewtonRaphson())
-    integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, u.retcode)
-    integrator.u = u
-end
+    # initial guess
+    if alg.extrapolant == :linear
+        @.. broadcast=false cache.z=integrator.u + dt * (integrator.u - integrator.uprev)
+    else # :constant
+        cache.z .= integrator.u
+    end
+    state = ImplicitDiscreteState(cache.z, p, t + dt)
 
-function initialize!(integrator, cache::IDSolveCache)
-    integrator.u isa AbstractVector && (cache.state.u .= integrator.u)
-    cache.state.p = integrator.p
-    cache.state.t_next = integrator.t
-    f = integrator.f
+    # nonlinear solve step
+    SciMLBase.reinit!(nlcache, p = state)
 
-    _f = if isinplace(f)
-        (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t_next)
-    else
-        (u_next, p) -> f(u_next, p.u, p.p, p.t_next)
+    # solve!(nlcache)
+    # The solve here is simply unrolled by hand to query the convergence rate estimates "manually" for now
+    if nlcache.retcode == ReturnCode.InitialFailure
+        integrator.force_stepfail = true
+        return
     end
-    u_len = isnothing(integrator.u) ? 0 : length(integrator.u)
-    nlls = !isnothing(f.resid_prototype) && (length(f.resid_prototype) != u_len)
 
-    prob = if nlls
-        NonlinearLeastSquaresProblem{isinplace(f)}(
-            NonlinearFunction(_f; resid_prototype = f.resid_prototype),
-            cache.state.u, cache.state)
-    else
-        NonlinearProblem{isinplace(f)}(_f, cache.state.u, cache.state)
+    resize!(Θks, 0)
+    residualnormprev = zero(eltype(u))
+    while NonlinearSolveBase.not_terminated(nlcache)
+        step!(nlcache)
+        residualnorm = NonlinearSolveBase.L2_NORM(nlcache.fu)
+        if nlcache.nsteps > 1
+            # Θk = min(residualnorm/residualnormprev, incrementnorm/incrementnormprev)
+            Θk = residualnorm / residualnormprev
+            if residualnormprev ≈ 0.0 #|| incrementnormprev ≈ 0.0
+                push!(Θks, 0.0)
+            else
+                push!(Θks, Θk)
+            end
+            # if nlcache.parameters.enforce_monotonic_convergence && Θk ≥ 1.0
+            #     @debug "Newton-Raphson diverged. Aborting. ||r|| = $residualnorm" _group=:nlsolve
+            #     return false
+            # end
+        end
+        residualnormprev = residualnorm
+    end
+
+    # The solver might have set a different `retcode`
+    if nlcache.retcode == ReturnCode.Default
+        nlcache.retcode = ifelse(
+            nlcache.nsteps ≥ nlcache.maxiters, ReturnCode.MaxIters, ReturnCode.Success
+        )
     end
-    cache.prob = prob
+
+    NonlinearSolveBase.update_from_termination_cache!(nlcache.termination_cache, nlcache)
+
+    NonlinearSolveBase.update_trace!(
+        nlcache.trace, nlcache.nsteps, NonlinearSolveBase.get_u(nlcache),
+        NonlinearSolveBase.get_fu(nlcache), nothing, nothing, nothing;
+        last = Val(true)
+    )
+
+    if nlcache.retcode != ReturnCode.Success
+        integrator.force_stepfail = true
+        return
+    end
+
+    # Accept step
+    u .= nlcache.u
+end
+
+function initialize!(integrator, cache::IDSolveCache)
+    integrator.u isa AbstractVector && (cache.z .= integrator.u)
 end
 
 function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem,
@@ -43,13 +77,14 @@ function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem,
         _initialize_dae!(integrator, prob,
             OverrideInit(atol), x)
     else
-        @unpack u, p, t, f = integrator
+        (; u, p, t, f) = integrator
+
         initstate = ImplicitDiscreteState(u, p, t)
 
         _f = if isinplace(f)
-            (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t_next)
+            (resid, u_next, p) -> f(resid, u_next, p.u, p.p, p.t)
         else
-            (u_next, p) -> f(u_next, p.u, p.p, p.t_next)
+            (u_next, p) -> f(u_next, p.u, p.p, p.t)
         end
 
         nlls = !isnothing(f.resid_prototype) &&
@@ -60,11 +95,12 @@ function _initialize_dae!(integrator, prob::ImplicitDiscreteProblem,
         else
             NonlinearProblem{isinplace(f)}(_f, u, initstate)
         end
-        sol = solve(prob, SimpleNewtonRaphson())
+        sol = solve(prob, integrator.alg.nlsolve)
         if sol.retcode == ReturnCode.Success
             integrator.u = sol
         else
-            integrator.sol = SciMLBase.solution_new_retcode(integrator.sol, ReturnCode.InitialFailure)
+            integrator.sol = SciMLBase.solution_new_retcode(
+                integrator.sol, ReturnCode.InitialFailure)
         end
     end
 end
diff --git a/lib/ImplicitDiscreteSolve/test/runtests.jl b/lib/ImplicitDiscreteSolve/test/runtests.jl
index 4c94621302..7837874636 100644
--- a/lib/ImplicitDiscreteSolve/test/runtests.jl
+++ b/lib/ImplicitDiscreteSolve/test/runtests.jl
@@ -9,20 +9,20 @@ using JET
 # Test implicit Euler using ImplicitDiscreteProblem
 @testset "Implicit Euler" begin
     function lotkavolterra(u, p, t)
-        [1.5*u[1] - u[1]*u[2], -3.0*u[2] + u[1]*u[2]]
+        [1.5 * u[1] - u[1] * u[2], -3.0 * u[2] + u[1] * u[2]]
     end
 
     function f!(resid, u_next, u, p, t)
         lv = lotkavolterra(u_next, p, t)
-        resid[1] = u_next[1] - u[1] - 0.01*lv[1]
-        resid[2] = u_next[2] - u[2] - 0.01*lv[2]
+        resid[1] = u_next[1] - u[1] - 0.01 * lv[1]
+        resid[2] = u_next[2] - u[2] - 0.01 * lv[2]
         nothing
     end
     u0 = [1.0, 1.0]
     tspan = (0.0, 0.5)
 
     idprob = ImplicitDiscreteProblem(f!, u0, tspan, []; dt = 0.01)
-    idsol = solve(idprob, IDSolve())
+    idsol = solve(idprob, IDSolve(); adaptive = false)
 
     oprob = ODEProblem(lotkavolterra, u0, tspan)
     osol = solve(oprob, ImplicitEuler())
@@ -37,15 +37,15 @@ using JET
 
     function g!(resid, u_next, u, p, t)
         f = ff(u_next, p, t)
-        resid[1] = u_next[1] - u[1] - 0.01*f[1]
-        resid[2] = u_next[2] - u[2] - 0.01*f[2]
+        resid[1] = u_next[1] - u[1] - 0.01 * f[1]
+        resid[2] = u_next[2] - u[2] - 0.01 * f[2]
         nothing
     end
     u0 = [10.0, 0.0]
     tspan = (0, 0.2)
 
     idprob = ImplicitDiscreteProblem(g!, u0, tspan, []; dt = 0.01)
-    idsol = solve(idprob, IDSolve())
+    idsol = solve(idprob, IDSolve(); adaptive = false)
 
     oprob = ODEProblem(ff, u0, tspan)
     osol = solve(oprob, ImplicitEuler())
@@ -55,7 +55,7 @@ end
 
 @testset "Solver initializes" begin
     function periodic!(resid, u_next, u, p, t)
-        resid[1] = u_next[1] - u[1] - sin(t*π/4)
+        resid[1] = u_next[1] - u[1] - sin(t * π / 4)
         resid[2] = 16 - u_next[2]^2 - u_next[1]^2
     end
 
@@ -71,15 +71,34 @@ end
     end
 end
 
+@testset "Hard problem" begin
+    function hard!(resid, u, u_prev, p, t)
+        resid[1] = tanh((u[1] - 10t)^2) / 2
+    end
+
+    u0 = [0.0]
+    idprob = ImplicitDiscreteProblem(hard!, u0, (0.0, 1.0), [])
+    integrator = init(idprob, IDSolve())
+    idsol = solve!(integrator)
+    @test idsol.retcode == ReturnCode.Success
+end
+
 @testset "Handle nothing in u0" begin
-    function empty(u_next, u, p, t)
+    function emptyiip(residual, u_next, u, p, t) # TODO OOP variant does not work yet
+        nothing
+    end
+    function emptyoop(u_next, u, p, t) # TODO OOP variant does not work yet
         nothing
     end
 
     tsteps = 5
     u0 = nothing
-    idprob = ImplicitDiscreteProblem(empty, u0, (0, tsteps), [])
+    idprob = ImplicitDiscreteProblem(emptyiip, u0, (0, tsteps), [])
     @test_nowarn integ = init(idprob, IDSolve())
+
+    idprob2 = ImplicitDiscreteProblem(emptyoop, u0, (0, tsteps), [])
+    @test_throws AssertionError("Empty u not supported with out of place functions yet.") integ=init(
+        idprob2, IDSolve())
 end
 
 @testset "Create NonlinearLeastSquaresProblem" begin
@@ -92,7 +111,7 @@ end
     idprob = ImplicitDiscreteProblem(
         ImplicitDiscreteFunction(over, resid_prototype = zeros(3)), u0, (0, tsteps), [])
     integ = init(idprob, IDSolve())
-    @test integ.cache.prob isa NonlinearLeastSquaresProblem
+    @test integ.cache.nlcache.prob isa NonlinearLeastSquaresProblem
 
     function under(u_next, u, p, t)
         [u_next[1] - u_next[2] - 1]
@@ -100,7 +119,7 @@ end
     idprob = ImplicitDiscreteProblem(
         ImplicitDiscreteFunction(under; resid_prototype = zeros(1)), u0, (0, tsteps), [])
     integ = init(idprob, IDSolve())
-    @test integ.cache.prob isa NonlinearLeastSquaresProblem
+    @test integ.cache.nlcache.prob isa NonlinearLeastSquaresProblem
 
     function full(u_next, u, p, t)
         [u_next[1]^2 - 3, u_next[2] - u[1]]
@@ -108,7 +127,7 @@ end
     idprob = ImplicitDiscreteProblem(
         ImplicitDiscreteFunction(full; resid_prototype = zeros(2)), u0, (0, tsteps), [])
     integ = init(idprob, IDSolve())
-    @test integ.cache.prob isa NonlinearProblem
+    @test integ.cache.nlcache.prob isa NonlinearProblem
 end
 
 @testset "InitialFailure thrown" begin