diff --git a/core/src/Ribasim.jl b/core/src/Ribasim.jl
index f7e67fa32..83a0ab56b 100644
--- a/core/src/Ribasim.jl
+++ b/core/src/Ribasim.jl
@@ -22,8 +22,6 @@ using OrdinaryDiffEqCore:
     AbstractNLSolver,
     calculate_residuals!
 using DiffEqBase: DiffEqBase
-using OrdinaryDiffEqNonlinearSolve: OrdinaryDiffEqNonlinearSolve, relax!, _compute_rhs!
-using DiffEqBase: remake
 
 # Interface for defining and solving the ODE problem of the physical layer.
 using SciMLBase:
diff --git a/core/src/config.jl b/core/src/config.jl
index 0c3b7d1d5..a49ad3ae7 100644
--- a/core/src/config.jl
+++ b/core/src/config.jl
@@ -16,7 +16,7 @@ using LineSearch: BackTracking
 using Logging: LogLevel, Debug, Info, Warn, Error
 using ..Ribasim: Ribasim, isnode, nodetype
 using OrdinaryDiffEqCore: OrdinaryDiffEqAlgorithm, OrdinaryDiffEqNewtonAdaptiveAlgorithm
-using OrdinaryDiffEqNonlinearSolve: NLNewton, NonlinearSolveAlg, NewtonRaphson
+using OrdinaryDiffEqNonlinearSolve: NonlinearSolveAlg, NewtonRaphson
 using OrdinaryDiffEqLowOrderRK: Euler, RK4
 using OrdinaryDiffEqTsit5: Tsit5
 using OrdinaryDiffEqSDIRK: ImplicitEuler, KenCarp4, TRBDF2
@@ -289,9 +289,6 @@ function algorithm(solver::Solver; u0 = [])::OrdinaryDiffEqAlgorithm
     kwargs = Dict{Symbol, Any}()
 
     if algotype <: OrdinaryDiffEqNewtonAdaptiveAlgorithm
-        # kwargs[:nlsolve] = NLNewton(;
-        #     relax = Ribasim.MonitoredBackTracking(; z_tmp = copy(u0), dz_tmp = copy(u0)),
-        # )
         kwargs[:nlsolve] = NonlinearSolveAlg(
             NewtonRaphson(; linesearch = BackTracking(), autodiff = AutoFiniteDiff()),
         )
diff --git a/core/src/util.jl b/core/src/util.jl
index d12498c9b..27051bd85 100644
--- a/core/src/util.jl
+++ b/core/src/util.jl
@@ -859,71 +859,6 @@ relaxed_root(x::GradientTracer, threshold::Real) = x
 get_level_from_storage(basin::Basin, state_idx::Int, storage::GradientTracer) = storage
 stop_declining_negative_storage!(du, u::ComponentVector{<:GradientTracer}) = nothing
 
-@kwdef struct MonitoredBackTracking{B, V}
-    linesearch::B = BackTracking()
-    dz_tmp::V = []
-    z_tmp::V = []
-end
-
-"""
-Compute the residual of the non-linear solver, i.e. a measure of the
-error in the solution to the implicit equation defined by the solver algorithm
-"""
-function residual(z, integrator, nlsolver, f)
-    (; uprev, t, p, dt, opts, isdae) = integrator
-    (; tmp, ztmp, γ, α, cache, method) = nlsolver
-    (; ustep, atmp, tstep, k, invγdt, tstep, k, invγdt) = cache
-    if isdae
-        _uprev = get_dae_uprev(integrator, uprev)
-        b, ustep2 =
-            _compute_rhs!(tmp, ztmp, ustep, α, tstep, k, invγdt, p, _uprev, f::TF, z)
-    else
-        b, ustep2 =
-            _compute_rhs!(tmp, ztmp, ustep, γ, α, tstep, k, invγdt, method, p, dt, f, z)
-    end
-    calculate_residuals!(
-        atmp,
-        b,
-        uprev,
-        ustep2,
-        opts.abstol,
-        opts.reltol,
-        opts.internalnorm,
-        t,
-    )
-    ndz = opts.internalnorm(atmp, t)
-    return ndz
-end
-
-"""
-MonitoredBackTracing is a thin wrapper of BackTracking, making sure that
-the BackTracking relaxation is rejected if it results in a residual increase
-"""
-function OrdinaryDiffEqNonlinearSolve.relax!(
-    dz,
-    nlsolver::AbstractNLSolver,
-    integrator::DEIntegrator,
-    f,
-    linesearch::MonitoredBackTracking,
-)
-    (; linesearch, dz_tmp, z_tmp) = linesearch
-
-    # Store step before relaxation
-    @. dz_tmp = dz
-
-    # Apply relaxation and measure the residual change
-    @. z_tmp = nlsolver.z + dz
-    resid_before = residual(z_tmp, integrator, nlsolver, f)
-    relax!(dz, nlsolver, integrator, f, linesearch)
-    @. z_tmp = nlsolver.z + dz
-    resid_after = residual(z_tmp, integrator, nlsolver, f)
-
-    # If the residual increased due to the relaxation, reject it
-    if resid_after > resid_before
-        @. dz = dz_tmp
-    end
-end
-
 function build_state_vector(p::Parameters)
     # It is assumed that the horizontal flow states come first in
     # p.state_inflow_edge and p.state_outflow_edge