diff --git a/.JuliaFormatter.toml b/.JuliaFormatter.toml
index 453925c..4d06911 100644
--- a/.JuliaFormatter.toml
+++ b/.JuliaFormatter.toml
@@ -1 +1,4 @@
-style = "sciml"
\ No newline at end of file
+style = "sciml"
+format_markdown = true
+annotate_untyped_fields_with_any = false
+format_docstrings = true
\ No newline at end of file
diff --git a/Project.toml b/Project.toml
index aa45b60..3c6f51a 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,37 +1,30 @@
 name = "SimpleNonlinearSolve"
 uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 authors = ["SciML"]
-version = "0.1.25"
+version = "0.2.0"
 
 [deps]
+ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
+ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MaybeInplace = "bb5d69b7-63fc-4a16-80bd-7e42200c7bdb"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 
-[weakdeps]
-NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
-
-[extensions]
-SimpleNonlinearSolveNNlibExt = "NNlib"
-
 [compat]
 ArrayInterface = "7"
 DiffEqBase = "6.126"
 FiniteDiff = "2"
 ForwardDiff = "0.10.3"
 LinearAlgebra = "1.9"
-NNlib = "0.8, 0.9"
 PrecompileTools = "1"
 Reexport = "1"
 SciMLBase = "2.7"
 StaticArraysCore = "1.4"
 julia = "1.9"
-
-[extras]
-NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
diff --git a/README.md b/README.md
index 53b32e3..6bba38f 100644
--- a/README.md
+++ b/README.md
@@ -6,12 +6,12 @@
 [![codecov](https://codecov.io/gh/SciML/SimpleNonlinearSolve.jl/branch/main/graph/badge.svg)](https://codecov.io/gh/SciML/SimpleNonlinearSolve.jl)
 [![Build Status](https://github.com/SciML/SimpleNonlinearSolve.jl/workflows/CI/badge.svg)](https://github.com/SciML/SimpleNonlinearSolve.jl/actions?query=workflow%3ACI)
 
-[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor's%20Guide-blueviolet)](https://github.com/SciML/ColPrac)
+[![ColPrac: Contributor's Guide on Collaborative Practices for Community Packages](https://img.shields.io/badge/ColPrac-Contributor%27s%20Guide-blueviolet)](https://github.com/SciML/ColPrac)
 [![SciML Code Style](https://img.shields.io/static/v1?label=code%20style&message=SciML&color=9558b2&labelColor=389826)](https://github.com/SciML/SciMLStyle)
 
 Fast implementations of root finding algorithms in Julia that satisfy the SciML common interface.
 SimpleNonlinearSolve.jl focuses on low-dependency implementations of very fast methods for
-very small and simple problems. For the full set of solvers, see 
+very small and simple problems. For the full set of solvers, see
 [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl), of which
 SimpleNonlinearSolve.jl is just one solver set.
 
@@ -25,7 +25,7 @@ the documentation which contains the unreleased features.
 ```julia
 using SimpleNonlinearSolve, StaticArrays
 
-f(u,p) = u .* u .- 2
+f(u, p) = u .* u .- 2
 u0 = @SVector[1.0, 1.0]
 probN = NonlinearProblem{false}(f, u0)
 solver = solve(probN, SimpleNewtonRaphson(), abstol = 1e-9)
@@ -39,3 +39,16 @@ sol = solve(probB, ITP())
 ```
 
 For more details on the bracketing methods, refer to the [Tutorials](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/nonlinear/#Using-Bracketing-Methods) and detailed [APIs](https://docs.sciml.ai/NonlinearSolve/stable/api/simplenonlinearsolve/#Solver-API)
+
+## Breaking Changes in v2
+
+  - Batched solvers have been removed in favor of `BatchedArrays.jl`. Stay tuned for detailed
+    tutorials on how to use `BatchedArrays.jl` with `NonlinearSolve` & `SimpleNonlinearSolve`
+    solvers.
+  - The old style of specifying autodiff with `chunksize`, `standardtag`, etc. has been
+    deprecated in favor of directly specifying the autodiff type, like `AutoForwardDiff`.
+  - `Broyden` and `Klement` have been renamed to `SimpleBroyden` and `SimpleKlement` to
+    avoid conflicts with `NonlinearSolve.jl`'s `GeneralBroyden` and `GeneralKlement`, which
+    will be renamed to `Broyden` and `Klement` in the future.
+  - `LBroyden` has been renamed to `SimpleLimitedMemoryBroyden` to make it consistent with
+    `NonlinearSolve.jl`'s `LimitedMemoryBroyden`.
diff --git a/ext/SimpleNonlinearSolveNNlibExt.jl b/ext/SimpleNonlinearSolveNNlibExt.jl
deleted file mode 100644
index 1132b64..0000000
--- a/ext/SimpleNonlinearSolveNNlibExt.jl
+++ /dev/null
@@ -1,81 +0,0 @@
-module SimpleNonlinearSolveNNlibExt
-
-using ArrayInterface, DiffEqBase, LinearAlgebra, NNlib, SimpleNonlinearSolve, SciMLBase
-import SimpleNonlinearSolve: _construct_batched_problem_structure,
-    _get_storage, _init_𝓙, _result_from_storage, _get_tolerance, @maybeinplace
-
-function __init__()
-    SimpleNonlinearSolve.NNlibExtLoaded[] = true
-    return
-end
-
-@views function SciMLBase.__solve(prob::NonlinearProblem,
-        alg::BatchedBroyden;
-        abstol = nothing,
-        reltol = nothing,
-        maxiters = 1000,
-        kwargs...)
-    iip = isinplace(prob)
-
-    u, f, reconstruct = _construct_batched_problem_structure(prob)
-    L, N = size(u)
-
-    tc = alg.termination_condition
-    mode = DiffEqBase.get_termination_mode(tc)
-
-    storage = _get_storage(mode, u)
-
-    xₙ, xₙ₋₁, δx, δf = ntuple(_ -> copy(u), 4)
-    T = eltype(u)
-
-    atol = _get_tolerance(abstol, tc.abstol, T)
-    rtol = _get_tolerance(reltol, tc.reltol, T)
-    termination_condition = tc(storage)
-
-    𝓙⁻¹ = _init_𝓙(xₙ)  # L × L × N
-    𝓙⁻¹f, xᵀ𝓙⁻¹δf, xᵀ𝓙⁻¹ = similar(𝓙⁻¹, L, N), similar(𝓙⁻¹, 1, N), similar(𝓙⁻¹, 1, L, N)
-
-    @maybeinplace iip fₙ₋₁=f(xₙ) u
-    iip && (fₙ = copy(fₙ₋₁))
-    for n in 1:maxiters
-        batched_mul!(reshape(𝓙⁻¹f, L, 1, N), 𝓙⁻¹, reshape(fₙ₋₁, L, 1, N))
-        xₙ .= xₙ₋₁ .- 𝓙⁻¹f
-
-        @maybeinplace iip fₙ=f(xₙ)
-        δx .= xₙ .- xₙ₋₁
-        δf .= fₙ .- fₙ₋₁
-
-        batched_mul!(reshape(𝓙⁻¹f, L, 1, N), 𝓙⁻¹, reshape(δf, L, 1, N))
-        δxᵀ = reshape(δx, 1, L, N)
-
-        batched_mul!(reshape(xᵀ𝓙⁻¹δf, 1, 1, N), δxᵀ, reshape(𝓙⁻¹f, L, 1, N))
-        batched_mul!(xᵀ𝓙⁻¹, δxᵀ, 𝓙⁻¹)
-        δx .= (δx .- 𝓙⁻¹f) ./ (xᵀ𝓙⁻¹δf .+ T(1e-5))
-        batched_mul!(𝓙⁻¹, reshape(δx, L, 1, N), xᵀ𝓙⁻¹, one(T), one(T))
-
-        if termination_condition(fₙ, xₙ, xₙ₋₁, atol, rtol)
-            retcode, xₙ, fₙ = _result_from_storage(storage, xₙ, fₙ, f, mode, iip)
-            return DiffEqBase.build_solution(prob,
-                alg,
-                reconstruct(xₙ),
-                reconstruct(fₙ);
-                retcode)
-        end
-
-        xₙ₋₁ .= xₙ
-        fₙ₋₁ .= fₙ
-    end
-
-    if mode ∈ DiffEqBase.SAFE_BEST_TERMINATION_MODES
-        xₙ = storage.u
-        @maybeinplace iip fₙ=f(xₙ)
-    end
-
-    return DiffEqBase.build_solution(prob,
-        alg,
-        reconstruct(xₙ),
-        reconstruct(fₙ);
-        retcode = ReturnCode.MaxIters)
-end
-
-end
diff --git a/src/SimpleNonlinearSolve.jl b/src/SimpleNonlinearSolve.jl
index 8c84c43..a723b0a 100644
--- a/src/SimpleNonlinearSolve.jl
+++ b/src/SimpleNonlinearSolve.jl
@@ -1,90 +1,101 @@
 module SimpleNonlinearSolve
 
-using Reexport
-using FiniteDiff, ForwardDiff
-using ForwardDiff: Dual
-using StaticArraysCore
-using LinearAlgebra
-import ArrayInterface
-using DiffEqBase
-
-@reexport using SciMLBase
+import PrecompileTools: @compile_workload, @setup_workload, @recompile_invalidations
+
+@recompile_invalidations begin
+    using ADTypes,
+        ArrayInterface, ConcreteStructs, DiffEqBase, Reexport, LinearAlgebra, SciMLBase
+
+    import DiffEqBase: AbstractNonlinearTerminationMode,
+        AbstractSafeNonlinearTerminationMode, AbstractSafeBestNonlinearTerminationMode,
+        NonlinearSafeTerminationReturnCode, get_termination_mode,
+        NONLINEARSOLVE_DEFAULT_NORM
+    using FiniteDiff, ForwardDiff
+    import ForwardDiff: Dual
+    import MaybeInplace: @bb, setindex_trait, CanSetindex, CannotSetindex
+    import SciMLBase: AbstractNonlinearAlgorithm, build_solution, isinplace
+    import StaticArraysCore: StaticArray, SVector, SMatrix, SArray, MArray, MMatrix, Size
+end
 
-const NNlibExtLoaded = Ref{Bool}(false)
+@reexport using ADTypes, SciMLBase
 
-abstract type AbstractSimpleNonlinearSolveAlgorithm <: SciMLBase.AbstractNonlinearAlgorithm end
+abstract type AbstractSimpleNonlinearSolveAlgorithm <: AbstractNonlinearAlgorithm end
 abstract type AbstractBracketingAlgorithm <: AbstractSimpleNonlinearSolveAlgorithm end
-abstract type AbstractNewtonAlgorithm{CS, AD, FDT} <: AbstractSimpleNonlinearSolveAlgorithm end
-abstract type AbstractImmutableNonlinearSolver <: AbstractSimpleNonlinearSolveAlgorithm end
-abstract type AbstractBatchedNonlinearSolveAlgorithm <:
-              AbstractSimpleNonlinearSolveAlgorithm end
+abstract type AbstractNewtonAlgorithm <: AbstractSimpleNonlinearSolveAlgorithm end
 
 include("utils.jl")
-include("bisection.jl")
-include("falsi.jl")
-include("raphson.jl")
-include("broyden.jl")
-include("lbroyden.jl")
-include("klement.jl")
-include("trustRegion.jl")
-include("ridder.jl")
-include("brent.jl")
-include("dfsane.jl")
-include("ad.jl")
-include("halley.jl")
-include("alefeld.jl")
-include("itp.jl")
 
-# Batched Solver Support
-include("batched/utils.jl")
-include("batched/raphson.jl")
-include("batched/dfsane.jl")
-include("batched/broyden.jl")
+## Nonlinear Solvers
+include("nlsolve/raphson.jl")
+include("nlsolve/broyden.jl")
+include("nlsolve/lbroyden.jl")
+include("nlsolve/klement.jl")
+include("nlsolve/trustRegion.jl")
+include("nlsolve/halley.jl")
+include("nlsolve/dfsane.jl")
+
+## Interval Nonlinear Solvers
+include("bracketing/bisection.jl")
+include("bracketing/falsi.jl")
+include("bracketing/ridder.jl")
+include("bracketing/brent.jl")
+include("bracketing/alefeld.jl")
+include("bracketing/itp.jl")
+
+# AD
+include("ad.jl")
 
 ## Default algorithm
 
 # Set the default bracketing method to ITP
-
 function SciMLBase.solve(prob::IntervalNonlinearProblem; kwargs...)
-    SciMLBase.solve(prob, ITP(); kwargs...)
+    return solve(prob, ITP(); kwargs...)
 end
 
 function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Nothing,
         args...; kwargs...)
-    SciMLBase.solve(prob, ITP(), args...; kwargs...)
+    return solve(prob, ITP(), args...; kwargs...)
 end
 
-import PrecompileTools
-
-PrecompileTools.@compile_workload begin
+@setup_workload begin
     for T in (Float32, Float64)
-        prob_no_brack = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2))
-        for alg in (SimpleNewtonRaphson, SimpleHalley, Broyden, Klement, SimpleTrustRegion,
-            SimpleDFSane)
-            solve(prob_no_brack, alg(), abstol = T(1e-2))
-        end
+        prob_no_brack_scalar = NonlinearProblem{false}((u, p) -> u .* u .- p, T(0.1), T(2))
+        prob_no_brack_iip = NonlinearProblem{true}((du, u, p) -> du .= u .* u .- p,
+            T.([1.0, 1.0, 1.0]), T(2))
+        prob_no_brack_oop = NonlinearProblem{false}((u, p) -> u .* u .- p,
+            T.([1.0, 1.0, 1.0]), T(2))
+
+        algs = [SimpleNewtonRaphson(), SimpleBroyden(), SimpleKlement(), SimpleDFSane(),
+            SimpleTrustRegion(), SimpleLimitedMemoryBroyden(; threshold = 2)]
+
+        algs_no_iip = [SimpleHalley()]
+
+        @compile_workload begin
+            for alg in algs
+                solve(prob_no_brack_scalar, alg, abstol = T(1e-2))
+                solve(prob_no_brack_iip, alg, abstol = T(1e-2))
+                solve(prob_no_brack_oop, alg, abstol = T(1e-2))
+            end
 
-        #=
-        for alg in (SimpleNewtonRaphson,)
-            for u0 in ([1., 1.], StaticArraysCore.SA[1.0, 1.0])
-                u0 = T.(.1)
-                probN = NonlinearProblem{false}((u,p) -> u .* u .- p, u0, T(2))
-                solve(probN, alg(), tol = T(1e-2))
+            for alg in algs_no_iip
+                solve(prob_no_brack_scalar, alg, abstol = T(1e-2))
+                solve(prob_no_brack_oop, alg, abstol = T(1e-2))
             end
         end
-        =#
 
         prob_brack = IntervalNonlinearProblem{false}((u, p) -> u * u - p,
-            T.((0.0, 2.0)),
-            T(2))
-        for alg in (Bisection, Falsi, Ridder, Brent, Alefeld, ITP)
-            solve(prob_brack, alg(), abstol = T(1e-2))
+            T.((0.0, 2.0)), T(2))
+        algs = [Bisection(), Falsi(), Ridder(), Brent(), Alefeld(), ITP()]
+        @compile_workload begin
+            for alg in algs
+                solve(prob_brack, alg, abstol = T(1e-2))
+            end
         end
     end
 end
 
-export Bisection, Brent, Broyden, LBroyden, SimpleDFSane, Falsi, SimpleHalley, Klement,
-    Ridder, SimpleNewtonRaphson, SimpleTrustRegion, Alefeld, ITP, SimpleGaussNewton
-export BatchedBroyden, BatchedSimpleNewtonRaphson, BatchedSimpleDFSane
+export SimpleBroyden, SimpleDFSane, SimpleGaussNewton, SimpleHalley, SimpleKlement,
+    SimpleLimitedMemoryBroyden, SimpleNewtonRaphson, SimpleTrustRegion
+export Alefeld, Bisection, Brent, Falsi, ITP, Ridder
 
 end # module
diff --git a/src/ad.jl b/src/ad.jl
index b0fd9f1..d4cbcf7 100644
--- a/src/ad.jl
+++ b/src/ad.jl
@@ -1,9 +1,8 @@
 function scalar_nlsolve_ad(prob, alg, args...; kwargs...)
     f = prob.f
     p = value(prob.p)
-
     if prob isa IntervalNonlinearProblem
-        tspan = value(prob.tspan)
+        tspan = value.(prob.tspan)
         newprob = IntervalNonlinearProblem(f, tspan, p; prob.kwargs...)
     else
         u0 = value(prob.u0)
@@ -13,66 +12,78 @@ function scalar_nlsolve_ad(prob, alg, args...; kwargs...)
     sol = solve(newprob, alg, args...; kwargs...)
 
     uu = sol.u
-    if p isa Number
-        f_p = ForwardDiff.derivative(Base.Fix1(f, uu), p)
-    else
-        f_p = ForwardDiff.gradient(Base.Fix1(f, uu), p)
-    end
+    f_p = scalar_nlsolve_∂f_∂p(f, uu, p)
+    f_x = scalar_nlsolve_∂f_∂u(f, uu, p)
+
+    z_arr = -inv(f_x) * f_p
 
-    f_x = ForwardDiff.derivative(Base.Fix2(f, p), uu)
     pp = prob.p
-    sumfun = let f_x′ = -f_x
-        ((fp, p),) -> (fp / f_x′) * ForwardDiff.partials(p)
+    sumfun = ((z, p),) -> map(zᵢ -> zᵢ * ForwardDiff.partials(p), z)
+    if uu isa Number
+        partials = sum(sumfun, zip(z_arr, pp))
+    elseif p isa Number
+        partials = sumfun((z_arr, pp))
+    else
+        partials = sum(sumfun, zip(eachcol(z_arr), pp))
     end
-    partials = sum(sumfun, zip(f_p, pp))
+
     return sol, partials
 end
 
-function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector},
-            iip,
-            <:Dual{T, V, P}},
-        alg::AbstractSimpleNonlinearSolveAlgorithm,
-        args...; kwargs...) where {iip, T, V, P}
+function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, SVector, <:AbstractArray},
+            false, <:Dual{T, V, P}}, alg::AbstractSimpleNonlinearSolveAlgorithm, args...;
+        kwargs...) where {T, V, P}
     sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
-    return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid;
-        retcode = sol.retcode)
+    dual_soln = scalar_nlsolve_dual_soln(sol.u, partials, prob.p)
+    return SciMLBase.build_solution(prob, alg, dual_soln, sol.resid; sol.retcode)
 end
-function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, StaticArraysCore.SVector},
-            iip,
-            <:AbstractArray{<:Dual{T, V, P}}},
-        alg::AbstractSimpleNonlinearSolveAlgorithm, args...;
-        kwargs...) where {iip, T, V, P}
+
+function SciMLBase.solve(prob::NonlinearProblem{<:Union{Number, SVector, <:AbstractArray},
+            false, <:AbstractArray{<:Dual{T, V, P}}},
+        alg::AbstractSimpleNonlinearSolveAlgorithm, args...; kwargs...) where {T, V, P}
     sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
-    return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials), sol.resid;
-        retcode = sol.retcode)
+    dual_soln = scalar_nlsolve_dual_soln(sol.u, partials, prob.p)
+    return SciMLBase.build_solution(prob, alg, dual_soln, sol.resid; sol.retcode)
+end
+
+function scalar_nlsolve_∂f_∂p(f, u, p)
+    ff = p isa Number ? ForwardDiff.derivative :
+         (u isa Number ? ForwardDiff.gradient : ForwardDiff.jacobian)
+    return ff(Base.Fix1(f, u), p)
+end
+
+function scalar_nlsolve_∂f_∂u(f, u, p)
+    ff = u isa Number ? ForwardDiff.derivative : ForwardDiff.jacobian
+    return ff(Base.Fix2(f, p), u)
+end
+
+function scalar_nlsolve_dual_soln(u::Number, partials,
+        ::Union{<:AbstractArray{<:Dual{T, V, P}}, Dual{T, V, P}}) where {T, V, P}
+    return Dual{T, V, P}(u, partials)
+end
+
+function scalar_nlsolve_dual_soln(u::AbstractArray, partials,
+        ::Union{<:AbstractArray{<:Dual{T, V, P}}, Dual{T, V, P}}) where {T, V, P}
+    return map(((uᵢ, pᵢ),) -> Dual{T, V, P}(uᵢ, pᵢ), zip(u, partials))
 end
 
 # avoid ambiguities
 for Alg in [Bisection]
     @eval function SciMLBase.solve(prob::IntervalNonlinearProblem{uType, iip,
-                <:Dual{T, V, P}},
-            alg::$Alg, args...;
-            kwargs...) where {uType, iip, T, V, P}
+                <:Dual{T, V, P}}, alg::$Alg, args...; kwargs...) where {uType, iip, T, V, P}
         sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
-        return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials),
-            sol.resid; retcode = sol.retcode,
+        dual_soln = scalar_nlsolve_dual_soln(sol.u, partials, prob.p)
+        return SciMLBase.build_solution(prob, alg, dual_soln, sol.resid; sol.retcode,
             left = Dual{T, V, P}(sol.left, partials),
             right = Dual{T, V, P}(sol.right, partials))
-        #return BracketingSolution(Dual{T,V,P}(sol.left, partials), Dual{T,V,P}(sol.right, partials), sol.retcode, sol.resid)
     end
     @eval function SciMLBase.solve(prob::IntervalNonlinearProblem{uType, iip,
-                <:AbstractArray{
-                    <:Dual{T,
-                        V,
-                        P},
-                }},
-            alg::$Alg, args...;
+                <:AbstractArray{<:Dual{T, V, P}}}, alg::$Alg, args...;
             kwargs...) where {uType, iip, T, V, P}
         sol, partials = scalar_nlsolve_ad(prob, alg, args...; kwargs...)
-        return SciMLBase.build_solution(prob, alg, Dual{T, V, P}(sol.u, partials),
-            sol.resid; retcode = sol.retcode,
+        dual_soln = scalar_nlsolve_dual_soln(sol.u, partials, prob.p)
+        return SciMLBase.build_solution(prob, alg, dual_soln, sol.resid; sol.retcode,
             left = Dual{T, V, P}(sol.left, partials),
             right = Dual{T, V, P}(sol.right, partials))
-        #return BracketingSolution(Dual{T,V,P}(sol.left, partials), Dual{T,V,P}(sol.right, partials), sol.retcode, sol.resid)
     end
 end
diff --git a/src/batched/broyden.jl b/src/batched/broyden.jl
deleted file mode 100644
index ed3cd5d..0000000
--- a/src/batched/broyden.jl
+++ /dev/null
@@ -1,6 +0,0 @@
-struct BatchedBroyden{TC <: NLSolveTerminationCondition} <:
-       AbstractBatchedNonlinearSolveAlgorithm
-    termination_condition::TC
-end
-
-# Implementation of solve using Package Extensions
diff --git a/src/batched/dfsane.jl b/src/batched/dfsane.jl
deleted file mode 100644
index 01b3b19..0000000
--- a/src/batched/dfsane.jl
+++ /dev/null
@@ -1,141 +0,0 @@
-Base.@kwdef struct BatchedSimpleDFSane{T, F, TC <: NLSolveTerminationCondition} <:
-                   AbstractBatchedNonlinearSolveAlgorithm
-    σₘᵢₙ::T = 1.0f-10
-    σₘₐₓ::T = 1.0f+10
-    σ₁::T = 1.0f0
-    M::Int = 10
-    γ::T = 1.0f-4
-    τₘᵢₙ::T = 0.1f0
-    τₘₐₓ::T = 0.5f0
-    nₑₓₚ::Int = 2
-    ηₛ::F = (f₍ₙₒᵣₘ₎₁, n, xₙ, fₙ) -> f₍ₙₒᵣₘ₎₁ ./ n .^ 2
-    termination_condition::TC = NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
-        abstol = nothing,
-        reltol = nothing)
-    max_inner_iterations::Int = 1000
-end
-
-function SciMLBase.__solve(prob::NonlinearProblem,
-        alg::BatchedSimpleDFSane,
-        args...;
-        abstol = nothing,
-        reltol = nothing,
-        maxiters = 100,
-        kwargs...)
-    iip = isinplace(prob)
-
-    u, f, reconstruct = _construct_batched_problem_structure(prob)
-    L, N = size(u)
-    T = eltype(u)
-
-    tc = alg.termination_condition
-    mode = DiffEqBase.get_termination_mode(tc)
-
-    storage = _get_storage(mode, u)
-
-    atol = _get_tolerance(abstol, tc.abstol, T)
-    rtol = _get_tolerance(reltol, tc.reltol, T)
-    termination_condition = tc(storage)
-
-    σₘᵢₙ, σₘₐₓ, γ, τₘᵢₙ, τₘₐₓ = T(alg.σₘᵢₙ), T(alg.σₘₐₓ), T(alg.γ), T(alg.τₘᵢₙ), T(alg.τₘₐₓ)
-    α₁ = one(T)
-    α₊, α₋ = similar(u, 1, N), similar(u, 1, N)
-    σₙ = fill(T(alg.σ₁), 1, N)
-    𝒹 = similar(σₙ, L, N)
-    M = alg.M
-    nₑₓₚ = alg.nₑₓₚ
-
-    xₙ, xₙ₋₁, f₍ₙₒᵣₘ₎ₙ₋₁, f₍ₙₒᵣₘ₎ₙ = copy(u), copy(u), similar(u, 1, N), similar(u, 1, N)
-
-    function ff!(fₓ, fₙₒᵣₘ, x)
-        f(fₓ, x)
-        sum!(abs2, fₙₒᵣₘ, fₓ)
-        fₙₒᵣₘ .^= (nₑₓₚ / 2)
-        return fₓ
-    end
-
-    function ff!(fₙₒᵣₘ, x)
-        fₓ = f(x)
-        sum!(abs2, fₙₒᵣₘ, fₓ)
-        fₙₒᵣₘ .^= (nₑₓₚ / 2)
-        return fₓ
-    end
-
-    @maybeinplace iip fₙ₋₁=ff!(f₍ₙₒᵣₘ₎ₙ₋₁, xₙ) xₙ
-    iip && (fₙ = similar(fₙ₋₁))
-    ℋ = repeat(f₍ₙₒᵣₘ₎ₙ₋₁, M, 1)
-    f̄ = similar(ℋ, 1, N)
-    ηₛ = (n, xₙ, fₙ) -> alg.ηₛ(f₍ₙₒᵣₘ₎ₙ₋₁, n, xₙ, fₙ)
-
-    for n in 1:maxiters
-        # Spectral parameter range check
-        @. σₙ = sign(σₙ) * clamp(abs(σₙ), σₘᵢₙ, σₘₐₓ)
-
-        # Line search direction
-        @. 𝒹 = -σₙ * fₙ₋₁
-
-        η = ηₛ(n, xₙ₋₁, fₙ₋₁)
-        maximum!(f̄, ℋ)
-        fill!(α₊, α₁)
-        fill!(α₋, α₁)
-        @. xₙ = xₙ₋₁ + α₊ * 𝒹
-
-        @maybeinplace iip fₙ=ff!(f₍ₙₒᵣₘ₎ₙ, xₙ)
-
-        for _ in 1:(alg.max_inner_iterations)
-            𝒸 = @. f̄ + η - γ * α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁
-
-            (sum(f₍ₙₒᵣₘ₎ₙ .≤ 𝒸) ≥ N ÷ 2) && break
-
-            @. α₊ = clamp(α₊^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₊ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁),
-                τₘᵢₙ * α₊,
-                τₘₐₓ * α₊)
-            @. xₙ = xₙ₋₁ - α₋ * 𝒹
-            @maybeinplace iip fₙ=ff!(f₍ₙₒᵣₘ₎ₙ, xₙ)
-
-            (sum(f₍ₙₒᵣₘ₎ₙ .≤ 𝒸) ≥ N ÷ 2) && break
-
-            @. α₋ = clamp(α₋^2 * f₍ₙₒᵣₘ₎ₙ₋₁ / (f₍ₙₒᵣₘ₎ₙ + (T(2) * α₋ - T(1)) * f₍ₙₒᵣₘ₎ₙ₋₁),
-                τₘᵢₙ * α₋,
-                τₘₐₓ * α₋)
-            @. xₙ = xₙ₋₁ + α₊ * 𝒹
-            @maybeinplace iip fₙ=ff!(f₍ₙₒᵣₘ₎ₙ, xₙ)
-        end
-
-        if termination_condition(fₙ, xₙ, xₙ₋₁, atol, rtol)
-            retcode, xₙ, fₙ = _result_from_storage(storage, xₙ, fₙ, f, mode, iip)
-            return DiffEqBase.build_solution(prob,
-                alg,
-                reconstruct(xₙ),
-                reconstruct(fₙ);
-                retcode)
-        end
-
-        # Update spectral parameter
-        @. xₙ₋₁ = xₙ - xₙ₋₁
-        @. fₙ₋₁ = fₙ - fₙ₋₁
-
-        sum!(abs2, α₊, xₙ₋₁)
-        sum!(α₋, xₙ₋₁ .* fₙ₋₁)
-        σₙ .= α₊ ./ (α₋ .+ T(1e-5))
-
-        # Take step
-        @. xₙ₋₁ = xₙ
-        @. fₙ₋₁ = fₙ
-        @. f₍ₙₒᵣₘ₎ₙ₋₁ = f₍ₙₒᵣₘ₎ₙ
-
-        # Update history
-        ℋ[n % M + 1, :] .= view(f₍ₙₒᵣₘ₎ₙ, 1, :)
-    end
-
-    if mode ∈ DiffEqBase.SAFE_BEST_TERMINATION_MODES
-        xₙ = storage.u
-        @maybeinplace iip fₙ=f(xₙ)
-    end
-
-    return DiffEqBase.build_solution(prob,
-        alg,
-        reconstruct(xₙ),
-        reconstruct(fₙ);
-        retcode = ReturnCode.MaxIters)
-end
diff --git a/src/batched/raphson.jl b/src/batched/raphson.jl
deleted file mode 100644
index 7bc7b8c..0000000
--- a/src/batched/raphson.jl
+++ /dev/null
@@ -1,92 +0,0 @@
-struct BatchedSimpleNewtonRaphson{CS, AD, FDT, TC <: NLSolveTerminationCondition} <:
-       AbstractBatchedNonlinearSolveAlgorithm
-    termination_condition::TC
-end
-
-alg_autodiff(alg::BatchedSimpleNewtonRaphson{CS, AD, FDT}) where {CS, AD, FDT} = AD
-diff_type(alg::BatchedSimpleNewtonRaphson{CS, AD, FDT}) where {CS, AD, FDT} = FDT
-
-function BatchedSimpleNewtonRaphson(; chunk_size = Val{0}(),
-        autodiff = Val{true}(),
-        diff_type = Val{:forward},
-        termination_condition = NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
-            abstol = nothing,
-            reltol = nothing))
-    return BatchedSimpleNewtonRaphson{SciMLBase._unwrap_val(chunk_size),
-        SciMLBase._unwrap_val(autodiff),
-        SciMLBase._unwrap_val(diff_type), typeof(termination_condition)}(termination_condition)
-end
-
-function SciMLBase.__solve(prob::NonlinearProblem, alg::BatchedSimpleNewtonRaphson;
-        abstol = nothing, reltol = nothing, maxiters = 1000, kwargs...)
-    iip = SciMLBase.isinplace(prob)
-    iip &&
-        @assert alg_autodiff(alg) "Inplace BatchedSimpleNewtonRaphson currently only supports autodiff."
-    u, f, reconstruct = _construct_batched_problem_structure(prob)
-
-    tc = alg.termination_condition
-    mode = DiffEqBase.get_termination_mode(tc)
-
-    storage = _get_storage(mode, u)
-
-    xₙ, xₙ₋₁ = copy(u), copy(u)
-    T = eltype(u)
-
-    atol = _get_tolerance(abstol, tc.abstol, T)
-    rtol = _get_tolerance(reltol, tc.reltol, T)
-    termination_condition = tc(storage)
-
-    if iip
-        𝓙 = similar(xₙ, length(xₙ), length(xₙ))
-        fₙ = similar(xₙ)
-        jac_cfg = ForwardDiff.JacobianConfig(f, fₙ, xₙ)
-    end
-
-    for i in 1:maxiters
-        if iip
-            value_derivative!(𝓙, fₙ, f, xₙ, jac_cfg)
-        else
-            if alg_autodiff(alg)
-                fₙ, 𝓙 = value_derivative(f, xₙ)
-            else
-                fₙ = f(xₙ)
-                𝓙 = FiniteDiff.finite_difference_jacobian(f,
-                    xₙ,
-                    diff_type(alg),
-                    eltype(xₙ),
-                    fₙ)
-            end
-        end
-
-        iszero(fₙ) && return DiffEqBase.build_solution(prob,
-            alg,
-            reconstruct(xₙ),
-            reconstruct(fₙ);
-            retcode = ReturnCode.Success)
-
-        δx = reshape(𝓙 \ vec(fₙ), size(xₙ))
-        xₙ .-= δx
-
-        if termination_condition(fₙ, xₙ, xₙ₋₁, atol, rtol)
-            retcode, xₙ, fₙ = _result_from_storage(storage, xₙ, fₙ, f, mode, iip)
-            return DiffEqBase.build_solution(prob,
-                alg,
-                reconstruct(xₙ),
-                reconstruct(fₙ);
-                retcode)
-        end
-
-        xₙ₋₁ .= xₙ
-    end
-
-    if mode ∈ DiffEqBase.SAFE_BEST_TERMINATION_MODES
-        xₙ = storage.u
-        @maybeinplace iip fₙ=f(xₙ)
-    end
-
-    return DiffEqBase.build_solution(prob,
-        alg,
-        reconstruct(xₙ),
-        reconstruct(fₙ);
-        retcode = ReturnCode.MaxIters)
-end
diff --git a/src/batched/utils.jl b/src/batched/utils.jl
deleted file mode 100644
index b8e66fe..0000000
--- a/src/batched/utils.jl
+++ /dev/null
@@ -1,79 +0,0 @@
-macro maybeinplace(iip::Symbol, expr::Expr, u0::Union{Symbol, Nothing} = nothing)
-    @assert expr.head == :(=)
-    x1, x2 = expr.args
-    @assert x2.head == :call
-    f, x... = x2.args
-    define_expr = u0 === nothing ? :() : :($(x1) = similar($(u0)))
-    return quote
-        if $(esc(iip))
-            $(esc(define_expr))
-            $(esc(f))($(esc(x1)), $(esc.(x)...))
-        else
-            $(esc(expr))
-        end
-    end
-end
-
-function _get_tolerance(η, tc_η, ::Type{T}) where {T}
-    fallback_η = real(oneunit(T)) * (eps(real(one(T))))^(4 // 5)
-    return ifelse(η !== nothing, η, ifelse(tc_η !== nothing, tc_η, fallback_η))
-end
-
-function _construct_batched_problem_structure(prob)
-    return _construct_batched_problem_structure(prob.u0,
-        prob.f,
-        prob.p,
-        Val(SciMLBase.isinplace(prob)))
-end
-
-function _construct_batched_problem_structure(u0::AbstractArray{T, N},
-        f,
-        p,
-        ::Val{iip}) where {T, N, iip}
-    # Reconstruct `u`
-    reconstruct = N == 2 ? identity : Base.Fix2(reshape, size(u0))
-    # Standardize `u`
-    standardize = N == 2 ? identity :
-                  (N == 1 ? Base.Fix2(reshape, (:, 1)) :
-                   Base.Fix2(reshape, (:, size(u0, ndims(u0)))))
-    # Updated Function
-    f_modified = if iip
-        function f_modified_iip(du, u)
-            f(reconstruct(du), reconstruct(u), p)
-            return standardize(du)
-        end
-    else
-        f_modified_oop(u) = standardize(f(reconstruct(u), p))
-    end
-    return standardize(u0), f_modified, reconstruct
-end
-
-@views function _init_𝓙(x::AbstractMatrix)
-    𝓙 = ArrayInterface.zeromatrix(x[:, 1])
-    if ismutable(x)
-        𝓙[diagind(𝓙)] .= one(eltype(x))
-    else
-        𝓙 .+= I
-    end
-    return repeat(𝓙, 1, 1, size(x, 2))
-end
-
-_result_from_storage(::Nothing, xₙ, fₙ, args...) = ReturnCode.Success, xₙ, fₙ
-function _result_from_storage(storage::NLSolveSafeTerminationResult, xₙ, fₙ, f, mode, iip)
-    if storage.return_code == DiffEqBase.NLSolveSafeTerminationReturnCode.Success
-        return ReturnCode.Success, xₙ, fₙ
-    else
-        if mode ∈ DiffEqBase.SAFE_BEST_TERMINATION_MODES
-            @maybeinplace iip fₙ=f(xₙ)
-            return ReturnCode.Terminated, storage.u, fₙ
-        else
-            return ReturnCode.Terminated, xₙ, fₙ
-        end
-    end
-end
-
-function _get_storage(mode, u)
-    return mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ?
-           NLSolveSafeTerminationResult(mode ∈ DiffEqBase.SAFE_BEST_TERMINATION_MODES ? u :
-                                        nothing) : nothing
-end
diff --git a/src/bisection.jl b/src/bisection.jl
deleted file mode 100644
index f7c98aa..0000000
--- a/src/bisection.jl
+++ /dev/null
@@ -1,92 +0,0 @@
-"""
-`Bisection(; exact_left = false, exact_right = false)`
-
-A common bisection method.
-
-### Keyword Arguments
-
-- `exact_left`: whether to enforce whether the left side of the interval must be exactly
-  zero for the returned result. Defaults to false.
-- `exact_right`: whether to enforce whether the right side of the interval must be exactly
-  zero for the returned result. Defaults to false.
-"""
-struct Bisection <: AbstractBracketingAlgorithm
-    exact_left::Bool
-    exact_right::Bool
-end
-
-function Bisection(; exact_left = false, exact_right = false)
-    Bisection(exact_left, exact_right)
-end
-
-function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args...;
-        maxiters = 1000, abstol = min(eps(prob.tspan[1]), eps(prob.tspan[2])),
-        kwargs...)
-    f = Base.Fix2(prob.f, prob.p)
-    left, right = prob.tspan
-    fl, fr = f(left), f(right)
-    if iszero(fl)
-        return SciMLBase.build_solution(prob, alg, left, fl;
-            retcode = ReturnCode.ExactSolutionLeft, left = left,
-            right = right)
-    end
-    if iszero(fr)
-        return SciMLBase.build_solution(prob, alg, right, fr;
-            retcode = ReturnCode.ExactSolutionRight, left = left,
-            right = right)
-    end
-
-    i = 1
-    if !iszero(fr)
-        while i < maxiters
-            mid = (left + right) / 2
-            (mid == left || mid == right) &&
-                return SciMLBase.build_solution(prob, alg, left, fl;
-                    retcode = ReturnCode.FloatingPointLimit,
-                    left = left, right = right)
-            fm = f(mid)
-            if abs((right - left) / 2) < abstol
-                return SciMLBase.build_solution(prob, alg, mid, fm;
-                    retcode = ReturnCode.Success,
-                    left = left, right = right)
-            end
-            if iszero(fm)
-                right = mid
-                break
-            end
-            if sign(fl) == sign(fm)
-                fl = fm
-                left = mid
-            else
-                fr = fm
-                right = mid
-            end
-            i += 1
-        end
-    end
-
-    while i < maxiters
-        mid = (left + right) / 2
-        (mid == left || mid == right) &&
-            return SciMLBase.build_solution(prob, alg, left, fl;
-                retcode = ReturnCode.FloatingPointLimit,
-                left = left, right = right)
-        fm = f(mid)
-        if abs((right - left) / 2) < abstol
-            return SciMLBase.build_solution(prob, alg, mid, fm;
-                retcode = ReturnCode.Success,
-                left = left, right = right)
-        end
-        if iszero(fm)
-            right = mid
-            fr = fm
-        else
-            left = mid
-            fl = fm
-        end
-        i += 1
-    end
-
-    return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters,
-        left = left, right = right)
-end
diff --git a/src/alefeld.jl b/src/bracketing/alefeld.jl
similarity index 68%
rename from src/alefeld.jl
rename to src/bracketing/alefeld.jl
index 0d4f561..39c984f 100644
--- a/src/alefeld.jl
+++ b/src/bracketing/alefeld.jl
@@ -1,31 +1,25 @@
 """
-`Alefeld()` 
+    Alefeld()
 
 An implementation of algorithm 4.2 from [Alefeld](https://dl.acm.org/doi/10.1145/210089.210111).
 
-The paper brought up two new algorithms. Here choose to implement algorithm 4.2 rather than 
+The paper brought up two new algorithms. Here choose to implement algorithm 4.2 rather than
 algorithm 4.1 because, in certain sense, the second algorithm(4.2) is an optimal procedure.
 """
 struct Alefeld <: AbstractBracketingAlgorithm end
 
-function SciMLBase.solve(prob::IntervalNonlinearProblem,
-        alg::Alefeld, args...; abstol = nothing,
-        reltol = nothing,
-        maxiters = 1000, kwargs...)
+function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Alefeld, args...;
+        maxiters = 1000, abstol = nothing, kwargs...)
     f = Base.Fix2(prob.f, prob.p)
     a, b = prob.tspan
     c = a - (b - a) / (f(b) - f(a)) * f(a)
 
     fc = f(c)
     (a == c || b == c) &&
-        return SciMLBase.build_solution(prob, alg, c, fc;
-            retcode = ReturnCode.FloatingPointLimit,
-            left = a,
-            right = b)
+        return build_solution(prob, alg, c, fc; retcode = ReturnCode.FloatingPointLimit,
+            left = a, right = b)
     iszero(fc) &&
-        return SciMLBase.build_solution(prob, alg, c, fc;
-            retcode = ReturnCode.Success,
-            left = a,
+        return build_solution(prob, alg, c, fc; retcode = ReturnCode.Success, left = a,
             right = b)
     a, b, d = _bracket(f, a, b, c)
     e = zero(a)   # Set e as 0 before iteration to avoid a non-value f(e)
@@ -44,15 +38,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
         end
         ē, fc = d, f(c)
         (a == c || b == c) &&
-            return SciMLBase.build_solution(prob, alg, c, fc;
-                retcode = ReturnCode.FloatingPointLimit,
-                left = a,
-                right = b)
+            return build_solution(prob, alg, c, fc; retcode = ReturnCode.FloatingPointLimit,
+                left = a, right = b)
         iszero(fc) &&
-            return SciMLBase.build_solution(prob, alg, c, fc;
-                retcode = ReturnCode.Success,
-                left = a,
-                right = b)
+            return build_solution(prob, alg, c, fc; retcode = ReturnCode.Success,
+                left = a, right = b)
         ā, b̄, d̄ = _bracket(f, a, b, c)
 
         # The second bracketing block
@@ -67,15 +57,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
         end
         fc = f(c)
         (ā == c || b̄ == c) &&
-            return SciMLBase.build_solution(prob, alg, c, fc;
-                retcode = ReturnCode.FloatingPointLimit,
-                left = ā,
-                right = b̄)
+            return build_solution(prob, alg, c, fc; retcode = ReturnCode.FloatingPointLimit,
+                left = ā, right = b̄)
         iszero(fc) &&
-            return SciMLBase.build_solution(prob, alg, c, fc;
-                retcode = ReturnCode.Success,
-                left = ā,
-                right = b̄)
+            return build_solution(prob, alg, c, fc; retcode = ReturnCode.Success,
+                left = ā, right = b̄)
         ā, b̄, d̄ = _bracket(f, ā, b̄, c)
 
         # The third bracketing block
@@ -90,15 +76,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
         end
         fc = f(c)
         (ā == c || b̄ == c) &&
-            return SciMLBase.build_solution(prob, alg, c, fc;
-                retcode = ReturnCode.FloatingPointLimit,
-                left = ā,
-                right = b̄)
+            return build_solution(prob, alg, c, fc; retcode = ReturnCode.FloatingPointLimit,
+                left = ā, right = b̄)
         iszero(fc) &&
-            return SciMLBase.build_solution(prob, alg, c, fc;
-                retcode = ReturnCode.Success,
-                left = ā,
-                right = b̄)
+            return build_solution(prob, alg, c, fc; retcode = ReturnCode.Success,
+                left = ā, right = b̄)
         ā, b̄, d = _bracket(f, ā, b̄, c)
 
         # The last bracketing block
@@ -109,15 +91,11 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
             c = 0.5 * (ā + b̄)
             fc = f(c)
             (ā == c || b̄ == c) &&
-                return SciMLBase.build_solution(prob, alg, c, fc;
-                    retcode = ReturnCode.FloatingPointLimit,
-                    left = ā,
-                    right = b̄)
+                return build_solution(prob, alg, c, fc;
+                    retcode = ReturnCode.FloatingPointLimit, left = ā, right = b̄)
             iszero(fc) &&
-                return SciMLBase.build_solution(prob, alg, c, fc;
-                    retcode = ReturnCode.Success,
-                    left = ā,
-                    right = b̄)
+                return build_solution(prob, alg, c, fc; retcode = ReturnCode.Success,
+                    left = ā, right = b̄)
             a, b, d = _bracket(f, ā, b̄, c)
         end
     end
@@ -131,7 +109,7 @@ function SciMLBase.solve(prob::IntervalNonlinearProblem,
     fc = f(c)
 
     # Reuturn solution when run out of max interation
-    return SciMLBase.build_solution(prob, alg, c, fc; retcode = ReturnCode.MaxIters,
+    return build_solution(prob, alg, c, fc; retcode = ReturnCode.MaxIters,
         left = a, right = b)
 end
 
diff --git a/src/bracketing/bisection.jl b/src/bracketing/bisection.jl
new file mode 100644
index 0000000..66418b3
--- /dev/null
+++ b/src/bracketing/bisection.jl
@@ -0,0 +1,107 @@
+"""
+    Bisection(; exact_left = false, exact_right = false)
+
+A common bisection method.
+
+### Keyword Arguments
+
+  - `exact_left`: whether to enforce whether the left side of the interval must be exactly
+    zero for the returned result. Defaults to false.
+  - `exact_right`: whether to enforce whether the right side of the interval must be exactly
+    zero for the returned result. Defaults to false.
+
+!!! warning
+
+    Currently, the keyword arguments are not implemented.
+"""
+@kwdef struct Bisection <: AbstractBracketingAlgorithm
+    exact_left::Bool = false
+    exact_right::Bool = false
+end
+
+function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Bisection, args...;
+        maxiters = 1000, abstol = nothing, kwargs...)
+    @assert !isinplace(prob) "`Bisection` only supports OOP problems."
+    f = Base.Fix2(prob.f, prob.p)
+    left, right = prob.tspan
+    fl, fr = f(left), f(right)
+
+    abstol = _get_tolerance(abstol,
+        promote_type(eltype(first(prob.tspan)), eltype(last(prob.tspan))))
+
+    if iszero(fl)
+        return build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft,
+            left, right)
+    end
+
+    if iszero(fr)
+        return build_solution(prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight,
+            left, right)
+    end
+
+    i = 1
+    if !iszero(fr)
+        while i < maxiters
+            mid = (left + right) / 2
+            (mid == left || mid == right) &&
+                return build_solution(prob, alg, left, fl; left, right,
+                    retcode = ReturnCode.FloatingPointLimit)
+            fm = f(mid)
+            if abs((right - left) / 2) < abstol
+                return build_solution(prob, alg, mid, fm; retcode = ReturnCode.Success,
+                    left, right)
+            end
+            if iszero(fm)
+                right = mid
+                break
+            end
+            if sign(fl) == sign(fm)
+                fl = fm
+                left = mid
+            else
+                fr = fm
+                right = mid
+            end
+            i += 1
+        end
+    end
+
+    sol, i, left, right, fl, fr = __bisection(left, right, fl, fr, f; abstol,
+        maxiters = maxiters - i, prob, alg)
+
+    sol !== nothing && return sol
+
+    return build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right)
+end
+
+function __bisection(left, right, fl, fr, f::F; abstol, maxiters, prob, alg) where {F}
+    i = 1
+    sol = nothing
+    while i < maxiters
+        mid = (left + right) / 2
+        if (mid == left || mid == right)
+            sol = build_solution(prob, alg, left, fl; left, right,
+                retcode = ReturnCode.FloatingPointLimit)
+            break
+        end
+
+        fm = f(mid)
+        if abs((right - left) / 2) < abstol
+            sol = build_solution(prob, alg, mid, fm; left, right,
+                retcode = ReturnCode.Success)
+            break
+        end
+
+        if iszero(fm)
+            right = mid
+            fr = fm
+        else
+            left = mid
+            fl = fm
+        end
+
+        i += 1
+    end
+
+    return sol, i, left, right, fl, fr
+end
diff --git a/src/bracketing/brent.jl b/src/bracketing/brent.jl
new file mode 100644
index 0000000..75497f3
--- /dev/null
+++ b/src/bracketing/brent.jl
@@ -0,0 +1,113 @@
+"""
+    Brent()
+
+left non-allocating Brent method.
+"""
+struct Brent <: AbstractBracketingAlgorithm end
+
+function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...;
+        maxiters = 1000, abstol = nothing, kwargs...)
+    @assert !isinplace(prob) "`Brent` only supports OOP problems."
+    f = Base.Fix2(prob.f, prob.p)
+    left, right = prob.tspan
+    fl, fr = f(left), f(right)
+    ϵ = eps(convert(typeof(fl), 1))
+
+    abstol = _get_tolerance(abstol,
+        promote_type(eltype(first(prob.tspan)), eltype(last(prob.tspan))))
+
+    if iszero(fl)
+        return build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft,
+            left, right)
+    end
+
+    if iszero(fr)
+        return build_solution(prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight,
+            left, right)
+    end
+
+    if abs(fl) < abs(fr)
+        c = right
+        right = left
+        left = c
+        tmp = fl
+        fl = fr
+        fr = tmp
+    end
+
+    c = left
+    d = c
+    i = 1
+    cond = true
+    if !iszero(fr)
+        while i < maxiters
+            fc = f(c)
+            if fl != fc && fr != fc
+                # Inverse quadratic interpolation
+                s = left * fr * fc / ((fl - fr) * (fl - fc)) +
+                    right * fl * fc / ((fr - fl) * (fr - fc)) +
+                    c * fl * fr / ((fc - fl) * (fc - fr))
+            else
+                # Secant method
+                s = right - fr * (right - left) / (fr - fl)
+            end
+            if (s < min((3 * left + right) / 4, right) ||
+                s > max((3 * left + right) / 4, right)) ||
+               (cond && abs(s - right) ≥ abs(right - c) / 2) ||
+               (!cond && abs(s - right) ≥ abs(c - d) / 2) ||
+               (cond && abs(right - c) ≤ ϵ) ||
+               (!cond && abs(c - d) ≤ ϵ)
+                # Bisection method
+                s = (left + right) / 2
+                (s == left || s == right) &&
+                    return SciMLBase.build_solution(prob, alg, left, fl;
+                        retcode = ReturnCode.FloatingPointLimit,
+                        left = left, right = right)
+                cond = true
+            else
+                cond = false
+            end
+            fs = f(s)
+            if abs((right - left) / 2) < abstol
+                return SciMLBase.build_solution(prob, alg, s, fs;
+                    retcode = ReturnCode.Success,
+                    left = left, right = right)
+            end
+            if iszero(fs)
+                if right < left
+                    left = right
+                    fl = fr
+                end
+                right = s
+                fr = fs
+                break
+            end
+            if fl * fs < 0
+                d = c
+                c = right
+                right = s
+                fr = fs
+            else
+                left = s
+                fl = fs
+            end
+            if abs(fl) < abs(fr)
+                d = c
+                c = right
+                right = left
+                left = c
+                fc = fr
+                fr = fl
+                fl = fc
+            end
+            i += 1
+        end
+    end
+
+    sol, i, left, right, fl, fr = __bisection(left, right, fl, fr, f; abstol,
+        maxiters = maxiters - i, prob, alg)
+
+    sol !== nothing && return sol
+
+    return build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right)
+end
diff --git a/src/bracketing/falsi.jl b/src/bracketing/falsi.jl
new file mode 100644
index 0000000..9db7d6c
--- /dev/null
+++ b/src/bracketing/falsi.jl
@@ -0,0 +1,70 @@
+"""
+    Falsi()
+
+A non-allocating regula falsi method.
+"""
+struct Falsi <: AbstractBracketingAlgorithm end
+
+function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Falsi, args...;
+        maxiters = 1000, abstol = nothing, kwargs...)
+    @assert !isinplace(prob) "`Falsi` only supports OOP problems."
+    f = Base.Fix2(prob.f, prob.p)
+    left, right = prob.tspan
+    fl, fr = f(left), f(right)
+
+    abstol = _get_tolerance(abstol,
+        promote_type(eltype(first(prob.tspan)), eltype(last(prob.tspan))))
+
+    if iszero(fl)
+        return build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft,
+            left, right)
+    end
+
+    if iszero(fr)
+        return build_solution(prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight,
+            left, right)
+    end
+
+    # Regula Falsi Steps
+    i = 0
+    if !iszero(fr)
+        while i < maxiters
+            if __nextfloat_tdir(left, prob.tspan...) == right
+                return build_solution(prob, alg, left, fl; left, right,
+                    retcode = ReturnCode.FloatingPointLimit)
+            end
+
+            mid = (fr * left - fl * right) / (fr - fl)
+            for _ in 1:10
+                mid = __max_tdir(left, __prevfloat_tdir(mid, prob.tspan...), prob.tspan...)
+            end
+
+            (mid == left || mid == right) && break
+
+            fm = f(mid)
+            if abs((right - left) / 2) < abstol
+                return build_solution(prob, alg, mid, fm; left, right,
+                    retcode = ReturnCode.Success)
+            end
+
+            if abs(fm) < abstol
+                right = mid
+                break
+            end
+
+            if sign(fl) == sign(fm)
+                fl, left = fm, mid
+            else
+                fr, right = fm, mid
+            end
+            i += 1
+        end
+    end
+
+    sol, i, left, right, fl, fr = __bisection(left, right, fl, fr, f; abstol,
+        maxiters = maxiters - i, prob, alg)
+    sol !== nothing && return sol
+
+    return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters,
+        left, right)
+end
diff --git a/src/bracketing/itp.jl b/src/bracketing/itp.jl
new file mode 100644
index 0000000..fd46de6
--- /dev/null
+++ b/src/bracketing/itp.jl
@@ -0,0 +1,143 @@
+"""
+    ITP(; k1::Real = 0.007, k2::Real = 1.5, n0::Int = 10)
+
+ITP (Interpolate Truncate & Project)
+
+Use the [ITP method](https://en.wikipedia.org/wiki/ITP_method) to find a root of a bracketed
+function, with a convergence rate between 1 and 1.62.
+
+This method was introduced in the paper "An Enhancement of the Bisection Method Average
+Performance Preserving Minmax Optimality" (https://doi.org/10.1145/3423597) by
+I. F. D. Oliveira and R. H. C. Takahashi.
+
+# Tuning Parameters
+
+The following keyword parameters are accepted.
+
+  - `n₀::Int = 1`, the 'slack'. Must not be negative. When n₀ = 0 the worst-case is
+    identical to that of bisection, but increacing n₀ provides greater oppotunity for
+    superlinearity.
+  - `κ₁::Float64 = 0.1`. Must not be negative. The recomended value is `0.2/(x₂ - x₁)`.
+    Lower values produce tighter asymptotic behaviour, while higher values improve the
+    steady-state behaviour when truncation is not helpful.
+  - `κ₂::Real = 2`. Must lie in [1, 1+ϕ ≈ 2.62). Higher values allow for a greater
+    convergence rate, but also make the method more succeptable to worst-case performance.
+    In practice, κ=1,2 seems to work well due to the computational simplicity, as κ₂ is used
+    as an exponent in the method.
+
+### Worst Case Performance
+
+n½ + `n₀` iterations, where n½ is the number of iterations using bisection
+(n½ = ⌈log2(Δx)/2`tol`⌉).
+
+### Asymptotic Performance
+
+If `f` is twice differentiable and the root is simple, then with `n₀` > 0 the convergence
+rate is √`κ₂`.
+"""
+struct ITP{T} <: AbstractBracketingAlgorithm
+    k1::T
+    k2::T
+    n0::Int
+    function ITP(; k1::Real = 0.007, k2::Real = 1.5, n0::Int = 10)
+        k1 < 0 && error("Hyper-parameter κ₁ should not be negative")
+        n0 < 0 && error("Hyper-parameter n₀ should not be negative")
+        if k2 < 1 || k2 > (1.5 + sqrt(5) / 2)
+            throw(ArgumentError("Hyper-parameter κ₂ should be between 1 and 1 + ϕ where \
+                                 ϕ ≈ 1.618... is the golden ratio"))
+        end
+        T = promote_type(eltype(k1), eltype(k2))
+        return new{T}(k1, k2, n0)
+    end
+end
+
+function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP, args...;
+        maxiters = 1000, abstol = nothing, kwargs...)
+    @assert !isinplace(prob) "`Bisection` only supports OOP problems."
+    f = Base.Fix2(prob.f, prob.p)
+    left, right = prob.tspan
+    fl, fr = f(left), f(right)
+
+    abstol = _get_tolerance(abstol,
+        promote_type(eltype(first(prob.tspan)), eltype(last(prob.tspan))))
+
+    if iszero(fl)
+        return build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft,
+            left, right)
+    end
+
+    if iszero(fr)
+        return build_solution(prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight,
+            left, right)
+    end
+    ϵ = abstol
+    #defining variables/cache
+    k1 = alg.k1
+    k2 = alg.k2
+    n0 = alg.n0
+    n_h = ceil(log2(abs(right - left) / (2 * ϵ)))
+    mid = (left + right) / 2
+    x_f = left + (right - left) * (fl / (fl - fr))
+    xt = left
+    xp = left
+    r = zero(left) #minmax radius
+    δ = zero(left) # truncation error
+    σ = 1.0
+    ϵ_s = ϵ * 2^(n_h + n0)
+    i = 0 #iteration
+    while i <= maxiters
+        span = abs(right - left)
+        r = ϵ_s - (span / 2)
+        δ = k1 * (span^k2)
+
+        ## Interpolation step ##
+        x_f = left + (right - left) * (fl / (fl - fr))
+
+        ## Truncation step ##
+        σ = sign(mid - x_f)
+        if δ <= abs(mid - x_f)
+            xt = x_f + (σ * δ)
+        else
+            xt = mid
+        end
+
+        ## Projection step ##
+        if abs(xt - mid) <= r
+            xp = xt
+        else
+            xp = mid - (σ * r)
+        end
+
+        if abs((left - right) / 2) < ϵ
+            return build_solution(prob, alg, mid, f(mid); retcode = ReturnCode.Success,
+                left, right)
+        end
+
+        ## Update ##
+        tmin, tmax = minmax(left, right)
+        xp >= tmax && (xp = prevfloat(tmax))
+        xp <= tmin && (xp = nextfloat(tmin))
+        yp = f(xp)
+        yps = yp * sign(fr)
+        if yps > 0
+            right = xp
+            fr = yp
+        elseif yps < 0
+            left = xp
+            fl = yp
+        else
+            return build_solution(prob, alg, xp, yps; retcode = ReturnCode.Success,
+                left = xp, right = xp)
+        end
+        i += 1
+        mid = (left + right) / 2
+        ϵ_s /= 2
+
+        if __nextfloat_tdir(left, prob.tspan...) == right
+            return build_solution(prob, alg, left, fl; left, right,
+                retcode = ReturnCode.FloatingPointLimit)
+        end
+    end
+
+    return build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters, left, right)
+end
diff --git a/src/bracketing/ridder.jl b/src/bracketing/ridder.jl
new file mode 100644
index 0000000..20e0db4
--- /dev/null
+++ b/src/bracketing/ridder.jl
@@ -0,0 +1,77 @@
+"""
+    Ridder()
+
+A non-allocating ridder method.
+"""
+struct Ridder <: AbstractBracketingAlgorithm end
+
+function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Ridder, args...;
+        maxiters = 1000, abstol = nothing, kwargs...)
+    @assert !isinplace(prob) "`Ridder` only supports OOP problems."
+    f = Base.Fix2(prob.f, prob.p)
+    left, right = prob.tspan
+    fl, fr = f(left), f(right)
+
+    abstol = _get_tolerance(abstol,
+        promote_type(eltype(first(prob.tspan)), eltype(last(prob.tspan))))
+
+    if iszero(fl)
+        return build_solution(prob, alg, left, fl; retcode = ReturnCode.ExactSolutionLeft,
+            left, right)
+    end
+
+    if iszero(fr)
+        return build_solution(prob, alg, right, fr; retcode = ReturnCode.ExactSolutionRight,
+            left, right)
+    end
+
+    xo = oftype(left, Inf)
+    i = 1
+    if !iszero(fr)
+        while i < maxiters
+            mid = (left + right) / 2
+            (mid == left || mid == right) &&
+                return build_solution(prob, alg, left, fl; left, right,
+                    retcode = ReturnCode.FloatingPointLimit)
+            fm = f(mid)
+            s = sqrt(fm^2 - fl * fr)
+            if iszero(s)
+                return build_solution(prob, alg, left, fl; left, right,
+                    retcode = ReturnCode.Failure)
+            end
+            x = mid + (mid - left) * sign(fl - fr) * fm / s
+            fx = f(x)
+            xo = x
+            if abs((right - left) / 2) < abstol
+                return build_solution(prob, alg, mid, fm; retcode = ReturnCode.Success,
+                    left, right)
+            end
+            if iszero(fx)
+                right = x
+                fr = fx
+                break
+            end
+            if sign(fx) != sign(fm)
+                left = mid
+                fl = fm
+                right = x
+                fr = fx
+            elseif sign(fx) != sign(fl)
+                right = x
+                fr = fx
+            else
+                @assert sign(fx) != sign(fr)
+                left = x
+                fl = fx
+            end
+            i += 1
+        end
+    end
+
+    sol, i, left, right, fl, fr = __bisection(left, right, fl, fr, f; abstol,
+        maxiters = maxiters - i, prob, alg)
+    sol !== nothing && return sol
+
+    return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters,
+        left, right)
+end
diff --git a/src/brent.jl b/src/brent.jl
deleted file mode 100644
index 7d7a6bc..0000000
--- a/src/brent.jl
+++ /dev/null
@@ -1,128 +0,0 @@
-"""
-`Brent()`
-
-A non-allocating Brent method
-
-"""
-struct Brent <: AbstractBracketingAlgorithm end
-
-function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Brent, args...;
-        maxiters = 1000, abstol = min(eps(prob.tspan[1]), eps(prob.tspan[2])),
-        kwargs...)
-    f = Base.Fix2(prob.f, prob.p)
-    a, b = prob.tspan
-    fa, fb = f(a), f(b)
-    ϵ = eps(convert(typeof(fa), 1.0))
-
-    if iszero(fa)
-        return SciMLBase.build_solution(prob, alg, a, fa;
-            retcode = ReturnCode.ExactSolutionLeft, left = a,
-            right = b)
-    elseif iszero(fb)
-        return SciMLBase.build_solution(prob, alg, b, fb;
-            retcode = ReturnCode.ExactSolutionRight, left = a,
-            right = b)
-    end
-    if abs(fa) < abs(fb)
-        c = b
-        b = a
-        a = c
-        tmp = fa
-        fa = fb
-        fb = tmp
-    end
-
-    c = a
-    d = c
-    i = 1
-    cond = true
-    if !iszero(fb)
-        while i < maxiters
-            fc = f(c)
-            if fa != fc && fb != fc
-                # Inverse quadratic interpolation
-                s = a * fb * fc / ((fa - fb) * (fa - fc)) +
-                    b * fa * fc / ((fb - fa) * (fb - fc)) +
-                    c * fa * fb / ((fc - fa) * (fc - fb))
-            else
-                # Secant method
-                s = b - fb * (b - a) / (fb - fa)
-            end
-            if (s < min((3 * a + b) / 4, b) || s > max((3 * a + b) / 4, b)) ||
-               (cond && abs(s - b) ≥ abs(b - c) / 2) ||
-               (!cond && abs(s - b) ≥ abs(c - d) / 2) ||
-               (cond && abs(b - c) ≤ ϵ) ||
-               (!cond && abs(c - d) ≤ ϵ)
-                # Bisection method
-                s = (a + b) / 2
-                (s == a || s == b) &&
-                    return SciMLBase.build_solution(prob, alg, a, fa;
-                        retcode = ReturnCode.FloatingPointLimit,
-                        left = a, right = b)
-                cond = true
-            else
-                cond = false
-            end
-            fs = f(s)
-            if abs((b - a) / 2) < abstol
-                return SciMLBase.build_solution(prob, alg, s, fs;
-                    retcode = ReturnCode.Success,
-                    left = a, right = b)
-            end
-            if iszero(fs)
-                if b < a
-                    a = b
-                    fa = fb
-                end
-                b = s
-                fb = fs
-                break
-            end
-            if fa * fs < 0
-                d = c
-                c = b
-                b = s
-                fb = fs
-            else
-                a = s
-                fa = fs
-            end
-            if abs(fa) < abs(fb)
-                d = c
-                c = b
-                b = a
-                a = c
-                fc = fb
-                fb = fa
-                fa = fc
-            end
-            i += 1
-        end
-    end
-
-    while i < maxiters
-        c = (a + b) / 2
-        if (c == a || c == b)
-            return SciMLBase.build_solution(prob, alg, a, fa;
-                retcode = ReturnCode.FloatingPointLimit,
-                left = a, right = b)
-        end
-        fc = f(c)
-        if abs((b - a) / 2) < abstol
-            return SciMLBase.build_solution(prob, alg, c, fc;
-                retcode = ReturnCode.Success,
-                left = a, right = b)
-        end
-        if iszero(fc)
-            b = c
-            fb = fc
-        else
-            a = c
-            fa = fc
-        end
-        i += 1
-    end
-
-    return SciMLBase.build_solution(prob, alg, a, fa; retcode = ReturnCode.MaxIters,
-        left = a, right = b)
-end
diff --git a/src/broyden.jl b/src/broyden.jl
deleted file mode 100644
index 07b2609..0000000
--- a/src/broyden.jl
+++ /dev/null
@@ -1,79 +0,0 @@
-"""
-    Broyden(; batched = false,
-        termination_condition = NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
-            abstol = nothing,
-            reltol = nothing))
-
-A low-overhead implementation of Broyden. This method is non-allocating on scalar
-and static array problems.
-
-!!! note
-
-    To use the `batched` version, remember to load `NNlib`, i.e., `using NNlib` or
-    `import NNlib` must be present in your code.
-"""
-struct Broyden{TC <: NLSolveTerminationCondition} <:
-       AbstractSimpleNonlinearSolveAlgorithm
-    termination_condition::TC
-end
-
-function Broyden(; batched = false,
-        termination_condition = NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
-            abstol = nothing,
-            reltol = nothing))
-    if batched
-        @assert NNlibExtLoaded[] "Please install and load `NNlib.jl` to use batched Broyden."
-        return BatchedBroyden(termination_condition)
-    end
-    return Broyden(termination_condition)
-end
-
-function SciMLBase.__solve(prob::NonlinearProblem, alg::Broyden, args...;
-        abstol = nothing, reltol = nothing, maxiters = 1000, kwargs...)
-    tc = alg.termination_condition
-    mode = DiffEqBase.get_termination_mode(tc)
-    f = Base.Fix2(prob.f, prob.p)
-    x = float(prob.u0)
-
-    fₙ = f(x)
-    T = eltype(x)
-    J⁻¹ = init_J(x)
-
-    if SciMLBase.isinplace(prob)
-        error("Broyden currently only supports out-of-place nonlinear problems")
-    end
-
-    atol = _get_tolerance(abstol, tc.abstol, T)
-    rtol = _get_tolerance(reltol, tc.reltol, T)
-
-    if mode ∈ DiffEqBase.SAFE_BEST_TERMINATION_MODES
-        error("Broyden currently doesn't support SAFE_BEST termination modes")
-    end
-
-    storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() :
-              nothing
-    termination_condition = tc(storage)
-
-    xₙ = x
-    xₙ₋₁ = x
-    fₙ₋₁ = fₙ
-    for _ in 1:maxiters
-        xₙ = xₙ₋₁ - _restructure(xₙ₋₁, J⁻¹ * _vec(fₙ₋₁))
-        fₙ = f(xₙ)
-        Δxₙ = xₙ - xₙ₋₁
-        Δfₙ = fₙ - fₙ₋₁
-        J⁻¹Δfₙ = _restructure(Δfₙ, J⁻¹ * _vec(Δfₙ))
-        J⁻¹ += _restructure(J⁻¹,
-            ((_vec(Δxₙ) .- _vec(J⁻¹Δfₙ)) ./ (_vec(Δxₙ)' * _vec(J⁻¹Δfₙ))) *
-            (_vec(Δxₙ)' * J⁻¹))
-
-        if termination_condition(fₙ, xₙ, xₙ₋₁, atol, rtol)
-            return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.Success)
-        end
-
-        xₙ₋₁ = xₙ
-        fₙ₋₁ = fₙ
-    end
-
-    return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.MaxIters)
-end
diff --git a/src/dfsane.jl b/src/dfsane.jl
deleted file mode 100644
index 49c50bc..0000000
--- a/src/dfsane.jl
+++ /dev/null
@@ -1,209 +0,0 @@
-"""
-    SimpleDFSane(; σ_min::Real = 1e-10, σ_max::Real = 1e10, σ_1::Real = 1.0,
-        M::Int = 10, γ::Real = 1e-4, τ_min::Real = 0.1, τ_max::Real = 0.5,
-        nexp::Int = 2, η_strategy::Function = (f_1, k, x, F) -> f_1 ./ k^2,
-        termination_condition = NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
-            abstol = nothing,
-            reltol = nothing),
-        batched::Bool = false,
-        max_inner_iterations::Int = 1000)
-
-A low-overhead implementation of the df-sane method for solving large-scale nonlinear
-systems of equations. For in depth information about all the parameters and the algorithm,
-see the paper: [W LaCruz, JM Martinez, and M Raydan (2006), Spectral residual mathod without
-gradient information for solving large-scale nonlinear systems of equations, Mathematics of
-Computation, 75, 1429-1448.](https://www.researchgate.net/publication/220576479_Spectral_Residual_Method_without_Gradient_Information_for_Solving_Large-Scale_Nonlinear_Systems_of_Equations)
-
-### Keyword Arguments
-
-- `σ_min`: the minimum value of the spectral coefficient `σ_k` which is related to the step
-  size in the algorithm. Defaults to `1e-10`.
-- `σ_max`: the maximum value of the spectral coefficient `σ_k` which is related to the step
-  size in the algorithm. Defaults to `1e10`.
-- `σ_1`: the initial value of the spectral coefficient `σ_k` which is related to the step
-  size in the algorithm.. Defaults to `1.0`.
-- `M`: The monotonicity of the algorithm is determined by a this positive integer.
-  A value of 1 for `M` would result in strict monotonicity in the decrease of the L2-norm
-  of the function `f`. However, higher values allow for more flexibility in this reduction.
-  Despite this, the algorithm still ensures global convergence through the use of a
-  non-monotone line-search algorithm that adheres to the Grippo-Lampariello-Lucidi
-  condition. Values in the range of 5 to 20 are usually sufficient, but some cases may call
-  for a higher value of `M`. The default setting is 10.
-- `γ`: a parameter that influences if a proposed step will be accepted. Higher value of `γ`
-  will make the algorithm more restrictive in accepting steps. Defaults to `1e-4`.
-- `τ_min`: if a step is rejected the new step size will get multiplied by factor, and this
-  parameter is the minimum value of that factor. Defaults to `0.1`.
-- `τ_max`: if a step is rejected the new step size will get multiplied by factor, and this
-  parameter is the maximum value of that factor. Defaults to `0.5`.
-- `nexp`: the exponent of the loss, i.e. ``f_k=||F(x_k)||^{nexp}``. The paper uses
-  `nexp ∈ {1,2}`. Defaults to `2`.
-- `η_strategy`:  function to determine the parameter `η_k`, which enables growth
-  of ``||F||^2``. Called as ``η_k = η_strategy(f_1, k, x, F)`` with `f_1` initialized as
-  ``f_1=||F(x_1)||^{nexp}``, `k` is the iteration number, `x` is the current `x`-value and
-  `F` the current residual. Should satisfy ``η_k > 0`` and ``∑ₖ ηₖ < ∞``. Defaults to
-  ``||F||^2 / k^2``.
-- `termination_condition`: a `NLSolveTerminationCondition` that determines when the solver
-  should terminate. Defaults to `NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
-  abstol = nothing, reltol = nothing)`.
-- `batched`: if `true`, the algorithm will use a batched version of the algorithm that treats each
-  column of `x` as a separate problem. This can be useful nonlinear problems involing neural
-  networks. Defaults to `false`.
-- `max_inner_iterations`: the maximum number of iterations allowed for the inner loop of the
-  algorithm. Used exclusively in `batched` mode. Defaults to `1000`.
-"""
-struct SimpleDFSane{T, TC} <: AbstractSimpleNonlinearSolveAlgorithm
-    σ_min::T
-    σ_max::T
-    σ_1::T
-    M::Int
-    γ::T
-    τ_min::T
-    τ_max::T
-    nexp::Int
-    η_strategy::Function
-    termination_condition::TC
-end
-
-function SimpleDFSane(; σ_min::Real = 1e-10, σ_max::Real = 1e10, σ_1::Real = 1.0,
-        M::Int = 10, γ::Real = 1e-4, τ_min::Real = 0.1, τ_max::Real = 0.5,
-        nexp::Int = 2, η_strategy::Function = (f_1, k, x, F) -> f_1 ./ k^2,
-        termination_condition = NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
-            abstol = nothing,
-            reltol = nothing),
-        batched::Bool = false,
-        max_inner_iterations = 1000)
-    if batched
-        return BatchedSimpleDFSane(; σₘᵢₙ = σ_min,
-            σₘₐₓ = σ_max,
-            σ₁ = σ_1,
-            M,
-            γ,
-            τₘᵢₙ = τ_min,
-            τₘₐₓ = τ_max,
-            nₑₓₚ = nexp,
-            ηₛ = η_strategy,
-            termination_condition,
-            max_inner_iterations)
-    end
-    return SimpleDFSane{typeof(σ_min), typeof(termination_condition)}(σ_min,
-        σ_max,
-        σ_1,
-        M,
-        γ,
-        τ_min,
-        τ_max,
-        nexp,
-        η_strategy,
-        termination_condition)
-end
-
-function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleDFSane,
-        args...; abstol = nothing, reltol = nothing, maxiters = 1000,
-        kwargs...)
-    tc = alg.termination_condition
-    mode = DiffEqBase.get_termination_mode(tc)
-
-    f = Base.Fix2(prob.f, prob.p)
-    x = float(prob.u0)
-
-    T = eltype(x)
-    σ_min = float(alg.σ_min)
-    σ_max = float(alg.σ_max)
-    σ_k = float(alg.σ_1)
-
-    M = alg.M
-    γ = float(alg.γ)
-    τ_min = float(alg.τ_min)
-    τ_max = float(alg.τ_max)
-    nexp = alg.nexp
-    η_strategy = alg.η_strategy
-
-    if SciMLBase.isinplace(prob)
-        error("SimpleDFSane currently only supports out-of-place nonlinear problems")
-    end
-
-    atol = _get_tolerance(abstol, tc.abstol, T)
-    rtol = _get_tolerance(reltol, tc.reltol, T)
-
-    if mode ∈ DiffEqBase.SAFE_BEST_TERMINATION_MODES
-        error("SimpleDFSane currently doesn't support SAFE_BEST termination modes")
-    end
-
-    storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() :
-              nothing
-    termination_condition = tc(storage)
-
-    function ff(x)
-        F = f(x)
-        f_k = norm(F)^nexp
-        return f_k, F
-    end
-
-    function generate_history(f_k, M)
-        return fill(f_k, M)
-    end
-
-    f_k, F_k = ff(x)
-    α_1 = convert(T, 1.0)
-    f_1 = f_k
-    history_f_k = generate_history(f_k, M)
-
-    for k in 1:maxiters
-        # Spectral parameter range check
-        σ_k = sign(σ_k) * clamp(abs(σ_k), σ_min, σ_max)
-
-        # Line search direction
-        d = -σ_k .* F_k
-
-        η = η_strategy(f_1, k, x, F_k)
-        f̄ = maximum(history_f_k)
-        α_p = α_1
-        α_m = α_1
-        x_new = @. x + α_p * d
-
-        f_new, F_new = ff(x_new)
-
-        inner_iterations = 0
-        while true
-            inner_iterations += 1
-
-            criteria = f̄ + η - γ * α_p^2 * f_k
-            f_new ≤ criteria && break
-
-            α_tp = @. α_p^2 * f_k / (f_new + (2 * α_p - 1) * f_k)
-            x_new = @. x - α_m * d
-            f_new, F_new = ff(x_new)
-
-            f_new ≤ criteria && break
-
-            α_tm = @. α_m^2 * f_k / (f_new + (2 * α_m - 1) * f_k)
-            α_p = @. clamp(α_tp, τ_min * α_p, τ_max * α_p)
-            α_m = @. clamp(α_tm, τ_min * α_m, τ_max * α_m)
-            x_new = @. x + α_p * d
-            f_new, F_new = ff(x_new)
-        end
-
-        if termination_condition(F_new, x_new, x, atol, rtol)
-            return SciMLBase.build_solution(prob,
-                alg,
-                x_new,
-                F_new;
-                retcode = ReturnCode.Success)
-        end
-
-        # Update spectral parameter
-        s_k = @. x_new - x
-        y_k = @. F_new - F_k
-
-        σ_k = (s_k' * s_k) / (s_k' * y_k)
-
-        # Take step
-        x = x_new
-        F_k = F_new
-        f_k = f_new
-
-        # Store function value
-        history_f_k[k % M + 1] = f_new
-    end
-    return SciMLBase.build_solution(prob, alg, x, F_k; retcode = ReturnCode.MaxIters)
-end
diff --git a/src/falsi.jl b/src/falsi.jl
deleted file mode 100644
index eb2ea1f..0000000
--- a/src/falsi.jl
+++ /dev/null
@@ -1,86 +0,0 @@
-"""
-`Falsi`: A non-allocating regula falsi method
-"""
-struct Falsi <: AbstractBracketingAlgorithm end
-
-function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Falsi, args...;
-        maxiters = 1000, abstol = min(eps(prob.tspan[1]), eps(prob.tspan[2])),
-        kwargs...)
-    f = Base.Fix2(prob.f, prob.p)
-    left, right = prob.tspan
-    fl, fr = f(left), f(right)
-
-    if iszero(fl)
-        return SciMLBase.build_solution(prob, alg, left, fl;
-            retcode = ReturnCode.ExactSolutionLeft, left = left,
-            right = right)
-    elseif iszero(fr)
-        return SciMLBase.build_solution(prob, alg, right, fr;
-            retcode = ReturnCode.ExactSolutionRight, left = left,
-            right = right)
-    end
-
-    i = 1
-    if !iszero(fr)
-        while i < maxiters
-            if nextfloat_tdir(left, prob.tspan...) == right
-                return SciMLBase.build_solution(prob, alg, left, fl;
-                    retcode = ReturnCode.FloatingPointLimit,
-                    left = left, right = right)
-            end
-            mid = (fr * left - fl * right) / (fr - fl)
-            for i in 1:10
-                mid = max_tdir(left, prevfloat_tdir(mid, prob.tspan...), prob.tspan...)
-            end
-            if mid == right || mid == left
-                break
-            end
-            fm = f(mid)
-            if abs((right - left) / 2) < abstol
-                return SciMLBase.build_solution(prob, alg, mid, fm;
-                    retcode = ReturnCode.Success,
-                    left = left, right = right)
-            end
-            if iszero(fm)
-                right = mid
-                break
-            end
-            if sign(fl) == sign(fm)
-                fl = fm
-                left = mid
-            else
-                fr = fm
-                right = mid
-            end
-            i += 1
-        end
-    end
-
-    while i < maxiters
-        mid = (left + right) / 2
-        (mid == left || mid == right) &&
-            return SciMLBase.build_solution(prob, alg, left, fl;
-                retcode = ReturnCode.FloatingPointLimit,
-                left = left, right = right)
-        fm = f(mid)
-        if abs((right - left) / 2) < abstol
-            return SciMLBase.build_solution(prob, alg, mid, fm;
-                retcode = ReturnCode.Success,
-                left = left, right = right)
-        end
-        if iszero(fm)
-            right = mid
-            fr = fm
-        elseif sign(fm) == sign(fl)
-            left = mid
-            fl = fm
-        else
-            right = mid
-            fr = fm
-        end
-        i += 1
-    end
-
-    return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters,
-        left = left, right = right)
-end
diff --git a/src/halley.jl b/src/halley.jl
deleted file mode 100644
index 8107dde..0000000
--- a/src/halley.jl
+++ /dev/null
@@ -1,118 +0,0 @@
-"""
-```julia
-SimpleHalley(; chunk_size = Val{0}(), autodiff = Val{true}(),
-                                 diff_type = Val{:forward})
-```
-
-A low-overhead implementation of SimpleHalley's Method. This method is non-allocating on scalar
-and static array problems.
-
-!!! note
-
-    As part of the decreased overhead, this method omits some of the higher level error
-    catching of the other methods. Thus, to see better error messages, use one of the other
-    methods like `NewtonRaphson`
-
-### Keyword Arguments
-
-- `chunk_size`: the chunk size used by the internal ForwardDiff.jl automatic differentiation
-  system. This allows for multiple derivative columns to be computed simultaneously,
-  improving performance. Defaults to `0`, which is equivalent to using ForwardDiff.jl's
-  default chunk size mechanism. For more details, see the documentation for
-  [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/).
-- `autodiff`: whether to use forward-mode automatic differentiation for the Jacobian.
-  Note that this argument is ignored if an analytical Jacobian is passed; as that will be
-  used instead. Defaults to `Val{true}`, which means ForwardDiff.jl is used by default.
-  If `Val{false}`, then FiniteDiff.jl is used for finite differencing.
-- `diff_type`: the type of finite differencing used if `autodiff = false`. Defaults to
-  `Val{:forward}` for forward finite differences. For more details on the choices, see the
-  [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl) documentation.
-"""
-struct SimpleHalley{CS, AD, FDT} <: AbstractNewtonAlgorithm{CS, AD, FDT}
-    function SimpleHalley(; chunk_size = Val{0}(), autodiff = Val{true}(),
-            diff_type = Val{:forward})
-        new{SciMLBase._unwrap_val(chunk_size), SciMLBase._unwrap_val(autodiff),
-            SciMLBase._unwrap_val(diff_type)}()
-    end
-end
-
-function SciMLBase.__solve(prob::NonlinearProblem,
-        alg::SimpleHalley, args...; abstol = nothing,
-        reltol = nothing,
-        maxiters = 1000, kwargs...)
-    f = Base.Fix2(prob.f, prob.p)
-    x = float(prob.u0)
-    fx = f(x)
-    if isa(x, AbstractArray)
-        n = length(x)
-    end
-    T = typeof(x)
-
-    if SciMLBase.isinplace(prob)
-        error("SimpleHalley currently only supports out-of-place nonlinear problems")
-    end
-
-    atol = abstol !== nothing ? abstol :
-           real(oneunit(eltype(T))) * (eps(real(one(eltype(T)))))^(4 // 5)
-    rtol = reltol !== nothing ? reltol : eps(real(one(eltype(T))))^(4 // 5)
-
-    if x isa Number
-        xo = oftype(one(eltype(x)), Inf)
-    else
-        xo = map(x -> oftype(one(eltype(x)), Inf), x)
-    end
-
-    for i in 1:maxiters
-        if alg_autodiff(alg)
-            if isa(x, Number)
-                fx = f(x)
-                dfx = ForwardDiff.derivative(f, x)
-                d2fx = ForwardDiff.derivative(x -> ForwardDiff.derivative(f, x), x)
-            else
-                fx = f(x)
-                dfx = ForwardDiff.jacobian(f, x)
-                d2fx = ForwardDiff.jacobian(x -> ForwardDiff.jacobian(f, x), x)
-                ai = -(dfx \ fx)
-                A = reshape(d2fx * ai, (n, n))
-                bi = (dfx) \ (A * ai)
-                ci = (ai .* ai) ./ (ai .+ (0.5 .* bi))
-            end
-        else
-            if isa(x, Number)
-                fx = f(x)
-                dfx = FiniteDiff.finite_difference_derivative(f, x, diff_type(alg),
-                    eltype(x))
-                d2fx = FiniteDiff.finite_difference_derivative(x -> FiniteDiff.finite_difference_derivative(f,
-                        x),
-                    x,
-                    diff_type(alg), eltype(x))
-            else
-                fx = f(x)
-                dfx = FiniteDiff.finite_difference_jacobian(f, x, diff_type(alg), eltype(x))
-                d2fx = FiniteDiff.finite_difference_jacobian(x -> FiniteDiff.finite_difference_jacobian(f,
-                        x),
-                    x,
-                    diff_type(alg), eltype(x))
-                ai = -(dfx \ fx)
-                A = reshape(d2fx * ai, (n, n))
-                bi = (dfx) \ (A * ai)
-                ci = (ai .* ai) ./ (ai .+ (0.5 .* bi))
-            end
-        end
-        iszero(fx) &&
-            return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
-        if isa(x, Number)
-            Δx = (2 * dfx^2 - fx * d2fx) \ (2fx * dfx)
-            x -= Δx
-        else
-            Δx = ci
-            x += Δx
-        end
-        if isapprox(x, xo, atol = atol, rtol = rtol)
-            return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
-        end
-        xo = x
-    end
-
-    return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)
-end
diff --git a/src/itp.jl b/src/itp.jl
deleted file mode 100644
index 3147c52..0000000
--- a/src/itp.jl
+++ /dev/null
@@ -1,149 +0,0 @@
-"""
-```julia
-ITP(; k1::Real = 0.007, k2::Real = 1.5, n0::Int = 10)
-```
-
-ITP (Interpolate Truncate & Project)
-
-Use the [ITP method](https://en.wikipedia.org/wiki/ITP_method) to find
-a root of a bracketed function, with a convergence rate between 1 and 1.62.
-
-This method was introduced in the paper "An Enhancement of the Bisection Method
-Average Performance Preserving Minmax Optimality"
-(https://doi.org/10.1145/3423597) by I. F. D. Oliveira and R. H. C. Takahashi.
-
-# Tuning Parameters
-
-The following keyword parameters are accepted.
-
-- `n₀::Int = 1`, the 'slack'. Must not be negative.\n
-  When n₀ = 0 the worst-case is identical to that of bisection,
-  but increacing n₀ provides greater oppotunity for superlinearity.
-- `κ₁::Float64 = 0.1`. Must not be negative.\n
-  The recomended value is `0.2/(x₂ - x₁)`.
-  Lower values produce tighter asymptotic behaviour, while higher values
-  improve the steady-state behaviour when truncation is not helpful.
-- `κ₂::Real = 2`. Must lie in [1, 1+ϕ ≈ 2.62).\n
-  Higher values allow for a greater convergence rate,
-  but also make the method more succeptable to worst-case performance.
-  In practice, κ=1,2 seems to work well due to the computational simplicity,
-  as κ₂ is used as an exponent in the method.
-
-### Worst Case Performance
-
-n½ + `n₀` iterations, where n½ is the number of iterations using bisection
-(n½ = ⌈log2(Δx)/2`tol`⌉).
-
-### Asymptotic Performance
-
-If `f` is twice differentiable and the root is simple,
-then with `n₀` > 0 the convergence rate is √`κ₂`.
-"""
-struct ITP{T} <: AbstractBracketingAlgorithm
-    k1::T
-    k2::T
-    n0::Int
-    function ITP(; k1::Real = 0.007, k2::Real = 1.5, n0::Int = 10)
-        if k1 < 0
-            error("Hyper-parameter κ₁ should not be negative")
-        end
-        if n0 < 0
-            error("Hyper-parameter n₀ should not be negative")
-        end
-        if k2 < 1 || k2 > (1.5 + sqrt(5) / 2)
-            ArgumentError("Hyper-parameter κ₂ should be between 1 and 1 + ϕ where ϕ ≈ 1.618... is the golden ratio")
-        end
-        T = promote_type(eltype(k1), eltype(k2))
-        return new{T}(k1, k2, n0)
-    end
-end
-
-function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::ITP,
-        args...; abstol = min(eps(prob.tspan[1]), eps(prob.tspan[2])),
-        maxiters = 1000, kwargs...)
-    f = Base.Fix2(prob.f, prob.p)
-    left, right = prob.tspan # a and b
-    fl, fr = f(left), f(right)
-    ϵ = abstol
-    if iszero(fl)
-        return SciMLBase.build_solution(prob, alg, left, fl;
-            retcode = ReturnCode.ExactSolutionLeft, left = left,
-            right = right)
-    elseif iszero(fr)
-        return SciMLBase.build_solution(prob, alg, right, fr;
-            retcode = ReturnCode.ExactSolutionRight, left = left,
-            right = right)
-    end
-    #defining variables/cache
-    k1 = alg.k1
-    k2 = alg.k2
-    n0 = alg.n0
-    n_h = ceil(log2(abs(right - left) / (2 * ϵ)))
-    mid = (left + right) / 2
-    x_f = left + (right - left) * (fl / (fl - fr))
-    xt = left
-    xp = left
-    r = zero(left) #minmax radius
-    δ = zero(left) # truncation error
-    σ = 1.0
-    ϵ_s = ϵ * 2^(n_h + n0)
-    i = 0 #iteration
-    while i <= maxiters
-        span = abs(right - left)
-        r = ϵ_s - (span / 2)
-        δ = k1 * (span^k2)
-
-        ## Interpolation step ##
-        x_f = left + (right - left) * (fl / (fl - fr))
-
-        ## Truncation step ##
-        σ = sign(mid - x_f)
-        if δ <= abs(mid - x_f)
-            xt = x_f + (σ * δ)
-        else
-            xt = mid
-        end
-
-        ## Projection step ##
-        if abs(xt - mid) <= r
-            xp = xt
-        else
-            xp = mid - (σ * r)
-        end
-
-        if abs((left - right) / 2) < ϵ
-            return SciMLBase.build_solution(prob, alg, mid, f(mid);
-                retcode = ReturnCode.Success,
-                left = left, right = right)
-        end
-
-        ## Update ##
-        tmin, tmax = minmax(left, right)
-        xp >= tmax && (xp = prevfloat(tmax))
-        xp <= tmin && (xp = nextfloat(tmin))
-        yp = f(xp)
-        yps = yp * sign(fr)
-        if yps > 0
-            right = xp
-            fr = yp
-        elseif yps < 0
-            left = xp
-            fl = yp
-        else
-            return SciMLBase.build_solution(prob, alg, xp, yps;
-                retcode = ReturnCode.Success, left = xp,
-                right = xp)
-        end
-        i += 1
-        mid = (left + right) / 2
-        ϵ_s /= 2
-
-        if nextfloat_tdir(left, prob.tspan...) == right
-            return SciMLBase.build_solution(prob, alg, left, fl;
-                retcode = ReturnCode.FloatingPointLimit, left = left,
-                right = right)
-        end
-    end
-    return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters,
-        left = left, right = right)
-end
diff --git a/src/klement.jl b/src/klement.jl
deleted file mode 100644
index e6a38ec..0000000
--- a/src/klement.jl
+++ /dev/null
@@ -1,106 +0,0 @@
-"""
-```julia
-Klement()
-```
-
-A low-overhead implementation of [Klement](https://jatm.com.br/jatm/article/view/373).
-This method is non-allocating on scalar problems.
-"""
-struct Klement <: AbstractSimpleNonlinearSolveAlgorithm end
-
-function SciMLBase.__solve(prob::NonlinearProblem,
-        alg::Klement, args...; abstol = nothing,
-        reltol = nothing,
-        maxiters = 1000, kwargs...)
-    f = Base.Fix2(prob.f, prob.p)
-    x = float(prob.u0)
-    fₙ = f(x)
-    T = eltype(x)
-    singular_tol = 1e-9
-
-    if SciMLBase.isinplace(prob)
-        error("Klement currently only supports out-of-place nonlinear problems")
-    end
-
-    atol = abstol !== nothing ? abstol :
-           real(oneunit(eltype(T))) * (eps(real(one(eltype(T)))))^(4 // 5)
-    rtol = reltol !== nothing ? reltol : eps(real(one(eltype(T))))^(4 // 5)
-
-    xₙ = x
-    xₙ₋₁ = x
-    fₙ₋₁ = fₙ
-
-    # x is scalar
-    if x isa Number
-        J = 1.0
-        for _ in 1:maxiters
-            xₙ = xₙ₋₁ - fₙ₋₁ / J
-            fₙ = f(xₙ)
-
-            iszero(fₙ) &&
-                return SciMLBase.build_solution(prob, alg, xₙ, fₙ;
-                    retcode = ReturnCode.Success)
-
-            if isapprox(xₙ, xₙ₋₁, atol = atol, rtol = rtol)
-                return SciMLBase.build_solution(prob, alg, xₙ, fₙ;
-                    retcode = ReturnCode.Success)
-            end
-
-            Δxₙ = xₙ - xₙ₋₁
-            Δfₙ = fₙ - fₙ₋₁
-
-            # Prevent division by 0
-            denominator = max(J^2 * Δxₙ^2, 1e-9)
-
-            k = (Δfₙ - J * Δxₙ) / denominator
-            J += (k * Δxₙ * J) * J
-
-            # Singularity test
-            if J < singular_tol
-                J = 1.0
-            end
-
-            xₙ₋₁ = xₙ
-            fₙ₋₁ = fₙ
-        end
-        # x is a vector
-    else
-        J = init_J(x)
-        for _ in 1:maxiters
-            F = lu(J, check = false)
-
-            # Singularity test
-            if any(abs.(F.U[diagind(F.U)]) .< singular_tol)
-                J = init_J(xₙ)
-                F = lu(J, check = false)
-            end
-
-            tmp = _restructure(fₙ₋₁, F \ _vec(fₙ₋₁))
-            xₙ = xₙ₋₁ - tmp
-            fₙ = f(xₙ)
-
-            iszero(fₙ) &&
-                return SciMLBase.build_solution(prob, alg, xₙ, fₙ;
-                    retcode = ReturnCode.Success)
-
-            if isapprox(xₙ, xₙ₋₁, atol = atol, rtol = rtol)
-                return SciMLBase.build_solution(prob, alg, xₙ, fₙ;
-                    retcode = ReturnCode.Success)
-            end
-
-            Δxₙ = xₙ - xₙ₋₁
-            Δfₙ = fₙ - fₙ₋₁
-
-            # Prevent division by 0
-            denominator = _restructure(Δxₙ, max.(J' .^ 2 * _vec(Δxₙ) .^ 2, 1e-9))
-
-            k = (Δfₙ - _restructure(Δxₙ, J * _vec(Δxₙ))) ./ denominator
-            J += (_vec(k) * _vec(Δxₙ)' .* J) * J
-
-            xₙ₋₁ = xₙ
-            fₙ₋₁ = fₙ
-        end
-    end
-
-    return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.MaxIters)
-end
diff --git a/src/lbroyden.jl b/src/lbroyden.jl
deleted file mode 100644
index 4820921..0000000
--- a/src/lbroyden.jl
+++ /dev/null
@@ -1,144 +0,0 @@
-"""
-    LBroyden(; batched = false,
-              termination_condition = NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
-                                                                  abstol = nothing, reltol = nothing),
-              threshold::Int = 27)
-
-A limited memory implementation of Broyden. This method applies the L-BFGS scheme to
-Broyden's method.
-
-!!! warn
-
-    This method is not very stable and can diverge even for very simple problems. This has mostly been
-    tested for neural networks in DeepEquilibriumNetworks.jl.
-"""
-struct LBroyden{batched, TC <: NLSolveTerminationCondition} <:
-       AbstractSimpleNonlinearSolveAlgorithm
-    termination_condition::TC
-    threshold::Int
-
-    function LBroyden(; batched = false, threshold::Int = 27,
-            termination_condition = NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
-                abstol = nothing,
-                reltol = nothing))
-        return new{batched, typeof(termination_condition)}(termination_condition, threshold)
-    end
-end
-
-@views function SciMLBase.__solve(prob::NonlinearProblem, alg::LBroyden{batched}, args...;
-        abstol = nothing, reltol = nothing, maxiters = 1000,
-        kwargs...) where {batched}
-    tc = alg.termination_condition
-    mode = DiffEqBase.get_termination_mode(tc)
-    threshold = min(maxiters, alg.threshold)
-    x = float(prob.u0)
-
-    batched && @assert ndims(x)==2 "Batched LBroyden only supports 2D arrays"
-
-    if x isa Number
-        restore_scalar = true
-        x = [x]
-        f = u -> prob.f(u[], prob.p)
-    else
-        f = Base.Fix2(prob.f, prob.p)
-        restore_scalar = false
-    end
-
-    fₙ = f(x)
-    T = eltype(x)
-
-    if SciMLBase.isinplace(prob)
-        error("LBroyden currently only supports out-of-place nonlinear problems")
-    end
-
-    U, Vᵀ = _init_lbroyden_state(batched, x, threshold)
-
-    atol = abstol !== nothing ? abstol :
-           (tc.abstol !== nothing ? tc.abstol :
-            real(oneunit(eltype(T))) * (eps(real(one(eltype(T)))))^(4 // 5))
-    rtol = reltol !== nothing ? reltol :
-           (tc.reltol !== nothing ? tc.reltol : eps(real(one(eltype(T))))^(4 // 5))
-
-    if mode ∈ DiffEqBase.SAFE_BEST_TERMINATION_MODES
-        error("LBroyden currently doesn't support SAFE_BEST termination modes")
-    end
-
-    storage = mode ∈ DiffEqBase.SAFE_TERMINATION_MODES ? NLSolveSafeTerminationResult() :
-              nothing
-    termination_condition = tc(storage)
-
-    xₙ = x
-    xₙ₋₁ = x
-    fₙ₋₁ = fₙ
-    update = fₙ
-    for i in 1:maxiters
-        xₙ = xₙ₋₁ .+ update
-        fₙ = f(xₙ)
-        Δxₙ = xₙ .- xₙ₋₁
-        Δfₙ = fₙ .- fₙ₋₁
-
-        if termination_condition(restore_scalar ? [fₙ] : fₙ, xₙ, xₙ₋₁, atol, rtol)
-            xₙ = restore_scalar ? xₙ[] : xₙ
-            return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.Success)
-        end
-
-        _U = selectdim(U, 1, 1:min(threshold, i))
-        _Vᵀ = selectdim(Vᵀ, 2, 1:min(threshold, i))
-
-        vᵀ = _rmatvec(_U, _Vᵀ, Δxₙ)
-        mvec = _matvec(_U, _Vᵀ, Δfₙ)
-        u = (Δxₙ .- mvec) ./ (sum(vᵀ .* Δfₙ) .+ convert(T, 1e-5))
-
-        selectdim(Vᵀ, 2, mod1(i, threshold)) .= vᵀ
-        selectdim(U, 1, mod1(i, threshold)) .= u
-
-        update = -_matvec(selectdim(U, 1, 1:min(threshold, i + 1)),
-            selectdim(Vᵀ, 2, 1:min(threshold, i + 1)), fₙ)
-
-        xₙ₋₁ = xₙ
-        fₙ₋₁ = fₙ
-    end
-
-    xₙ = restore_scalar ? xₙ[] : xₙ
-    return SciMLBase.build_solution(prob, alg, xₙ, fₙ; retcode = ReturnCode.MaxIters)
-end
-
-function _init_lbroyden_state(batched::Bool, x, threshold)
-    T = eltype(x)
-    if batched
-        U = fill!(similar(x, (threshold, size(x, 1), size(x, 2))), zero(T))
-        Vᵀ = fill!(similar(x, (size(x, 1), threshold, size(x, 2))), zero(T))
-    else
-        U = fill!(similar(x, (threshold, length(x))), zero(T))
-        Vᵀ = fill!(similar(x, (length(x), threshold)), zero(T))
-    end
-    return U, Vᵀ
-end
-
-function _rmatvec(U::AbstractMatrix, Vᵀ::AbstractMatrix,
-        x::Union{<:AbstractVector, <:Number})
-    length(U) == 0 && return x
-    return -x .+ vec((x' * Vᵀ) * U)
-end
-
-function _rmatvec(U::AbstractArray{T1, 3}, Vᵀ::AbstractArray{T2, 3},
-        x::AbstractMatrix) where {T1, T2}
-    length(U) == 0 && return x
-    Vᵀx = sum(Vᵀ .* reshape(x, size(x, 1), 1, size(x, 2)); dims = 1)
-    return -x .+ _drdims_sum(U .* permutedims(Vᵀx, (2, 1, 3)); dims = 1)
-end
-
-function _matvec(U::AbstractMatrix, Vᵀ::AbstractMatrix,
-        x::Union{<:AbstractVector, <:Number})
-    length(U) == 0 && return x
-    return -x .+ vec(Vᵀ * (U * x))
-end
-
-function _matvec(U::AbstractArray{T1, 3}, Vᵀ::AbstractArray{T2, 3},
-        x::AbstractMatrix) where {T1, T2}
-    length(U) == 0 && return x
-    xUᵀ = sum(reshape(x, size(x, 1), 1, size(x, 2)) .* permutedims(U, (2, 1, 3)); dims = 1)
-    return -x .+ _drdims_sum(xUᵀ .* Vᵀ; dims = 2)
-end
-
-_drdims_sum(args...; dims = :) = dropdims(sum(args...; dims); dims)
diff --git a/src/nlsolve/broyden.jl b/src/nlsolve/broyden.jl
new file mode 100644
index 0000000..1e544a4
--- /dev/null
+++ b/src/nlsolve/broyden.jl
@@ -0,0 +1,56 @@
+"""
+    SimpleBroyden()
+
+A low-overhead implementation of Broyden. This method is non-allocating on scalar
+and static array problems.
+"""
+struct SimpleBroyden <: AbstractSimpleNonlinearSolveAlgorithm end
+
+function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleBroyden, args...;
+        abstol = nothing, reltol = nothing, maxiters = 1000,
+        termination_condition = nothing, kwargs...)
+    @bb x = copy(float(prob.u0))
+    fx = _get_fx(prob, x)
+
+    @bb xo = copy(x)
+    @bb δx = copy(x)
+    @bb δf = copy(fx)
+    @bb fprev = copy(fx)
+
+    J⁻¹ = __init_identity_jacobian(fx, x)
+    @bb J⁻¹δf = copy(x)
+    @bb xᵀJ⁻¹ = copy(x)
+    @bb δJ⁻¹n = copy(x)
+    @bb δJ⁻¹ = copy(J⁻¹)
+
+    abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fx, x,
+        termination_condition)
+
+    for _ in 1:maxiters
+        @bb δx = J⁻¹ × vec(fprev)
+        @bb @. x = xo - δx
+        fx = __eval_f(prob, fx, x)
+        @bb @. δf = fx - fprev
+
+        # Termination Checks
+        tc_sol = check_termination(tc_cache, fx, x, xo, prob, alg)
+        tc_sol !== nothing && return tc_sol
+
+        @bb J⁻¹δf = J⁻¹ × vec(δf)
+        @bb δx .*= -1
+        d = dot(δx, J⁻¹δf)
+        @bb xᵀJ⁻¹ = transpose(J⁻¹) × vec(δx)
+
+        @bb @. δJ⁻¹n = (δx - J⁻¹δf) / d
+
+        δJ⁻¹n_ = _vec(δJ⁻¹n)
+        xᵀJ⁻¹_ = _vec(xᵀJ⁻¹)
+        @bb δJ⁻¹ = δJ⁻¹n_ × transpose(xᵀJ⁻¹_)
+        @bb J⁻¹ .+= δJ⁻¹
+
+        @bb copyto!(xo, x)
+        @bb copyto!(fprev, fx)
+    end
+
+    return build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)
+end
diff --git a/src/nlsolve/dfsane.jl b/src/nlsolve/dfsane.jl
new file mode 100644
index 0000000..2cbbd16
--- /dev/null
+++ b/src/nlsolve/dfsane.jl
@@ -0,0 +1,146 @@
+"""
+    SimpleDFSane(; σ_min::Real = 1e-10, σ_max::Real = 1e10, σ_1::Real = 1.0,
+        M::Int = 10, γ::Real = 1e-4, τ_min::Real = 0.1, τ_max::Real = 0.5,
+        nexp::Int = 2, η_strategy::Function = (f_1, k, x, F) -> f_1 ./ k^2)
+
+A low-overhead implementation of the df-sane method for solving large-scale nonlinear
+systems of equations. For in depth information about all the parameters and the algorithm,
+see the paper: [W LaCruz, JM Martinez, and M Raydan (2006), Spectral residual mathod without
+gradient information for solving large-scale nonlinear systems of equations, Mathematics of
+Computation, 75, 1429-1448.](https://www.researchgate.net/publication/220576479_Spectral_Residual_Method_without_Gradient_Information_for_Solving_Large-Scale_Nonlinear_Systems_of_Equations)
+
+### Keyword Arguments
+
+  - `σ_min`: the minimum value of the spectral coefficient `σ_k` which is related to the step
+    size in the algorithm. Defaults to `1e-10`.
+  - `σ_max`: the maximum value of the spectral coefficient `σ_k` which is related to the step
+    size in the algorithm. Defaults to `1e10`.
+  - `σ_1`: the initial value of the spectral coefficient `σ_k` which is related to the step
+    size in the algorithm.. Defaults to `1.0`.
+  - `M`: The monotonicity of the algorithm is determined by a this positive integer.
+    A value of 1 for `M` would result in strict monotonicity in the decrease of the L2-norm
+    of the function `f`. However, higher values allow for more flexibility in this reduction.
+    Despite this, the algorithm still ensures global convergence through the use of a
+    non-monotone line-search algorithm that adheres to the Grippo-Lampariello-Lucidi
+    condition. Values in the range of 5 to 20 are usually sufficient, but some cases may call
+    for a higher value of `M`. The default setting is 10.
+  - `γ`: a parameter that influences if a proposed step will be accepted. Higher value of `γ`
+    will make the algorithm more restrictive in accepting steps. Defaults to `1e-4`.
+  - `τ_min`: if a step is rejected the new step size will get multiplied by factor, and this
+    parameter is the minimum value of that factor. Defaults to `0.1`.
+  - `τ_max`: if a step is rejected the new step size will get multiplied by factor, and this
+    parameter is the maximum value of that factor. Defaults to `0.5`.
+  - `nexp`: the exponent of the loss, i.e. ``f_k=||F(x_k)||^{nexp}``. The paper uses
+    `nexp ∈ {1,2}`. Defaults to `2`.
+  - `η_strategy`:  function to determine the parameter `η_k`, which enables growth
+    of ``||F||^2``. Called as ``η_k = η_strategy(f_1, k, x, F)`` with `f_1` initialized as
+    ``f_1=||F(x_1)||^{nexp}``, `k` is the iteration number, `x` is the current `x`-value and
+    `F` the current residual. Should satisfy ``η_k > 0`` and ``∑ₖ ηₖ < ∞``. Defaults to
+    ``||F||^2 / k^2``.
+"""
+@kwdef @concrete struct SimpleDFSane <: AbstractSimpleNonlinearSolveAlgorithm
+    σ_min = 1e-10
+    σ_max = 1e10
+    σ_1 = 1.0
+    M::Int = 10
+    γ = 1e-4
+    τ_min = 0.1
+    τ_max = 0.5
+    nexp::Int = 2
+    η_strategy = (f_1, k, x, F) -> f_1 ./ k^2
+end
+
+function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleDFSane, args...;
+        abstol = nothing, reltol = nothing, maxiters = 1000,
+        termination_condition = nothing, kwargs...)
+    x = float(copy(prob.u0))
+    fx = _get_fx(prob, x)
+    T = eltype(x)
+
+    σ_min = T(alg.σ_min)
+    σ_max = T(alg.σ_max)
+    σ_k = T(alg.σ_1)
+
+    (; M, nexp, η_strategy) = alg
+    γ = T(alg.γ)
+    τ_min = T(alg.τ_min)
+    τ_max = T(alg.τ_max)
+
+    abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fx, x,
+        termination_condition)
+
+    fx_norm = NONLINEARSOLVE_DEFAULT_NORM(fx)^nexp
+    α_1 = one(T)
+    f_1 = fx_norm
+
+    history_f_k = fill(fx_norm, M)
+
+    # Generate the cache
+    @bb d = copy(x)
+    @bb xo = copy(x)
+    @bb x_cache = copy(x)
+    @bb δx = copy(x)
+    @bb fxo = copy(fx)
+    @bb δf = copy(fx)
+
+    k = 0
+    while k < maxiters
+        # Spectral parameter range check
+        σ_k = sign(σ_k) * clamp(abs(σ_k), σ_min, σ_max)
+
+        # Line search direction
+        @bb @. d = -σ_k * fx
+
+        η = η_strategy(f_1, k, x, fx)
+        f_bar = maximum(history_f_k)
+        α_p = α_1
+        α_m = α_1
+
+        @bb @. x += α_p * d
+
+        fx = __eval_f(prob, fx, x)
+        fx_norm_new = NONLINEARSOLVE_DEFAULT_NORM(fx)^nexp
+
+        while k < maxiters
+            fx_norm_new ≤ (f_bar + η - γ * α_p^2 * fx_norm) && break
+
+            α_p = α_p^2 * fx_norm / (fx_norm_new + (T(2) * α_p - T(1)) * fx_norm)
+            @bb @. x -= α_m * d
+
+            fx = __eval_f(prob, fx, x)
+            fx_norm_new = NONLINEARSOLVE_DEFAULT_NORM(fx)^nexp
+
+            fx_norm_new ≤ (f_bar + η - γ * α_m^2 * fx_norm) && break
+
+            α_tm = α_m^2 * fx_norm / (fx_norm_new + (T(2) * α_m - T(1)) * fx_norm)
+            α_p = clamp(α_p, τ_min * α_p, τ_max * α_p)
+            α_m = clamp(α_tm, τ_min * α_m, τ_max * α_m)
+            @bb @. x += α_p * d
+
+            fx = __eval_f(prob, fx, x)
+            fx_norm_new = NONLINEARSOLVE_DEFAULT_NORM(fx)^nexp
+
+            k += 1
+        end
+
+        tc_sol = check_termination(tc_cache, fx, x, xo, prob, alg)
+        tc_sol !== nothing && return tc_sol
+
+        # Update spectral parameter
+        @bb @. δx = x - xo
+        @bb @. δf = fx - fxo
+
+        σ_k = dot(δx, δx) / dot(δx, δf)
+
+        # Take step
+        @bb copyto!(xo, x)
+        @bb copyto!(fxo, fx)
+        fx_norm = fx_norm_new
+
+        # Store function value
+        history_f_k[mod1(k, M)] = fx_norm_new
+        k += 1
+    end
+
+    return build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)
+end
diff --git a/src/nlsolve/halley.jl b/src/nlsolve/halley.jl
new file mode 100644
index 0000000..161abae
--- /dev/null
+++ b/src/nlsolve/halley.jl
@@ -0,0 +1,81 @@
+"""
+    SimpleHalley(autodiff)
+    SimpleHalley(; autodiff = AutoForwardDiff())
+
+A low-overhead implementation of Halley's Method.
+
+!!! note
+
+    As part of the decreased overhead, this method omits some of the higher level error
+    catching of the other methods. Thus, to see better error messages, use one of the other
+    methods like `NewtonRaphson`
+
+### Keyword Arguments
+
+  - `autodiff`: determines the backend used for the Hessian. Defaults to
+    `AutoForwardDiff()`. Valid choices are `AutoForwardDiff()` or `AutoFiniteDiff()`.
+"""
+@kwdef @concrete struct SimpleHalley <: AbstractNewtonAlgorithm
+    autodiff = AutoForwardDiff()
+end
+
+function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleHalley, args...;
+        abstol = nothing, reltol = nothing, maxiters = 1000,
+        termination_condition = nothing, kwargs...)
+    isinplace(prob) &&
+        error("SimpleHalley currently only supports out-of-place nonlinear problems")
+
+    x = copy(float(prob.u0))
+    fx = _get_fx(prob, x)
+    T = eltype(x)
+
+    abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fx, x,
+        termination_condition)
+
+    @bb xo = copy(x)
+
+    if setindex_trait(x) === CanSetindex()
+        A = similar(x, length(x), length(x))
+        Aaᵢ = similar(x, length(x))
+        cᵢ = similar(x)
+    else
+        A = x
+        Aaᵢ = x
+        cᵢ = x
+    end
+
+    for i in 1:maxiters
+        # Hessian Computation is unfortunately type unstable
+        fx, dfx, d2fx = compute_jacobian_and_hessian(alg.autodiff, prob, fx, x)
+        setindex_trait(x) === CannotSetindex() && (A = dfx)
+
+        aᵢ = dfx \ _vec(fx)
+        A_ = _vec(A)
+        @bb A_ = d2fx × aᵢ
+        A = _restructure(A, A_)
+
+        @bb Aaᵢ = A × aᵢ
+        @bb A .*= -1
+        bᵢ = dfx \ Aaᵢ
+
+        cᵢ_ = _vec(cᵢ)
+        @bb @. cᵢ_ = (aᵢ * aᵢ) / (-aᵢ + (T(0.5) * bᵢ))
+        cᵢ = _restructure(cᵢ, cᵢ_)
+
+        if i == 1
+            if iszero(fx)
+                return build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
+            end
+        else
+            # Termination Checks
+            tc_sol = check_termination(tc_cache, fx, x, xo, prob, alg)
+            tc_sol !== nothing && return tc_sol
+        end
+
+        @bb @. x += cᵢ
+
+        @bb copyto!(xo, x)
+    end
+
+    return build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)
+end
diff --git a/src/nlsolve/klement.jl b/src/nlsolve/klement.jl
new file mode 100644
index 0000000..5041dc4
--- /dev/null
+++ b/src/nlsolve/klement.jl
@@ -0,0 +1,90 @@
+"""
+    SimpleKlement()
+
+A low-overhead implementation of [Klement](https://jatm.com.br/jatm/article/view/373). This
+method is non-allocating on scalar and static array problems.
+"""
+struct SimpleKlement <: AbstractSimpleNonlinearSolveAlgorithm end
+
+function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleKlement, args...;
+        abstol = nothing, reltol = nothing, maxiters = 1000,
+        termination_condition = nothing, kwargs...)
+    x = float(prob.u0)
+    T = eltype(x)
+    fx = _get_fx(prob, x)
+
+    singular_tol = eps(T)^(2 // 3)
+
+    abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fx, x,
+        termination_condition)
+
+    @bb δx = copy(x)
+    @bb fprev = copy(fx)
+    @bb xo = copy(x)
+    @bb δf = copy(fx)
+    @bb d = copy(x)
+
+    J = __init_identity_jacobian(fx, x)
+    @bb J_cache = copy(J)
+    @bb δx² = copy(x)
+    @bb J_cache2 = copy(J)
+    @bb F = copy(J)
+
+    for _ in 1:maxiters
+        if x isa Number
+            J < singular_tol && (J = __init_identity_jacobian!!(J))
+            F_ = J
+        else
+            @bb copyto!(F, J)
+            if setindex_trait(F) === CanSetindex()
+                F_ = lu!(F; check = false)
+            else
+                F_ = lu(F; check = false)
+            end
+
+            # Singularity test
+            if !issuccess(F_)
+                J = __init_identity_jacobian!!(J)
+                if setindex_trait(J) === CanSetindex()
+                    lu!(J; check = false)
+                else
+                    J = lu(J; check = false)
+                end
+            end
+        end
+
+        @bb copyto!(δx, fprev)
+        if setindex_trait(δx) === CanSetindex()
+            ldiv!(F_, _vec(δx))
+        else
+            δx = _restructure(δx, F_ \ _vec(δx))
+        end
+        @bb @. x = xo - δx
+        fx = __eval_f(prob, fx, x)
+
+        # Termination Checks
+        tc_sol = check_termination(tc_cache, fx, x, xo, prob, alg)
+        tc_sol !== nothing && return tc_sol
+
+        @bb δx .*= -1
+        @bb @. δf = fx - fprev
+
+        # Prevent division by 0
+        @bb @. δx² = δx^2
+        @bb @. J_cache = J^2
+        @bb d = transpose(J_cache) × vec(δx²)
+        @bb @. d = max(d, singular_tol)
+
+        @bb δx² = J × vec(δx)
+        @bb @. δf = (δf - δx²) / d
+
+        _vδf, _vδx = _vec(δf), _vec(δx)
+        @bb J_cache = _vδf × transpose(_vδx)
+        @bb @. J_cache *= J
+        @bb J_cache2 = J_cache × J
+
+        @bb @. J += J_cache2
+    end
+
+    return build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)
+end
diff --git a/src/nlsolve/lbroyden.jl b/src/nlsolve/lbroyden.jl
new file mode 100644
index 0000000..4cc8ee0
--- /dev/null
+++ b/src/nlsolve/lbroyden.jl
@@ -0,0 +1,119 @@
+"""
+    SimpleLimitedMemoryBroyden(; threshold::Int = 27)
+    SimpleLimitedMemoryBroyden(; threshold::Val = Val(27))
+
+A limited memory implementation of Broyden. This method applies the L-BFGS scheme to
+Broyden's method. This Alogrithm unfortunately cannot non-allocating for StaticArrays
+without compromising on the "simple" aspect.
+
+If the threshold is larger than the problem size, then this method will use `SimpleBroyden`.
+
+!!! warning
+
+    This method is not very stable and can diverge even for very simple problems. This has
+    mostly been tested for neural networks in DeepEquilibriumNetworks.jl.
+"""
+struct SimpleLimitedMemoryBroyden{threshold} <: AbstractSimpleNonlinearSolveAlgorithm end
+
+__get_threshold(::SimpleLimitedMemoryBroyden{threshold}) where {threshold} = Val(threshold)
+
+function SimpleLimitedMemoryBroyden(; threshold::Union{Val, Int} = Val(27))
+    return SimpleLimitedMemoryBroyden{SciMLBase._unwrap_val(threshold)}()
+end
+
+@views function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleLimitedMemoryBroyden,
+        args...; abstol = nothing, reltol = nothing, maxiters = 1000,
+        termination_condition = nothing, kwargs...)
+    @bb x = copy(float(prob.u0))
+    threshold = __get_threshold(alg)
+    η = min(SciMLBase._unwrap_val(threshold), maxiters)
+
+    # For scalar problems / if the threshold is larger than problem size just use Broyden
+    if x isa Number || length(x) ≤ η
+        return SciMLBase.__solve(prob, SimpleBroyden(), args...;
+            abstol, reltol, maxiters, termination_condition, kwargs...)
+    end
+
+    fx = _get_fx(prob, x)
+
+    U, Vᵀ = __init_low_rank_jacobian(x, fx, threshold)
+
+    abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fx, x,
+        termination_condition)
+
+    @bb xo = copy(x)
+    @bb δx = copy(fx)
+    @bb δx .*= -1
+    @bb fo = copy(fx)
+    @bb δf = copy(fx)
+
+    @bb vᵀ_cache = copy(x)
+    Tcache = __lbroyden_threshold_cache(x, threshold)
+    @bb mat_cache = copy(x)
+
+    for i in 1:maxiters
+        @bb @. x = xo + δx
+        fx = __eval_f(prob, fx, x)
+        @bb @. δf = fx - fo
+
+        # Termination Checks
+        tc_sol = check_termination(tc_cache, fx, x, xo, prob, alg)
+        tc_sol !== nothing && return tc_sol
+
+        _U = selectdim(U, 2, 1:min(η, i - 1))
+        _Vᵀ = selectdim(Vᵀ, 1, 1:min(η, i - 1))
+
+        vᵀ = _rmatvec!!(vᵀ_cache, Tcache, _U, _Vᵀ, δx)
+        mvec = _matvec!!(mat_cache, Tcache, _U, _Vᵀ, δf)
+        d = dot(vᵀ, δf)
+        @bb @. δx = (δx - mvec) / d
+
+        selectdim(U, 2, mod1(i, η)) .= δx
+        selectdim(Vᵀ, 1, mod1(i, η)) .= vᵀ
+
+        _U = selectdim(U, 2, 1:min(η, i))
+        _Vᵀ = selectdim(Vᵀ, 1, 1:min(η, i))
+        δx = _matvec!!(δx, Tcache, _U, _Vᵀ, fx)
+        @bb @. δx *= -1
+
+        @bb copyto!(xo, x)
+        @bb copyto!(fo, fx)
+    end
+
+    return build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)
+end
+
+function _rmatvec!!(y, xᵀU, U, Vᵀ, x)
+    # xᵀ × (-I + UVᵀ)
+    η = size(U, 2)
+    if η == 0
+        @bb @. y = -x
+        return y
+    end
+    x_ = vec(x)
+    xᵀU_ = view(xᵀU, 1:η)
+    @bb xᵀU_ = transpose(U) × x_
+    @bb y = transpose(Vᵀ) × xᵀU_
+    @bb @. y -= x
+    return y
+end
+
+function _matvec!!(y, Vᵀx, U, Vᵀ, x)
+    # (-I + UVᵀ) × x
+    η = size(U, 2)
+    if η == 0
+        @bb @. y = -x
+        return y
+    end
+    x_ = vec(x)
+    Vᵀx_ = view(Vᵀx, 1:η)
+    @bb Vᵀx_ = Vᵀ × x_
+    @bb y = U × Vᵀx_
+    @bb @. y -= x
+    return y
+end
+
+__lbroyden_threshold_cache(x, ::Val{threshold}) where {threshold} = similar(x, threshold)
+function __lbroyden_threshold_cache(x::SArray, ::Val{threshold}) where {threshold}
+    return SArray{Tuple{threshold}, eltype(x)}(ntuple(_ -> zero(eltype(x)), threshold))
+end
diff --git a/src/nlsolve/raphson.jl b/src/nlsolve/raphson.jl
new file mode 100644
index 0000000..22f7fba
--- /dev/null
+++ b/src/nlsolve/raphson.jl
@@ -0,0 +1,55 @@
+"""
+    SimpleNewtonRaphson(autodiff)
+    SimpleNewtonRaphson(; autodiff = AutoForwardDiff())
+
+A low-overhead implementation of Newton-Raphson. This method is non-allocating on scalar
+and static array problems.
+
+!!! note
+
+    As part of the decreased overhead, this method omits some of the higher level error
+    catching of the other methods. Thus, to see better error messages, use one of the other
+    methods like `NewtonRaphson`.
+
+### Keyword Arguments
+
+  - `autodiff`: determines the backend used for the Jacobian. Defaults to
+    `AutoForwardDiff()`. Valid choices are `AutoForwardDiff()` or `AutoFiniteDiff()`.
+"""
+@kwdef @concrete struct SimpleNewtonRaphson <: AbstractNewtonAlgorithm
+    autodiff = AutoForwardDiff()
+end
+
+const SimpleGaussNewton = SimpleNewtonRaphson
+
+function SciMLBase.__solve(prob::Union{NonlinearProblem, NonlinearLeastSquaresProblem},
+        alg::SimpleNewtonRaphson, args...; abstol = nothing, reltol = nothing,
+        maxiters = 1000, termination_condition = nothing, kwargs...)
+    @bb x = copy(float(prob.u0))
+    fx = _get_fx(prob, x)
+    @bb xo = copy(x)
+    J, jac_cache = jacobian_cache(alg.autodiff, prob.f, fx, x, prob.p)
+
+    abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fx, x,
+        termination_condition)
+
+    for i in 1:maxiters
+        fx, dfx = value_and_jacobian(alg.autodiff, prob.f, fx, x, prob.p, jac_cache; J)
+
+        if i == 1
+            if iszero(fx)
+                return build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
+            end
+        else
+            # Termination Checks
+            tc_sol = check_termination(tc_cache, fx, x, xo, prob, alg)
+            tc_sol !== nothing && return tc_sol
+        end
+
+        @bb copyto!(xo, x)
+        δx = _restructure(x, dfx \ _vec(fx))
+        @bb x .-= δx
+    end
+
+    return build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)
+end
diff --git a/src/nlsolve/trustRegion.jl b/src/nlsolve/trustRegion.jl
new file mode 100644
index 0000000..b4db396
--- /dev/null
+++ b/src/nlsolve/trustRegion.jl
@@ -0,0 +1,161 @@
+"""
+    SimpleTrustRegion(; autodiff = AutoForwardDiff(), max_trust_radius::Real = 0.0,
+        initial_trust_radius::Real = 0.0, step_threshold::Real = 0.1,
+        shrink_threshold::Real = 0.25, expand_threshold::Real = 0.75,
+        shrink_factor::Real = 0.25, expand_factor::Real = 2.0, max_shrink_times::Int = 32)
+
+A low-overhead implementation of a trust-region solver. This method is non-allocating on
+scalar and static array problems.
+
+### Keyword Arguments
+
+  - `autodiff`: determines the backend used for the Jacobian. Defaults to
+    `AutoForwardDiff()`. Valid choices are `AutoForwardDiff()` or `AutoFiniteDiff()`.
+  - `max_trust_radius`: the maximum radius of the trust region. Defaults to
+    `max(norm(f(u0)), maximum(u0) - minimum(u0))`.
+  - `initial_trust_radius`: the initial trust region radius. Defaults to
+    `max_trust_radius / 11`.
+  - `step_threshold`: the threshold for taking a step. In every iteration, the threshold is
+    compared with a value `r`, which is the actual reduction in the objective function divided
+    by the predicted reduction. If `step_threshold > r` the model is not a good approximation,
+    and the step is rejected. Defaults to `0.1`. For more details, see
+    [Rahpeymaii, F.](https://link.springer.com/article/10.1007/s40096-020-00339-4)
+  - `shrink_threshold`: the threshold for shrinking the trust region radius. In every
+    iteration, the threshold is compared with a value `r` which is the actual reduction in the
+    objective function divided by the predicted reduction. If `shrink_threshold > r` the trust
+    region radius is shrunk by `shrink_factor`. Defaults to `0.25`. For more details, see
+    [Rahpeymaii, F.](https://link.springer.com/article/10.1007/s40096-020-00339-4)
+  - `expand_threshold`: the threshold for expanding the trust region radius. If a step is
+    taken, i.e `step_threshold < r` (with `r` defined in `shrink_threshold`), a check is also
+    made to see if `expand_threshold < r`. If that is true, the trust region radius is
+    expanded by `expand_factor`. Defaults to `0.75`.
+  - `shrink_factor`: the factor to shrink the trust region radius with if
+    `shrink_threshold > r` (with `r` defined in `shrink_threshold`). Defaults to `0.25`.
+  - `expand_factor`: the factor to expand the trust region radius with if
+    `expand_threshold < r` (with `r` defined in `shrink_threshold`). Defaults to `2.0`.
+  - `max_shrink_times`: the maximum number of times to shrink the trust region radius in a
+    row, `max_shrink_times` is exceeded, the algorithm returns. Defaults to `32`.
+"""
+@kwdef @concrete struct SimpleTrustRegion <: AbstractNewtonAlgorithm
+    autodiff = AutoForwardDiff()
+    max_trust_radius = 0.0
+    initial_trust_radius = 0.0
+    step_threshold = 0.0001
+    shrink_threshold = 0.25
+    expand_threshold = 0.75
+    shrink_factor = 0.25
+    expand_factor = 2.0
+    max_shrink_times::Int = 32
+end
+
+function SciMLBase.__solve(prob::NonlinearProblem, alg::SimpleTrustRegion, args...;
+        abstol = nothing, reltol = nothing, maxiters = 1000,
+        termination_condition = nothing, kwargs...)
+    @bb x = copy(float(prob.u0))
+    T = eltype(real(x))
+    Δₘₐₓ = T(alg.max_trust_radius)
+    Δ = T(alg.initial_trust_radius)
+    η₁ = T(alg.step_threshold)
+    η₂ = T(alg.shrink_threshold)
+    η₃ = T(alg.expand_threshold)
+    t₁ = T(alg.shrink_factor)
+    t₂ = T(alg.expand_factor)
+    max_shrink_times = alg.max_shrink_times
+
+    fx = _get_fx(prob, x)
+    @bb xo = copy(x)
+    J, jac_cache = jacobian_cache(alg.autodiff, prob.f, fx, x, prob.p)
+    fx, ∇f = value_and_jacobian(alg.autodiff, prob.f, fx, x, prob.p, jac_cache; J)
+
+    abstol, reltol, tc_cache = init_termination_cache(abstol, reltol, fx, x,
+        termination_condition)
+
+    # Set default trust region radius if not specified by user.
+    Δₘₐₓ == 0 && (Δₘₐₓ = max(norm(fx), maximum(x) - minimum(x)))
+    Δ == 0 && (Δ = Δₘₐₓ / 11)
+
+    fₖ = 0.5 * norm(fx)^2
+    H = ∇f' * ∇f
+    g = _restructure(x, ∇f' * _vec(fx))
+    shrink_counter = 0
+
+    @bb δsd = copy(x)
+    @bb δN_δsd = copy(x)
+    @bb δN = copy(x)
+    @bb Hδ = copy(x)
+    dogleg_cache = (; δsd, δN_δsd, δN)
+
+    F = fx
+    for k in 1:maxiters
+        # Solve the trust region subproblem.
+        δ = dogleg_method!!(dogleg_cache, ∇f, fx, g, Δ)
+        @bb @. x = xo + δ
+
+        fx = __eval_f(prob, fx, x)
+
+        fₖ₊₁ = norm(fx)^2 / T(2)
+
+        # Compute the ratio of the actual to predicted reduction.
+        @bb Hδ = H × vec(δ)
+        r = (fₖ₊₁ - fₖ) / (dot(δ', g) + dot(δ', Hδ) / T(2))
+
+        # Update the trust region radius.
+        if r < η₂
+            Δ = t₁ * Δ
+            shrink_counter += 1
+            shrink_counter > max_shrink_times && return build_solution(prob, alg, x, fx;
+                retcode = ReturnCode.ConvergenceFailure)
+        else
+            shrink_counter = 0
+        end
+
+        if r > η₁
+            # Termination Checks
+            tc_sol = check_termination(tc_cache, fx, x, xo, prob, alg)
+            tc_sol !== nothing && return tc_sol
+
+            # Take the step.
+            @bb @. xo = x
+
+            fx, ∇f = value_and_jacobian(alg.autodiff, prob.f, fx, x, prob.p, jac_cache; J)
+
+            # Update the trust region radius.
+            (r > η₃) && (norm(δ) ≈ Δ) && (Δ = min(t₂ * Δ, Δₘₐₓ))
+            fₖ = fₖ₊₁
+
+            @bb H = transpose(∇f) × ∇f
+            @bb g = transpose(∇f) × vec(fx)
+        end
+    end
+
+    return build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)
+end
+
+function dogleg_method!!(cache, J, f, g, Δ)
+    (; δsd, δN_δsd, δN) = cache
+
+    # Compute the Newton step.
+    @bb δN .= _restructure(δN, J \ _vec(f))
+    @bb δN .*= -1
+    # Test if the full step is within the trust region.
+    (norm(δN) ≤ Δ) && return δN
+
+    # Calcualte Cauchy point, optimum along the steepest descent direction.
+    @bb δsd .= g
+    @bb @. δsd *= -1
+    norm_δsd = norm(δsd)
+    if (norm_δsd ≥ Δ)
+        @bb @. δsd *= Δ / norm_δsd
+        return δsd
+    end
+
+    # Find the intersection point on the boundary.
+    @bb @. δN_δsd = δN - δsd
+    dot_δN_δsd = dot(δN_δsd, δN_δsd)
+    dot_δsd_δN_δsd = dot(δsd, δN_δsd)
+    dot_δsd = dot(δsd, δsd)
+    fact = dot_δsd_δN_δsd^2 - dot_δN_δsd * (dot_δsd - Δ^2)
+    tau = (-dot_δsd_δN_δsd + sqrt(fact)) / dot_δN_δsd
+    @bb @. δsd += tau * δN_δsd
+    return δsd
+end
diff --git a/src/raphson.jl b/src/raphson.jl
deleted file mode 100644
index 138e672..0000000
--- a/src/raphson.jl
+++ /dev/null
@@ -1,125 +0,0 @@
-"""
-    SimpleNewtonRaphson(; batched = false,
-        chunk_size = Val{0}(),
-        autodiff = Val{true}(),
-        diff_type = Val{:forward},
-        termination_condition = missing)
-
-A low-overhead implementation of Newton-Raphson. This method is non-allocating on scalar
-and static array problems.
-
-!!! note
-
-    As part of the decreased overhead, this method omits some of the higher level error
-    catching of the other methods. Thus, to see better error messages, use one of the other
-    methods like `NewtonRaphson`
-
-### Keyword Arguments
-
-- `chunk_size`: the chunk size used by the internal ForwardDiff.jl automatic differentiation
-  system. This allows for multiple derivative columns to be computed simultaneously,
-  improving performance. Defaults to `0`, which is equivalent to using ForwardDiff.jl's
-  default chunk size mechanism. For more details, see the documentation for
-  [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/).
-- `autodiff`: whether to use forward-mode automatic differentiation for the Jacobian.
-  Note that this argument is ignored if an analytical Jacobian is passed; as that will be
-  used instead. Defaults to `Val{true}`, which means ForwardDiff.jl is used by default.
-  If `Val{false}`, then FiniteDiff.jl is used for finite differencing.
-- `diff_type`: the type of finite differencing used if `autodiff = false`. Defaults to
-  `Val{:forward}` for forward finite differences. For more details on the choices, see the
-  [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl) documentation.
-- `termination_condition`: control the termination of the algorithm. (Only works for batched
-  problems)
-"""
-struct SimpleNewtonRaphson{CS, AD, FDT} <: AbstractNewtonAlgorithm{CS, AD, FDT} end
-
-function SimpleNewtonRaphson(; batched = false,
-        chunk_size = Val{0}(),
-        autodiff = Val{true}(),
-        diff_type = Val{:forward},
-        termination_condition = missing)
-    if !ismissing(termination_condition) && !batched
-        throw(ArgumentError("`termination_condition` is currently only supported for batched problems"))
-    end
-    if batched
-        # @assert ADLinearSolveFDExtLoaded[] "Please install and load `LinearSolve.jl`, `FiniteDifferences.jl` and `AbstractDifferentiation.jl` to use batched Newton-Raphson."
-        termination_condition = ismissing(termination_condition) ?
-                                NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
-            abstol = nothing,
-            reltol = nothing) :
-                                termination_condition
-        return BatchedSimpleNewtonRaphson(; chunk_size,
-            autodiff,
-            diff_type,
-            termination_condition)
-        return SimpleNewtonRaphson{SciMLBase._unwrap_val(chunk_size),
-            SciMLBase._unwrap_val(autodiff),
-            SciMLBase._unwrap_val(diff_type)}()
-    end
-    return SimpleNewtonRaphson{SciMLBase._unwrap_val(chunk_size),
-        SciMLBase._unwrap_val(autodiff),
-        SciMLBase._unwrap_val(diff_type)}()
-end
-
-const SimpleGaussNewton = SimpleNewtonRaphson
-
-function SciMLBase.__solve(prob::Union{NonlinearProblem, NonlinearLeastSquaresProblem},
-        alg::SimpleNewtonRaphson, args...; abstol = nothing,
-        reltol = nothing,
-        maxiters = 1000, kwargs...)
-    f = Base.Fix2(prob.f, prob.p)
-    x = float(prob.u0)
-    fx = float(prob.u0)
-    T = typeof(x)
-
-    if SciMLBase.isinplace(prob)
-        error("SimpleNewtonRaphson currently only supports out-of-place nonlinear problems")
-    end
-
-    if prob isa NonlinearLeastSquaresProblem &&
-       !(typeof(prob.u0) <: Union{Number, AbstractVector})
-        error("SimpleGaussNewton only supports Number and AbstactVector types. Please convert any problem of AbstractArray into one with u0 as AbstractVector")
-    end
-
-    atol = abstol !== nothing ? abstol :
-           real(oneunit(eltype(T))) * (eps(real(one(eltype(T)))))^(4 // 5)
-    rtol = reltol !== nothing ? reltol : eps(real(one(eltype(T))))^(4 // 5)
-
-    if x isa Number
-        xo = oftype(one(eltype(x)), Inf)
-    else
-        xo = map(x -> oftype(one(eltype(x)), Inf), x)
-    end
-
-    for i in 1:maxiters
-        if DiffEqBase.has_jac(prob.f)
-            dfx = prob.f.jac(x, prob.p)
-            fx = f(x)
-        elseif alg_autodiff(alg)
-            fx, dfx = value_derivative(f, x)
-        elseif x isa AbstractArray
-            fx = f(x)
-            dfx = FiniteDiff.finite_difference_jacobian(f, x, diff_type(alg), eltype(x), fx)
-        else
-            fx = f(x)
-            dfx = FiniteDiff.finite_difference_derivative(f, x, diff_type(alg), eltype(x),
-                fx)
-        end
-        iszero(fx) &&
-            return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
-
-        if prob isa NonlinearProblem
-            Δx = _restructure(fx, dfx \ _vec(fx))
-        else
-            Δx = dfx \ fx
-        end
-
-        x -= Δx
-        if isapprox(x, xo, atol = atol, rtol = rtol)
-            return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
-        end
-        xo = x
-    end
-
-    return SciMLBase.build_solution(prob, alg, x, fx; retcode = ReturnCode.MaxIters)
-end
diff --git a/src/ridder.jl b/src/ridder.jl
deleted file mode 100644
index eabd7b2..0000000
--- a/src/ridder.jl
+++ /dev/null
@@ -1,95 +0,0 @@
-"""
-`Ridder()`
-
-A non-allocating ridder method
-
-"""
-struct Ridder <: AbstractBracketingAlgorithm end
-
-function SciMLBase.solve(prob::IntervalNonlinearProblem, alg::Ridder, args...;
-        maxiters = 1000, abstol = min(eps(prob.tspan[1]), eps(prob.tspan[2])),
-        kwargs...)
-    f = Base.Fix2(prob.f, prob.p)
-    left, right = prob.tspan
-    fl, fr = f(left), f(right)
-
-    if iszero(fl)
-        return SciMLBase.build_solution(prob, alg, left, fl;
-            retcode = ReturnCode.ExactSolutionLeft, left = left,
-            right = right)
-    elseif iszero(fr)
-        return SciMLBase.build_solution(prob, alg, right, fr;
-            retcode = ReturnCode.ExactSolutionRight, left = left,
-            right = right)
-    end
-
-    xo = oftype(left, Inf)
-    i = 1
-    if !iszero(fr)
-        while i < maxiters
-            mid = (left + right) / 2
-            (mid == left || mid == right) &&
-                return SciMLBase.build_solution(prob, alg, left, fl;
-                    retcode = ReturnCode.FloatingPointLimit,
-                    left = left, right = right)
-            fm = f(mid)
-            s = sqrt(fm^2 - fl * fr)
-            iszero(s) &&
-                return SciMLBase.build_solution(prob, alg, left, fl;
-                    retcode = ReturnCode.Failure,
-                    left = left, right = right)
-            x = mid + (mid - left) * sign(fl - fr) * fm / s
-            fx = f(x)
-            xo = x
-            if abs((right - left) / 2) < abstol
-                return SciMLBase.build_solution(prob, alg, mid, fm;
-                    retcode = ReturnCode.Success,
-                    left = left, right = right)
-            end
-            if iszero(fx)
-                right = x
-                fr = fx
-                break
-            end
-            if sign(fx) != sign(fm)
-                left = mid
-                fl = fm
-                right = x
-                fr = fx
-            elseif sign(fx) != sign(fl)
-                right = x
-                fr = fx
-            else
-                @assert sign(fx) != sign(fr)
-                left = x
-                fl = fx
-            end
-            i += 1
-        end
-    end
-
-    while i < maxiters
-        mid = (left + right) / 2
-        (mid == left || mid == right) &&
-            return SciMLBase.build_solution(prob, alg, left, fl;
-                retcode = ReturnCode.FloatingPointLimit,
-                left = left, right = right)
-        fm = f(mid)
-        if abs((right - left) / 2) < abstol
-            return SciMLBase.build_solution(prob, alg, mid, fm;
-                retcode = ReturnCode.Success,
-                left = left, right = right)
-        end
-        if iszero(fm)
-            right = mid
-            fr = fm
-        else
-            left = mid
-            fl = fm
-        end
-        i += 1
-    end
-
-    return SciMLBase.build_solution(prob, alg, left, fl; retcode = ReturnCode.MaxIters,
-        left = left, right = right)
-end
diff --git a/src/trustRegion.jl b/src/trustRegion.jl
deleted file mode 100644
index 0fba7b1..0000000
--- a/src/trustRegion.jl
+++ /dev/null
@@ -1,196 +0,0 @@
-"""
-```julia
-SimpleTrustRegion(; chunk_size = Val{0}(),
-                    autodiff = Val{true}(),
-                    diff_type = Val{:forward},
-                    max_trust_radius::Real = 0.0,
-                    initial_trust_radius::Real = 0.0,
-                    step_threshold::Real = 0.1,
-                    shrink_threshold::Real = 0.25,
-                    expand_threshold::Real = 0.75,
-                    shrink_factor::Real = 0.25,
-                    expand_factor::Real = 2.0,
-                    max_shrink_times::Int = 32
-```
-
-A low-overhead implementation of a trust-region solver.
-
-### Keyword Arguments
-
-- `chunk_size`: the chunk size used by the internal ForwardDiff.jl automatic differentiation
-  system. This allows for multiple derivative columns to be computed simultaneously,
-  improving performance. Defaults to `0`, which is equivalent to using ForwardDiff.jl's
-  default chunk size mechanism. For more details, see the documentation for
-  [ForwardDiff.jl](https://juliadiff.org/ForwardDiff.jl/stable/).
-- `autodiff`: whether to use forward-mode automatic differentiation for the Jacobian.
-  Note that this argument is ignored if an analytical Jacobian is passed; as that will be
-  used instead. Defaults to `Val{true}`, which means ForwardDiff.jl is used by default.
-  If `Val{false}`, then FiniteDiff.jl is used for finite differencing.
-- `diff_type`: the type of finite differencing used if `autodiff = false`. Defaults to
-  `Val{:forward}` for forward finite differences. For more details on the choices, see the
-  [FiniteDiff.jl](https://github.com/JuliaDiff/FiniteDiff.jl) documentation.
-- `max_trust_radius`: the maximum radius of the trust region. Defaults to
-  `max(norm(f(u0)), maximum(u0) - minimum(u0))`.
-- `initial_trust_radius`: the initial trust region radius. Defaults to
-  `max_trust_radius / 11`.
-- `step_threshold`: the threshold for taking a step. In every iteration, the threshold is
-  compared with a value `r`, which is the actual reduction in the objective function divided
-  by the predicted reduction. If `step_threshold > r` the model is not a good approximation,
-  and the step is rejected. Defaults to `0.1`. For more details, see
-  [Rahpeymaii, F.](https://link.springer.com/article/10.1007/s40096-020-00339-4)
-- `shrink_threshold`: the threshold for shrinking the trust region radius. In every
-  iteration, the threshold is compared with a value `r` which is the actual reduction in the
-  objective function divided by the predicted reduction. If `shrink_threshold > r` the trust
-  region radius is shrunk by `shrink_factor`. Defaults to `0.25`. For more details, see
-  [Rahpeymaii, F.](https://link.springer.com/article/10.1007/s40096-020-00339-4)
-- `expand_threshold`: the threshold for expanding the trust region radius. If a step is
-  taken, i.e `step_threshold < r` (with `r` defined in `shrink_threshold`), a check is also
-  made to see if `expand_threshold < r`. If that is true, the trust region radius is
-  expanded by `expand_factor`. Defaults to `0.75`.
-- `shrink_factor`: the factor to shrink the trust region radius with if
-  `shrink_threshold > r` (with `r` defined in `shrink_threshold`). Defaults to `0.25`.
-- `expand_factor`: the factor to expand the trust region radius with if
-  `expand_threshold < r` (with `r` defined in `shrink_threshold`). Defaults to `2.0`.
-- `max_shrink_times`: the maximum number of times to shrink the trust region radius in a
-  row, `max_shrink_times` is exceeded, the algorithm returns. Defaults to `32`.
-"""
-struct SimpleTrustRegion{T, CS, AD, FDT} <: AbstractNewtonAlgorithm{CS, AD, FDT}
-    max_trust_radius::T
-    initial_trust_radius::T
-    step_threshold::T
-    shrink_threshold::T
-    expand_threshold::T
-    shrink_factor::T
-    expand_factor::T
-    max_shrink_times::Int
-    function SimpleTrustRegion(; chunk_size = Val{0}(),
-            autodiff = Val{true}(),
-            diff_type = Val{:forward},
-            max_trust_radius::Real = 0.0,
-            initial_trust_radius::Real = 0.0,
-            step_threshold::Real = 0.0001,
-            shrink_threshold::Real = 0.25,
-            expand_threshold::Real = 0.75,
-            shrink_factor::Real = 0.25,
-            expand_factor::Real = 2.0,
-            max_shrink_times::Int = 32)
-        new{typeof(initial_trust_radius),
-            SciMLBase._unwrap_val(chunk_size),
-            SciMLBase._unwrap_val(autodiff),
-            SciMLBase._unwrap_val(diff_type)}(max_trust_radius,
-            initial_trust_radius,
-            step_threshold,
-            shrink_threshold,
-            expand_threshold,
-            shrink_factor,
-            expand_factor,
-            max_shrink_times)
-    end
-end
-
-function SciMLBase.__solve(prob::NonlinearProblem,
-        alg::SimpleTrustRegion, args...; abstol = nothing,
-        reltol = nothing,
-        maxiters = 1000, kwargs...)
-    f = Base.Fix2(prob.f, prob.p)
-    x = float(prob.u0)
-    T = typeof(x)
-    Δₘₐₓ = float(alg.max_trust_radius)
-    Δ = float(alg.initial_trust_radius)
-    η₁ = float(alg.step_threshold)
-    η₂ = float(alg.shrink_threshold)
-    η₃ = float(alg.expand_threshold)
-    t₁ = float(alg.shrink_factor)
-    t₂ = float(alg.expand_factor)
-    max_shrink_times = alg.max_shrink_times
-
-    if SciMLBase.isinplace(prob)
-        error("SimpleTrustRegion currently only supports out-of-place nonlinear problems")
-    end
-
-    atol = abstol !== nothing ? abstol :
-           real(oneunit(eltype(T))) * (eps(real(one(eltype(T)))))^(4 // 5)
-    rtol = reltol !== nothing ? reltol : eps(real(one(eltype(T))))^(4 // 5)
-
-    if DiffEqBase.has_jac(prob.f)
-        ∇f = prob.f.jac(x, prob.p)
-        F = f(x)
-    elseif alg_autodiff(alg)
-        F, ∇f = value_derivative(f, x)
-    elseif x isa AbstractArray
-        F = f(x)
-        ∇f = FiniteDiff.finite_difference_jacobian(f, x, diff_type(alg), eltype(x), F)
-    else
-        F = f(x)
-        ∇f = FiniteDiff.finite_difference_derivative(f, x, diff_type(alg), eltype(x), F)
-    end
-
-    # Set default trust region radius if not specified by user.
-    if Δₘₐₓ == 0.0
-        Δₘₐₓ = max(norm(F), maximum(x) - minimum(x))
-    end
-    if Δ == 0.0
-        Δ = Δₘₐₓ / 11
-    end
-
-    fₖ = 0.5 * norm(F)^2
-    H = ∇f' * ∇f
-    g = ∇f' * F
-    shrink_counter = 0
-
-    for k in 1:maxiters
-        # Solve the trust region subproblem.
-        δ = dogleg_method(∇f, F, g, Δ)
-        xₖ₊₁ = x + δ
-        Fₖ₊₁ = f(xₖ₊₁)
-        fₖ₊₁ = 0.5 * norm(Fₖ₊₁)^2
-
-        # Compute the ratio of the actual to predicted reduction.
-        model = -(δ' * g + 0.5 * δ' * H * δ)
-        r = (fₖ - fₖ₊₁) / model
-
-        # Update the trust region radius.
-        if r < η₂
-            Δ = t₁ * Δ
-            shrink_counter += 1
-            if shrink_counter > max_shrink_times
-                return SciMLBase.build_solution(prob, alg, x, F;
-                    retcode = ReturnCode.Success)
-            end
-        else
-            shrink_counter = 0
-        end
-        if r > η₁
-            if isapprox(xₖ₊₁, x, atol = atol, rtol = rtol)
-                return SciMLBase.build_solution(prob, alg, xₖ₊₁, Fₖ₊₁;
-                    retcode = ReturnCode.Success)
-            end
-            # Take the step.
-            x = xₖ₊₁
-            F = Fₖ₊₁
-            if alg_autodiff(alg)
-                F, ∇f = value_derivative(f, x)
-            elseif x isa AbstractArray
-                ∇f = FiniteDiff.finite_difference_jacobian(f, x, diff_type(alg), eltype(x),
-                    F)
-            else
-                ∇f = FiniteDiff.finite_difference_derivative(f, x, diff_type(alg),
-                    eltype(x),
-                    F)
-            end
-
-            iszero(F) &&
-                return SciMLBase.build_solution(prob, alg, x, F;
-                    retcode = ReturnCode.Success)
-
-            # Update the trust region radius.
-            if r > η₃ && norm(δ) ≈ Δ
-                Δ = min(t₂ * Δ, Δₘₐₓ)
-            end
-            fₖ = fₖ₊₁
-            H = ∇f' * ∇f
-            g = ∇f' * F
-        end
-    end
-    return SciMLBase.build_solution(prob, alg, x, F; retcode = ReturnCode.MaxIters)
-end
diff --git a/src/utils.jl b/src/utils.jl
index af66f63..444128b 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,86 +1,238 @@
+struct SimpleNonlinearSolveTag end
+
+function ForwardDiff.checktag(::Type{<:ForwardDiff.Tag{<:SimpleNonlinearSolveTag, <:T}},
+        f::F, x::AbstractArray{T}) where {T, F}
+    return true
+end
+
 """
-  prevfloat_tdir(x, x0, x1)
+    __prevfloat_tdir(x, x0, x1)
 
 Move `x` one floating point towards x0.
 """
-function prevfloat_tdir(x, x0, x1)
-    x1 > x0 ? prevfloat(x) : nextfloat(x)
-end
+__prevfloat_tdir(x, x0, x1) = ifelse(x1 > x0, prevfloat(x), nextfloat(x))
+
+"""
+    __nextfloat_tdir(x, x0, x1)
+
+Move `x` one floating point towards x1.
+"""
+__nextfloat_tdir(x, x0, x1) = ifelse(x1 > x0, nextfloat(x), prevfloat(x))
 
-function nextfloat_tdir(x, x0, x1)
-    x1 > x0 ? nextfloat(x) : prevfloat(x)
+"""
+    __max_tdir(a, b, x0, x1)
+
+Return the maximum of `a` and `b` if `x1 > x0`, otherwise return the minimum.
+"""
+__max_tdir(a, b, x0, x1) = ifelse(x1 > x0, max(a, b), min(a, b))
+
+__cvt_real(::Type{T}, ::Nothing) where {T} = nothing
+__cvt_real(::Type{T}, x) where {T} = real(T(x))
+
+_get_tolerance(η, ::Type{T}) where {T} = __cvt_real(T, η)
+function _get_tolerance(::Nothing, ::Type{T}) where {T}
+    η = real(oneunit(T)) * (eps(real(one(T))))^(4 // 5)
+    return _get_tolerance(η, T)
 end
 
-function max_tdir(a, b, x0, x1)
-    x1 > x0 ? max(a, b) : min(a, b)
+__standard_tag(::Nothing, x) = ForwardDiff.Tag(SimpleNonlinearSolveTag(), eltype(x))
+__standard_tag(tag::ForwardDiff.Tag, _) = tag
+__standard_tag(tag, x) = ForwardDiff.Tag(tag, eltype(x))
+
+__pick_forwarddiff_chunk(x) = ForwardDiff.Chunk(length(x))
+function __pick_forwarddiff_chunk(x::StaticArray)
+    L = prod(Size(x))
+    if L ≤ ForwardDiff.DEFAULT_CHUNK_THRESHOLD
+        return ForwardDiff.Chunk{L}()
+    else
+        return ForwardDiff.Chunk{ForwardDiff.DEFAULT_CHUNK_THRESHOLD}()
+    end
 end
 
-alg_autodiff(alg::AbstractNewtonAlgorithm{CS, AD, FDT}) where {CS, AD, FDT} = AD
-diff_type(alg::AbstractNewtonAlgorithm{CS, AD, FDT}) where {CS, AD, FDT} = FDT
+function __get_jacobian_config(ad::AutoForwardDiff{CS}, f, x) where {CS}
+    ck = (CS === nothing || CS ≤ 0) ? __pick_forwarddiff_chunk(x) : ForwardDiff.Chunk{CS}()
+    tag = __standard_tag(ad.tag, x)
+    return ForwardDiff.JacobianConfig(f, x, ck, tag)
+end
+function __get_jacobian_config(ad::AutoForwardDiff{CS}, f!, y, x) where {CS}
+    ck = (CS === nothing || CS ≤ 0) ? __pick_forwarddiff_chunk(x) : ForwardDiff.Chunk{CS}()
+    tag = __standard_tag(ad.tag, x)
+    return ForwardDiff.JacobianConfig(f!, y, x, ck, tag)
+end
 
 """
-  value_derivative(f, x)
+    value_and_jacobian(ad, f, y, x, p, cache; J = nothing)
 
-Compute `f(x), d/dx f(x)` in the most efficient way.
+Compute `f(x), d/dx f(x)` in the most efficient way based on `ad`. None of the arguments
+except `cache` (& `J` if not nothing) are mutated.
 """
-function value_derivative(f::F, x::R) where {F, R}
-    T = typeof(ForwardDiff.Tag(f, R))
-    out = f(ForwardDiff.Dual{T}(x, one(x)))
-    ForwardDiff.value(out), ForwardDiff.extract_derivative(T, out)
+function value_and_jacobian(ad, f::F, y, x::X, p, cache; J = nothing) where {F, X}
+    if isinplace(f)
+        _f = (du, u) -> f(du, u, p)
+        if DiffEqBase.has_jac(f)
+            f.jac(J, x, p)
+            _f(y, x)
+            return y, J
+        elseif ad isa AutoForwardDiff
+            res = DiffResults.DiffResult(y, J)
+            ForwardDiff.jacobian!(res, _f, y, x, cache)
+            return DiffResults.value(res), DiffResults.jacobian(res)
+        elseif ad isa AutoFiniteDiff
+            FiniteDiff.finite_difference_jacobian!(J, _f, x, cache)
+            _f(y, x)
+            return y, J
+        else
+            throw(ArgumentError("Unsupported AD method: $(ad)"))
+        end
+    else
+        _f = Base.Fix2(f, p)
+        if DiffEqBase.has_jac(f)
+            return _f(x), f.jac(x, p)
+        elseif ad isa AutoForwardDiff
+            if ArrayInterface.can_setindex(x)
+                res = DiffResults.DiffResult(y, J)
+                ForwardDiff.jacobian!(res, _f, x, cache)
+                return DiffResults.value(res), DiffResults.jacobian(res)
+            else
+                J_fd = ForwardDiff.jacobian(_f, x, cache)
+                return _f(x), J_fd
+            end
+        elseif ad isa AutoFiniteDiff
+            J_fd = FiniteDiff.finite_difference_jacobian(_f, x, cache)
+            return _f(x), J_fd
+        else
+            throw(ArgumentError("Unsupported AD method: $(ad)"))
+        end
+    end
+end
+
+function value_and_jacobian(ad, f::F, y, x::Number, p, cache; J = nothing) where {F}
+    if DiffEqBase.has_jac(f)
+        return f(x, p), f.jac(x, p)
+    elseif ad isa AutoForwardDiff
+        T = typeof(__standard_tag(ad.tag, x))
+        out = f(ForwardDiff.Dual{T}(x, one(x)), p)
+        return ForwardDiff.value(out), ForwardDiff.extract_derivative(T, out)
+    elseif ad isa AutoFiniteDiff
+        _f = Base.Fix2(f, p)
+        return _f(x), FiniteDiff.finite_difference_derivative(_f, x, ad.fdtype)
+    else
+        throw(ArgumentError("Unsupported AD method: $(ad)"))
+    end
 end
-value_derivative(f::F, x::AbstractArray) where {F} = f(x), ForwardDiff.jacobian(f, x)
 
 """
-    value_derivative!(J, y, f!, x, cfg = JacobianConfig(f!, y, x))
+    jacobian_cache(ad, f, y, x, p) --> J, cache
 
-Inplace version of [`SimpleNonlinearSolve.value_derivative`](@ref).
+Returns a Jacobian Matrix and a cache for the Jacobian computation.
 """
-function value_derivative!(J::AbstractMatrix,
-        y::AbstractArray,
-        f!::F,
-        x::AbstractArray,
-        cfg::ForwardDiff.JacobianConfig = ForwardDiff.JacobianConfig(f!, y, x)) where {F}
-    ForwardDiff.jacobian!(J, f!, y, x, cfg)
-    return y, J
-end
-
-value(x) = x
-value(x::Dual) = ForwardDiff.value(x)
-value(x::AbstractArray{<:Dual}) = map(ForwardDiff.value, x)
-
-function init_J(x)
-    J = ArrayInterface.zeromatrix(x)
-    if ismutable(x)
-        J[diagind(J)] .= one(eltype(x))
+function jacobian_cache(ad, f::F, y, x::X, p) where {F, X <: AbstractArray}
+    if isinplace(f)
+        _f = (du, u) -> f(du, u, p)
+        J = similar(y, length(y), length(x))
+        if DiffEqBase.has_jac(f)
+            return J, nothing
+        elseif ad isa AutoForwardDiff
+            return J, __get_jacobian_config(ad, _f, y, x)
+        elseif ad isa AutoFiniteDiff
+            return J, FiniteDiff.JacobianCache(copy(x), copy(y), copy(y), ad.fdtype)
+        else
+            throw(ArgumentError("Unsupported AD method: $(ad)"))
+        end
     else
-        J += I
+        _f = Base.Fix2(f, p)
+        if DiffEqBase.has_jac(f)
+            return nothing, nothing
+        elseif ad isa AutoForwardDiff
+            J = ArrayInterface.can_setindex(x) ? similar(y, length(y), length(x)) : nothing
+            return J, __get_jacobian_config(ad, _f, x)
+        elseif ad isa AutoFiniteDiff
+            return nothing, FiniteDiff.JacobianCache(copy(x), copy(y), copy(y), ad.fdtype)
+        else
+            throw(ArgumentError("Unsupported AD method: $(ad)"))
+        end
     end
-    return J
 end
 
-function dogleg_method(J, f, g, Δ)
-    # Compute the Newton step.
-    δN = J \ (-f)
-    # Test if the full step is within the trust region.
-    if norm(δN) ≤ Δ
-        return δN
+jacobian_cache(ad, f::F, y, x::Number, p) where {F} = nothing, nothing
+
+function compute_jacobian_and_hessian(ad::AutoForwardDiff, prob, _, x::Number)
+    fx = prob.f(x, prob.p)
+    J_fn = Base.Fix1(ForwardDiff.derivative, Base.Fix2(prob.f, prob.p))
+    dfx = J_fn(x)
+    d2fx = ForwardDiff.derivative(J_fn, x)
+    return fx, dfx, d2fx
+end
+
+function compute_jacobian_and_hessian(ad::AutoForwardDiff, prob, fx, x)
+    if isinplace(prob)
+        error("Inplace version for Nested ForwardDiff Not Implemented Yet!")
+    else
+        f = Base.Fix2(prob.f, prob.p)
+        fx = f(x)
+        J_fn = Base.Fix1(ForwardDiff.jacobian, f)
+        dfx = J_fn(x)
+        d2fx = ForwardDiff.jacobian(J_fn, x)
+        return fx, dfx, d2fx
     end
+end
 
-    # Calcualte Cauchy point, optimum along the steepest descent direction.
-    δsd = -g
-    norm_δsd = norm(δsd)
-    if norm_δsd ≥ Δ
-        return δsd .* Δ / norm_δsd
+function compute_jacobian_and_hessian(ad::AutoFiniteDiff, prob, _, x::Number)
+    fx = prob.f(x, prob.p)
+    J_fn = x -> FiniteDiff.finite_difference_derivative(Base.Fix2(prob.f, prob.p), x,
+        ad.fdtype)
+    dfx = J_fn(x)
+    d2fx = FiniteDiff.finite_difference_derivative(J_fn, x, ad.fdtype)
+    return fx, dfx, d2fx
+end
+
+function compute_jacobian_and_hessian(ad::AutoFiniteDiff, prob, fx, x)
+    if isinplace(prob)
+        error("Inplace version for Nested FiniteDiff Not Implemented Yet!")
+    else
+        f = Base.Fix2(prob.f, prob.p)
+        fx = f(x)
+        J_fn = x -> FiniteDiff.finite_difference_jacobian(f, x, ad.fdtype)
+        dfx = J_fn(x)
+        d2fx = FiniteDiff.finite_difference_jacobian(J_fn, x, ad.fdtype)
+        return fx, dfx, d2fx
     end
+end
 
-    # Find the intersection point on the boundary.
-    δN_δsd = δN - δsd
-    dot_δN_δsd = dot(δN_δsd, δN_δsd)
-    dot_δsd_δN_δsd = dot(δsd, δN_δsd)
-    dot_δsd = dot(δsd, δsd)
-    fact = dot_δsd_δN_δsd^2 - dot_δN_δsd * (dot_δsd - Δ^2)
-    tau = (-dot_δsd_δN_δsd + sqrt(fact)) / dot_δN_δsd
-    return δsd + tau * δN_δsd
+__init_identity_jacobian(u::Number, _) = one(u)
+__init_identity_jacobian!!(J::Number) = one(J)
+function __init_identity_jacobian(u, fu)
+    J = similar(u, promote_type(eltype(u), eltype(fu)), length(fu), length(u))
+    fill!(J, zero(eltype(J)))
+    J[diagind(J)] .= one(eltype(J))
+    return J
+end
+function __init_identity_jacobian!!(J)
+    fill!(J, zero(eltype(J)))
+    J[diagind(J)] .= one(eltype(J))
+    return J
+end
+function __init_identity_jacobian(u::StaticArray, fu)
+    S1, S2 = length(fu), length(u)
+    J = SMatrix{S1, S2, eltype(u)}(I)
+    return J
+end
+function __init_identity_jacobian!!(J::SMatrix{S1, S2}) where {S1, S2}
+    return SMatrix{S1, S2, eltype(J)}(I)
+end
+
+function __init_low_rank_jacobian(u::StaticArray{S1, T1}, fu::StaticArray{S2, T2},
+        ::Val{threshold}) where {S1, S2, T1, T2, threshold}
+    T = promote_type(T1, T2)
+    fuSize, uSize = Size(fu), Size(u)
+    Vᵀ = MArray{Tuple{threshold, prod(uSize)}, T}(undef)
+    U = MArray{Tuple{prod(fuSize), threshold}, T}(undef)
+    return U, Vᵀ
+end
+function __init_low_rank_jacobian(u, fu, ::Val{threshold}) where {threshold}
+    Vᵀ = similar(u, threshold, length(u))
+    U = similar(u, length(fu), threshold)
+    return U, Vᵀ
 end
 
 @inline _vec(v) = vec(v)
@@ -89,3 +241,97 @@ end
 
 @inline _restructure(y::Number, x::Number) = x
 @inline _restructure(y, x) = ArrayInterface.restructure(y, x)
+
+@inline function _get_fx(prob::NonlinearLeastSquaresProblem, x)
+    isinplace(prob) && prob.f.resid_prototype === nothing &&
+        error("Inplace NonlinearLeastSquaresProblem requires a `resid_prototype`")
+    return _get_fx(prob.f, x, prob.p)
+end
+@inline _get_fx(prob::NonlinearProblem, x) = _get_fx(prob.f, x, prob.p)
+@inline function _get_fx(f::NonlinearFunction, x, p)
+    if isinplace(f)
+        if f.resid_prototype !== nothing
+            T = eltype(x)
+            return T.(f.resid_prototype)
+        else
+            fx = similar(x)
+            f(fx, x, p)
+            return fx
+        end
+    else
+        return f(x, p)
+    end
+end
+
+# Termination Conditions Support
+# Taken directly from NonlinearSolve.jl
+# The default here is different from NonlinearSolve since the userbases are assumed to be
+# different. NonlinearSolve is more for robust / cached solvers while SimpleNonlinearSolve
+# is meant for low overhead solvers, users can opt into the other termination modes but the
+# default is to use the least overhead version.
+function init_termination_cache(abstol, reltol, du, u, ::Nothing)
+    return init_termination_cache(abstol, reltol, du, u, AbsNormTerminationMode())
+end
+function init_termination_cache(abstol, reltol, du, u, tc::AbstractNonlinearTerminationMode)
+    T = promote_type(eltype(du), eltype(u))
+    abstol !== nothing && (abstol = T(abstol))
+    reltol !== nothing && (reltol = T(reltol))
+    tc_cache = init(du, u, tc; abstol, reltol)
+    return DiffEqBase.get_abstol(tc_cache), DiffEqBase.get_reltol(tc_cache), tc_cache
+end
+
+function check_termination(tc_cache, fx, x, xo, prob, alg)
+    return check_termination(tc_cache, fx, x, xo, prob, alg,
+        DiffEqBase.get_termination_mode(tc_cache))
+end
+function check_termination(tc_cache, fx, x, xo, prob, alg,
+        ::AbstractNonlinearTerminationMode)
+    if tc_cache(fx, x, xo)
+        return build_solution(prob, alg, x, fx; retcode = ReturnCode.Success)
+    end
+    return nothing
+end
+function check_termination(tc_cache, fx, x, xo, prob, alg,
+        ::AbstractSafeNonlinearTerminationMode)
+    if tc_cache(fx, x, xo)
+        if tc_cache.retcode == NonlinearSafeTerminationReturnCode.Success
+            retcode = ReturnCode.Success
+        elseif tc_cache.retcode == NonlinearSafeTerminationReturnCode.PatienceTermination
+            retcode = ReturnCode.ConvergenceFailure
+        elseif tc_cache.retcode == NonlinearSafeTerminationReturnCode.ProtectiveTermination
+            retcode = ReturnCode.Unstable
+        else
+            error("Unknown termination code: $(tc_cache.retcode)")
+        end
+        return build_solution(prob, alg, x, fx; retcode)
+    end
+    return nothing
+end
+function check_termination(tc_cache, fx, x, xo, prob, alg,
+        ::AbstractSafeBestNonlinearTerminationMode)
+    if tc_cache(fx, x, xo)
+        if tc_cache.retcode == NonlinearSafeTerminationReturnCode.Success
+            retcode = ReturnCode.Success
+        elseif tc_cache.retcode == NonlinearSafeTerminationReturnCode.PatienceTermination
+            retcode = ReturnCode.ConvergenceFailure
+        elseif tc_cache.retcode == NonlinearSafeTerminationReturnCode.ProtectiveTermination
+            retcode = ReturnCode.Unstable
+        else
+            error("Unknown termination code: $(tc_cache.retcode)")
+        end
+        if isinplace(prob)
+            prob.f(fx, x, prob.p)
+        else
+            fx = prob.f(x, prob.p)
+        end
+        return build_solution(prob, alg, tc_cache.u, fx; retcode)
+    end
+    return nothing
+end
+
+@inline value(x) = x
+@inline value(x::Dual) = ForwardDiff.value(x)
+@inline value(x::AbstractArray{<:Dual}) = map(ForwardDiff.value, x)
+
+@inline __eval_f(prob, fx, x) = isinplace(prob) ? (prob.f(fx, x, prob.p); fx) :
+                                prob.f(x, prob.p)
diff --git a/test/23_test_problems.jl b/test/23_test_problems.jl
new file mode 100644
index 0000000..a0d5bc6
--- /dev/null
+++ b/test/23_test_problems.jl
@@ -0,0 +1,88 @@
+using SimpleNonlinearSolve, LinearAlgebra, NonlinearProblemLibrary, DiffEqBase, Test
+
+problems = NonlinearProblemLibrary.problems
+dicts = NonlinearProblemLibrary.dicts
+
+function test_on_library(problems, dicts, alg_ops, broken_tests, ϵ = 1e-4;
+        skip_tests = nothing)
+    for (idx, (problem, dict)) in enumerate(zip(problems, dicts))
+        x = dict["start"]
+        res = similar(x)
+        nlprob = NonlinearProblem(problem, copy(x))
+        @testset "$idx: $(dict["title"])" begin
+            for alg in alg_ops
+                try
+                    sol = solve(nlprob, alg;
+                        termination_condition = AbsNormTerminationMode())
+                    problem(res, sol.u, nothing)
+
+                    skip = skip_tests !== nothing && idx in skip_tests[alg]
+                    if skip
+                        @test_skip norm(res) ≤ ϵ
+                        continue
+                    end
+                    broken = idx in broken_tests[alg] ? true : false
+                    @test norm(res)≤ϵ broken=broken
+                catch e
+                    @error e
+                    broken = idx in broken_tests[alg] ? true : false
+                    if broken
+                        @test false broken=true
+                    else
+                        @test 1 == 2
+                    end
+                end
+            end
+        end
+    end
+end
+
+@testset "SimpleNewtonRaphson 23 Test Problems" begin
+    alg_ops = (SimpleNewtonRaphson(),)
+
+    # dictionary with indices of test problems where method does not converge to small residual
+    broken_tests = Dict(alg => Int[] for alg in alg_ops)
+    broken_tests[alg_ops[1]] = [6]
+
+    test_on_library(problems, dicts, alg_ops, broken_tests)
+end
+
+@testset "SimpleTrustRegion 23 Test Problems" begin
+    alg_ops = (SimpleTrustRegion(),)
+
+    # dictionary with indices of test problems where method does not converge to small residual
+    broken_tests = Dict(alg => Int[] for alg in alg_ops)
+    broken_tests[alg_ops[1]] = [3, 6, 15, 16, 21]
+
+    test_on_library(problems, dicts, alg_ops, broken_tests)
+end
+
+@testset "SimpleDFSane 23 Test Problems" begin
+    alg_ops = (SimpleDFSane(),)
+
+    broken_tests = Dict(alg => Int[] for alg in alg_ops)
+    broken_tests[alg_ops[1]] = [1, 2, 3, 4, 5, 6, 7, 9, 11, 12, 13, 15, 16, 17, 21, 22]
+
+    test_on_library(problems, dicts, alg_ops, broken_tests)
+end
+
+@testset "SimpleBroyden 23 Test Problems" begin
+    alg_ops = (SimpleBroyden(),)
+
+    broken_tests = Dict(alg => Int[] for alg in alg_ops)
+    broken_tests[alg_ops[1]] = [1, 4, 5, 6, 11, 12, 13, 14]
+
+    skip_tests = Dict(alg => Int[] for alg in alg_ops)
+    skip_tests[alg_ops[1]] = [2, 22]
+
+    test_on_library(problems, dicts, alg_ops, broken_tests; skip_tests)
+end
+
+@testset "SimpleKlement 23 Test Problems" begin
+    alg_ops = (SimpleKlement(),)
+
+    broken_tests = Dict(alg => Int[] for alg in alg_ops)
+    broken_tests[alg_ops[1]] = [1, 2, 4, 5, 6, 7, 9, 10, 11, 12, 13, 19, 21, 22]
+
+    test_on_library(problems, dicts, alg_ops, broken_tests)
+end
diff --git a/test/Project.toml b/test/Project.toml
index 469f302..230ab90 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -1,10 +1,13 @@
 [deps]
+AllocCheck = "9b6a8646-10ed-4001-bbdc-1d2f46dfbb1a"
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-NNlib = "872c559c-99b0-510c-b3b7-b6c96a88d5cd"
+LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
diff --git a/test/basictests.jl b/test/basictests.jl
index 027f766..6b70403 100644
--- a/test/basictests.jl
+++ b/test/basictests.jl
@@ -1,578 +1,366 @@
-using SimpleNonlinearSolve,
-    StaticArrays, BenchmarkTools, DiffEqBase, LinearAlgebra, Test,
-    NNlib
-
-const BATCHED_BROYDEN_SOLVERS = []
-const BROYDEN_SOLVERS = []
-const BATCHED_LBROYDEN_SOLVERS = []
-const LBROYDEN_SOLVERS = []
-const BATCHED_DFSANE_SOLVERS = []
-const DFSANE_SOLVERS = []
-const BATCHED_RAPHSON_SOLVERS = []
-
-for mode in instances(NLSolveTerminationMode.T)
-    if mode ∈
-       (NLSolveTerminationMode.SteadyStateDefault, NLSolveTerminationMode.RelSafeBest,
-        NLSolveTerminationMode.AbsSafeBest)
-        continue
-    end
-
-    termination_condition = NLSolveTerminationCondition(mode; abstol = nothing,
-        reltol = nothing)
-    push!(BROYDEN_SOLVERS, Broyden(; batched = false, termination_condition))
-    push!(BATCHED_BROYDEN_SOLVERS, Broyden(; batched = true, termination_condition))
-    push!(LBROYDEN_SOLVERS, LBroyden(; batched = false, termination_condition))
-    push!(BATCHED_LBROYDEN_SOLVERS, LBroyden(; batched = true, termination_condition))
-    push!(DFSANE_SOLVERS, SimpleDFSane(; batched = false, termination_condition))
-    push!(BATCHED_DFSANE_SOLVERS, SimpleDFSane(; batched = true, termination_condition))
-    push!(BATCHED_RAPHSON_SOLVERS,
-        SimpleNewtonRaphson(; batched = true,
-            termination_condition))
-    push!(BATCHED_RAPHSON_SOLVERS,
-        SimpleNewtonRaphson(; batched = true, autodiff = false,
-            termination_condition))
-end
-
-# SimpleNewtonRaphson
-function benchmark_scalar(f, u0)
-    probN = NonlinearProblem{false}(f, u0)
-    sol = (solve(probN, SimpleNewtonRaphson()))
-end
-
-function ff(u, p)
-    u .* u .- 2
-end
-const cu0 = @SVector[1.0, 1.0]
-function sf(u, p)
-    u * u - 2
-end
-const csu0 = 1.0
-
-sol = benchmark_scalar(sf, csu0)
-@test sol.retcode === ReturnCode.Success
-@test sol.u * sol.u - 2 < 1e-9
-
-if VERSION >= v"1.7"
-    @test (@ballocated benchmark_scalar(sf, csu0)) == 0
-end
-
-# SimpleHalley
-function benchmark_scalar(f, u0)
-    probN = NonlinearProblem{false}(f, u0)
-    sol = (solve(probN, SimpleHalley()))
-end
+using AllocCheck, BenchmarkTools, LinearSolve, SimpleNonlinearSolve, StaticArrays, Random,
+    LinearAlgebra, Test, ForwardDiff, DiffEqBase
+
+_nameof(x) = applicable(nameof, x) ? nameof(x) : _nameof(typeof(x))
+
+quadratic_f(u, p) = u .* u .- p
+quadratic_f!(du, u, p) = (du .= u .* u .- p)
+quadratic_f2(u, p) = @. p[1] * u * u - p[2]
+
+function newton_fails(u, p)
+    return 0.010000000000000002 .+
+           10.000000000000002 ./ (1 .+
+            (0.21640425613334457 .+
+             216.40425613334457 ./ (1 .+
+              (0.21640425613334457 .+
+               216.40425613334457 ./
+               (1 .+ 0.0006250000000000001(u .^ 2.0))) .^ 2.0)) .^ 2.0) .-
+           0.0011552453009332421u .- p
+end
+
+const TERMINATION_CONDITIONS = [
+    NormTerminationMode(), RelTerminationMode(), RelNormTerminationMode(),
+    AbsTerminationMode(), AbsNormTerminationMode(), RelSafeTerminationMode(),
+    AbsSafeTerminationMode(), RelSafeBestTerminationMode(), AbsSafeBestTerminationMode(),
+]
 
-function ff(u, p)
-    u .* u .- 2
-end
-const cu0 = @SVector[1.0, 1.0]
-function sf(u, p)
-    u * u - 2
-end
-const csu0 = 1.0
+# --- SimpleNewtonRaphson tests ---
 
-sol = benchmark_scalar(sf, csu0)
-@test sol.retcode === ReturnCode.Success
-@test sol.u * sol.u - 2 < 1e-9
+@testset "$(alg)" for alg in (SimpleNewtonRaphson, SimpleTrustRegion)
+    # Eval else the alg is type unstable
+    @eval begin
+        function benchmark_nlsolve_oop(f, u0, p = 2.0; autodiff = AutoForwardDiff())
+            prob = NonlinearProblem{false}(f, u0, p)
+            return solve(prob, $(alg)(; autodiff), abstol = 1e-9)
+        end
 
-sol = benchmark_scalar(ff, cu0)
-@test sol.retcode === ReturnCode.Success
-@test sol.u .* sol.u .- 2 < [1e-9, 1e-9]
+        function benchmark_nlsolve_iip(f, u0, p = 2.0; autodiff = AutoForwardDiff())
+            prob = NonlinearProblem{true}(f, u0, p)
+            return solve(prob, $(alg)(; autodiff), abstol = 1e-9)
+        end
+    end
 
-if VERSION >= v"1.7"
-    @test (@ballocated benchmark_scalar(sf, csu0)) == 0
-end
+    @testset "AutoDiff: $(_nameof(autodiff))" for autodiff in (AutoFiniteDiff(),
+        AutoForwardDiff())
+        @testset "[OOP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; autodiff)
+            @test SciMLBase.successful_retcode(sol)
+            @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+        end
 
-# Broyden
-function benchmark_scalar(f, u0, alg)
-    probN = NonlinearProblem{false}(f, u0)
-    sol = (solve(probN, alg))
-end
+        @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
+            sol = benchmark_nlsolve_iip(quadratic_f!, u0; autodiff)
+            @test SciMLBase.successful_retcode(sol)
+            @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+        end
+    end
 
-for alg in BROYDEN_SOLVERS
-    sol = benchmark_scalar(sf, csu0, alg)
-    @test sol.retcode === ReturnCode.Success
-    @test sol.u * sol.u - 2 < 1e-9
-    # FIXME: Termination Condition Implementation is allocating. Not sure how to fix it.
-    # if VERSION >= v"1.7"
-    #     @test (@ballocated benchmark_scalar($sf, $csu0, $termination_condition)) == 0
-    # end
-end
+    @testset "Allocations: Static Array and Scalars" begin
+        @test (@ballocated $(benchmark_nlsolve_oop)($quadratic_f, $(@SVector[1.0, 1.0]),
+            2.0; autodiff = AutoForwardDiff())) < 200
+        @test (@ballocated $(benchmark_nlsolve_oop)($quadratic_f, 1.0, 2.0;
+            autodiff = AutoForwardDiff())) == 0
+    end
 
-# Klement
-function benchmark_scalar(f, u0)
-    probN = NonlinearProblem{false}(f, u0)
-    sol = (solve(probN, Klement()))
-end
+    @testset "[OOP] Immutable AD" begin
+        for p in [1.0, 100.0]
+            @test begin
+                res = benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p)
+                res_true = sqrt(p)
+                all(res.u .≈ res_true)
+            end
+            @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                @SVector[1.0, 1.0], p).u[end], p) ≈ 1 / (2 * sqrt(p))
+        end
+    end
 
-sol = benchmark_scalar(sf, csu0)
-@test sol.retcode === ReturnCode.Success
-@test sol.u * sol.u - 2 < 1e-9
-if VERSION >= v"1.7"
-    @test (@ballocated benchmark_scalar(sf, csu0)) == 0
-end
+    @testset "[OOP] Scalar AD" begin
+        for p in 1.0:0.1:100.0
+            @test begin
+                res = benchmark_nlsolve_oop(quadratic_f, 1.0, p)
+                res_true = sqrt(p)
+                res.u ≈ res_true
+            end
+            @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, 1.0, p).u,
+                p) ≈ 1 / (2 * sqrt(p))
+        end
+    end
 
-# SimpleTrustRegion
-function benchmark_scalar(f, u0)
-    probN = NonlinearProblem{false}(f, u0)
-    sol = (solve(probN, SimpleTrustRegion()))
-end
+    t = (p) -> [sqrt(p[2] / p[1])]
+    p = [0.9, 50.0]
+    @test benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u ≈ sqrt(p[2] / p[1])
+    @test ForwardDiff.jacobian(p -> [benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u],
+        p) ≈ ForwardDiff.jacobian(t, p)
 
-sol = benchmark_scalar(sf, csu0)
-@test sol.retcode === ReturnCode.Success
-@test sol.u * sol.u - 2 < 1e-9
+    @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS,
+        u0 in (1.0, [1.0, 1.0], @SVector[1.0, 1.0])
 
-# SimpleDFSane
-function benchmark_scalar(f, u0)
-    probN = NonlinearProblem{false}(f, u0)
-    sol = (solve(probN, SimpleDFSane()))
+        probN = NonlinearProblem(quadratic_f, u0, 2.0)
+        @test all(solve(probN, alg(); termination_condition).u .≈ sqrt(2.0))
+    end
 end
 
-sol = benchmark_scalar(sf, csu0)
-@test sol.retcode === ReturnCode.Success
-@test sol.u * sol.u - 2 < 1e-9
-
-# AD Tests
-using ForwardDiff
-
-# Immutable
-f, u0 = (u, p) -> u .* u .- p, @SVector[1.0, 1.0]
+# --- SimpleHalley tests ---
 
-for alg in (SimpleNewtonRaphson(), LBroyden(), Klement(), SimpleTrustRegion(),
-    SimpleDFSane(), SimpleHalley(), BROYDEN_SOLVERS...)
-    g = function (p)
-        probN = NonlinearProblem{false}(f, csu0, p)
-        sol = solve(probN, alg, abstol = 1e-9)
-        return sol.u[end]
+@testset "SimpleHalley" begin
+    function benchmark_nlsolve_oop(f, u0, p = 2.0; autodiff = AutoForwardDiff())
+        prob = NonlinearProblem{false}(f, u0, p)
+        return solve(prob, SimpleHalley(; autodiff), abstol = 1e-9)
     end
 
-    for p in 1.1:0.1:100.0
-        res = abs.(g(p))
-        # Not surprising if LBrouden fails to converge
-        if any(x -> isnan(x) || x <= 1e-5 || x >= 1e5, res) && alg isa LBroyden
-            @test_broken res ≈ sqrt(p)
-            @test_broken abs.(ForwardDiff.derivative(g, p)) ≈ 1 / (2 * sqrt(p))
-        else
-            @test res ≈ sqrt(p)
-            @test abs.(ForwardDiff.derivative(g, p)) ≈ 1 / (2 * sqrt(p))
+    @testset "AutoDiff: $(_nameof(autodiff))" for autodiff in (AutoFiniteDiff(),
+        AutoForwardDiff())
+        @testset "[OOP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
+            sol = benchmark_nlsolve_oop(quadratic_f, u0; autodiff)
+            @test SciMLBase.successful_retcode(sol)
+            @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
         end
     end
-end
 
-# Scalar
-f, u0 = (u, p) -> u * u - p, 1.0
-for alg in (SimpleNewtonRaphson(), Klement(), SimpleTrustRegion(),
-    SimpleDFSane(), SimpleHalley(), BROYDEN_SOLVERS..., LBROYDEN_SOLVERS...)
-    g = function (p)
-        probN = NonlinearProblem{false}(f, oftype(p, u0), p)
-        sol = solve(probN, alg)
-        return sol.u
+    @testset "Allocations: Static Array and Scalars" begin
+        @test (@ballocated $(benchmark_nlsolve_oop)($quadratic_f, 1.0, 2.0;
+            autodiff = AutoForwardDiff())) == 0
     end
 
-    for p in 1.1:0.1:100.0
-        res = abs.(g(p))
-        # Not surprising if LBrouden fails to converge
-        if any(x -> isnan(x) || x <= 1e-5 || x >= 1e5, res) && alg isa LBroyden
-            @test_broken res ≈ sqrt(p)
-            @test_broken abs.(ForwardDiff.derivative(g, p)) ≈ 1 / (2 * sqrt(p))
-        else
-            @test res ≈ sqrt(p)
-            @test abs.(ForwardDiff.derivative(g, p)) ≈ 1 / (2 * sqrt(p))
+    @testset "[OOP] Immutable AD" begin
+        for p in [1.0, 100.0]
+            @test begin
+                res = benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p)
+                res_true = sqrt(p)
+                all(res.u .≈ res_true)
+            end
+            @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                @SVector[1.0, 1.0], p).u[end], p) ≈ 1 / (2 * sqrt(p))
         end
     end
-end
-
-tspan = (1.0, 20.0)
-# Falsi
-g = function (p)
-    probN = IntervalNonlinearProblem{false}(f, typeof(p).(tspan), p)
-    sol = solve(probN, Falsi())
-    return sol.left
-end
-
-for p in 1.1:0.1:100.0
-    @test g(p) ≈ sqrt(p)
-    @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
-end
 
-# Ridder
-g = function (p)
-    probN = IntervalNonlinearProblem{false}(f, typeof(p).(tspan), p)
-    sol = solve(probN, Ridder())
-    return sol.left
-end
+    @testset "[OOP] Scalar AD" begin
+        for p in 1.0:0.1:100.0
+            @test begin
+                res = benchmark_nlsolve_oop(quadratic_f, 1.0, p)
+                res_true = sqrt(p)
+                res.u ≈ res_true
+            end
+            @test ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f, 1.0, p).u,
+                p) ≈ 1 / (2 * sqrt(p))
+        end
+    end
 
-for p in 1.1:0.1:100.0
-    @test g(p) ≈ sqrt(p)
-    @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
-end
+    t = (p) -> [sqrt(p[2] / p[1])]
+    p = [0.9, 50.0]
+    @test benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u ≈ sqrt(p[2] / p[1])
+    @test ForwardDiff.jacobian(p -> [benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u],
+        p) ≈ ForwardDiff.jacobian(t, p)
 
-# Brent
-g = function (p)
-    probN = IntervalNonlinearProblem{false}(f, typeof(p).(tspan), p)
-    sol = solve(probN, Brent())
-    return sol.left
-end
+    @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS,
+        u0 in (1.0, [1.0, 1.0], @SVector[1.0, 1.0])
 
-for p in 1.1:0.1:100.0
-    @test g(p) ≈ sqrt(p)
-    @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
+        probN = NonlinearProblem(quadratic_f, u0, 2.0)
+        @test all(solve(probN, SimpleHalley(); termination_condition).u .≈ sqrt(2.0))
+    end
 end
 
-# ITP
-g = function (p)
-    probN = IntervalNonlinearProblem{false}(f, typeof(p).(tspan), p)
-    sol = solve(probN, ITP())
-    return sol.u
-end
+# --- SimpleBroyden / SimpleKlement / SimpleLimitedMemoryBroyden tests ---
 
-for p in 1.1:0.1:100.0
-    @test g(p) ≈ sqrt(p)
-    @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
-end
+@testset "$(_nameof(alg))" for alg in [SimpleBroyden(), SimpleKlement(), SimpleDFSane(),
+    SimpleLimitedMemoryBroyden()]
+    function benchmark_nlsolve_oop(f, u0, p = 2.0)
+        prob = NonlinearProblem{false}(f, u0, p)
+        return solve(prob, alg, abstol = 1e-9)
+    end
 
-# Alefeld
-g = function (p)
-    probN = IntervalNonlinearProblem{false}(f, typeof(p).(tspan), p)
-    sol = solve(probN, Alefeld())
-    return sol.u
-end
+    function benchmark_nlsolve_iip(f, u0, p = 2.0)
+        prob = NonlinearProblem{true}(f, u0, p)
+        return solve(prob, alg, abstol = 1e-9)
+    end
 
-for p in 1.1:0.1:100.0
-    @test g(p) ≈ sqrt(p)
-    @test ForwardDiff.derivative(g, p) ≈ 1 / (2 * sqrt(p))
-end
+    @testset "[OOP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0], @SVector[1.0, 1.0], 1.0)
+        sol = benchmark_nlsolve_oop(quadratic_f, u0)
+        @test SciMLBase.successful_retcode(sol)
+        @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+    end
 
-f, tspan = (u, p) -> p[1] * u * u - p[2], (1.0, 100.0)
-t = (p) -> [sqrt(p[2] / p[1])]
-p = [0.9, 50.0]
-g = function (p)
-    probN = IntervalNonlinearProblem{false}(f, tspan, p)
-    sol = solve(probN, Alefeld())
-    return [sol.u]
-end
+    @testset "[IIP] u0: $(typeof(u0))" for u0 in ([1.0, 1.0],)
+        sol = benchmark_nlsolve_iip(quadratic_f!, u0)
+        @test SciMLBase.successful_retcode(sol)
+        @test all(abs.(sol.u .* sol.u .- 2) .< 1e-9)
+    end
 
-@test g(p) ≈ [sqrt(p[2] / p[1])]
-@test ForwardDiff.jacobian(g, p) ≈ ForwardDiff.jacobian(t, p)
-
-f, tspan = (u, p) -> p[1] * u * u - p[2], (1.0, 100.0)
-t = (p) -> [sqrt(p[2] / p[1])]
-p = [0.9, 50.0]
-for alg in [Bisection(), Falsi(), Ridder(), Brent(), ITP()]
-    global g, p
-    g = function (p)
-        probN = IntervalNonlinearProblem{false}(f, tspan, p)
-        sol = solve(probN, alg)
-        return [sol.left]
+    @testset "Allocations: Static Array and Scalars" begin
+        @test (@ballocated $(benchmark_nlsolve_oop)($quadratic_f, $(@SVector[1.0, 1.0]),
+            2.0)) < 200
+        allocs = alg isa SimpleDFSane ? 144 : 0
+        @test (@ballocated $(benchmark_nlsolve_oop)($quadratic_f, 1.0, 2.0)) == allocs
     end
 
-    @test g(p) ≈ [sqrt(p[2] / p[1])]
-    @test ForwardDiff.jacobian(g, p) ≈ ForwardDiff.jacobian(t, p)
-end
+    @testset "[OOP] Immutable AD" begin
+        for p in [1.0, 100.0]
+            res = benchmark_nlsolve_oop(quadratic_f, @SVector[1.0, 1.0], p)
+
+            if any(x -> isnan(x) || x <= 1e-5 || x >= 1e5, res)
+                @test_broken all(res .≈ sqrt(p))
+                @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                    @SVector[1.0, 1.0], p).u[end], p)) ≈ 1 / (2 * sqrt(p))
+            else
+                @test all(res .≈ sqrt(p))
+                @test isapprox(abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                        @SVector[1.0, 1.0], p).u[end], p)), 1 / (2 * sqrt(p)))
+            end
+        end
+    end
 
-for alg in (SimpleNewtonRaphson(), Klement(), SimpleTrustRegion(),
-    SimpleDFSane(), SimpleHalley(), BROYDEN_SOLVERS..., LBROYDEN_SOLVERS...)
-    global g, p
-    g = function (p)
-        probN = NonlinearProblem{false}(f, 0.5, p)
-        sol = solve(probN, alg)
-        return [abs(sol.u)]
+    @testset "[OOP] Scalar AD" begin
+        for p in 1.0:0.1:100.0
+            res = benchmark_nlsolve_oop(quadratic_f, 1.0, p)
+
+            if any(x -> isnan(x) || x <= 1e-5 || x >= 1e5, res)
+                @test_broken all(res .≈ sqrt(p))
+                @test_broken abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                        1.0, p).u, p)) ≈ 1 / (2 * sqrt(p))
+            else
+                @test all(res .≈ sqrt(p))
+                @test isapprox(abs.(ForwardDiff.derivative(p -> benchmark_nlsolve_oop(quadratic_f,
+                            1.0, p).u, p)), 1 / (2 * sqrt(p)))
+            end
+        end
     end
-    @test g(p) ≈ [sqrt(p[2] / p[1])]
-    @test ForwardDiff.jacobian(g, p) ≈ ForwardDiff.jacobian(t, p)
-end
 
-# Error Checks
-f, u0 = (u, p) -> u .* u .- 2.0, @SVector[1.0, 1.0]
-probN = NonlinearProblem(f, u0)
+    t = (p) -> [sqrt(p[2] / p[1])]
+    p = [0.9, 50.0]
+    @test benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u ≈ sqrt(p[2] / p[1])
+    @test ForwardDiff.jacobian(p -> [benchmark_nlsolve_oop(quadratic_f2, 0.5, p).u],
+        p) ≈ ForwardDiff.jacobian(t, p)
 
-for alg in (SimpleNewtonRaphson(), SimpleNewtonRaphson(; autodiff = false),
-    SimpleTrustRegion(),
-    SimpleTrustRegion(; autodiff = false), SimpleHalley(), SimpleHalley(; autodiff = false),
-    Klement(), SimpleDFSane(),
-    BROYDEN_SOLVERS..., LBROYDEN_SOLVERS...)
-    sol = solve(probN, alg)
+    @testset "Termination condition: $(termination_condition) u0: $(_nameof(u0))" for termination_condition in TERMINATION_CONDITIONS,
+        u0 in (1.0, [1.0, 1.0], @SVector[1.0, 1.0])
 
-    @test sol.retcode == ReturnCode.Success
-    @test sol.u[end] ≈ sqrt(2.0)
+        probN = NonlinearProblem(quadratic_f, u0, 2.0)
+        @test all(solve(probN, alg; termination_condition).u .≈ sqrt(2.0))
+    end
 end
 
-for u0 in [1.0, [1, 1.0]]
-    local f, probN, sol
-    f = (u, p) -> u .* u .- 2.0
-    probN = NonlinearProblem(f, u0)
-    sol = sqrt(2) * u0
+@testset "Newton Fails" begin
+    function benchmark_nlsolve_oop(f, u0, p, alg)
+        prob = NonlinearProblem{false}(f, u0, p)
+        return solve(prob, alg; abstol = 1e-9)
+    end
 
-    for alg in (SimpleNewtonRaphson(), SimpleNewtonRaphson(; autodiff = false),
-        SimpleTrustRegion(), SimpleTrustRegion(; autodiff = false), Klement(),
-        SimpleDFSane(), BROYDEN_SOLVERS..., LBROYDEN_SOLVERS...)
-        sol2 = solve(probN, alg)
+    u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
+    p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
 
-        @test sol2.retcode == ReturnCode.Success
-        @test sol2.u ≈ sol
+    for alg in (SimpleDFSane(), SimpleTrustRegion(), SimpleHalley())
+        sol = benchmark_nlsolve_oop(newton_fails, u0, p, alg)
+        @test SciMLBase.successful_retcode(sol)
+        @test all(abs.(newton_fails(sol.u, p)) .< 1e-9)
     end
 end
 
-# Bisection Tests
-f, tspan = (u, p) -> u .* u .- 2.0, (1.0, 2.0)
-probB = IntervalNonlinearProblem(f, tspan)
-
-# Falsi
-sol = solve(probB, Falsi())
-@test sol.left ≈ sqrt(2.0)
-
-# Bisection
-sol = solve(probB, Bisection())
-@test sol.left ≈ sqrt(2.0)
-
-# Ridder
-sol = solve(probB, Ridder())
-@test sol.left ≈ sqrt(2.0)
-tspan = (sqrt(2.0), 10.0)
-probB = IntervalNonlinearProblem(f, tspan)
-sol = solve(probB, Ridder())
-@test sol.left ≈ sqrt(2.0)
-tspan = (0.0, sqrt(2.0))
-probB = IntervalNonlinearProblem(f, tspan)
-sol = solve(probB, Ridder())
-@test sol.left ≈ sqrt(2.0)
-
-# Brent
-sol = solve(probB, Brent())
-@test sol.left ≈ sqrt(2.0)
-tspan = (sqrt(2.0), 10.0)
-probB = IntervalNonlinearProblem(f, tspan)
-sol = solve(probB, Brent())
-@test sol.left ≈ sqrt(2.0)
-tspan = (0.0, sqrt(2.0))
-probB = IntervalNonlinearProblem(f, tspan)
-sol = solve(probB, Brent())
-@test sol.left ≈ sqrt(2.0)
-
-# Alefeld
-sol = solve(probB, Alefeld())
-@test sol.u ≈ sqrt(2.0)
-tspan = (sqrt(2.0), 10.0)
-probB = IntervalNonlinearProblem(f, tspan)
-sol = solve(probB, Alefeld())
-@test sol.u ≈ sqrt(2.0)
-tspan = (0.0, sqrt(2.0))
-probB = IntervalNonlinearProblem(f, tspan)
-sol = solve(probB, Alefeld())
-@test sol.u ≈ sqrt(2.0)
-
-# ITP
-sol = solve(probB, ITP())
-@test sol.u ≈ sqrt(2.0)
-tspan = (sqrt(2.0), 10.0)
-probB = IntervalNonlinearProblem(f, tspan)
-sol = solve(probB, ITP())
-@test sol.u ≈ sqrt(2.0)
-tspan = (0.0, sqrt(2.0))
-probB = IntervalNonlinearProblem(f, tspan)
-sol = solve(probB, ITP())
-@test sol.u ≈ sqrt(2.0)
-
-# Tolerance tests for Interval methods
-f, tspan = (u, p) -> u .* u .- 2.0, (1.0, 10.0)
-probB = IntervalNonlinearProblem(f, tspan)
-tols = [0.1, 0.01, 0.001, 0.0001, 1e-5, 1e-6, 1e-7]
-ϵ = eps(1.0) #least possible tol for all methods
-
-for atol in tols
-    sol = solve(probB, Bisection(), abstol = atol)
-    @test abs(sol.u - sqrt(2)) < atol
-    @test abs(sol.u - sqrt(2)) > ϵ #test that the solution is not calculated upto max precision
-    sol = solve(probB, Falsi(), abstol = atol)
-    @test abs(sol.u - sqrt(2)) < atol
-    @test abs(sol.u - sqrt(2)) > ϵ
-    sol = solve(probB, ITP(), abstol = atol)
-    @test abs(sol.u - sqrt(2)) < atol
-    @test abs(sol.u - sqrt(2)) > ϵ
-end
+# --- Allocation Checks ---
 
-tols = [0.1] # Ridder and Brent converge rapidly so as we lower tolerance below 0.01, it converges with max precision to the solution
-for atol in tols
-    sol = solve(probB, Ridder(), abstol = atol)
-    @test abs(sol.u - sqrt(2)) < atol
-    @test abs(sol.u - sqrt(2)) > ϵ
-    sol = solve(probB, Brent(), abstol = atol)
-    @test abs(sol.u - sqrt(2)) < atol
-    @test abs(sol.u - sqrt(2)) > ϵ
-end
+## SimpleDFSane needs to allocate a history vector
+@testset "Allocation Checks: $(_nameof(alg))" for alg in (SimpleNewtonRaphson(),
+    SimpleHalley(), SimpleBroyden(), SimpleKlement(), SimpleLimitedMemoryBroyden(),
+    SimpleTrustRegion())
+    @check_allocs nlsolve(prob, alg) = DiffEqBase.__solve(prob, alg; abstol = 1e-9)
 
-# Garuntee Tests for Bisection
-f = function (u, p)
-    if u < 2.0
-        return u - 2.0
-    elseif u > 3.0
-        return u - 3.0
-    else
-        return 0.0
+    nlprob_scalar = NonlinearProblem{false}(quadratic_f, 1.0, 2.0)
+    nlprob_sa = NonlinearProblem{false}(quadratic_f, @SVector[1.0, 1.0], 2.0)
+
+    try
+        nlsolve(nlprob_scalar, alg)
+        @test true
+    catch e
+        @error e
+        @test false
     end
-end
-probB = IntervalNonlinearProblem(f, (0.0, 4.0))
-
-sol = solve(probB, Bisection(; exact_left = true))
-@test f(sol.left, nothing) < 0.0
-@test f(nextfloat(sol.left), nothing) >= 0.0
-
-sol = solve(probB, Bisection(; exact_right = true))
-@test f(sol.right, nothing) >= 0.0
-@test f(prevfloat(sol.right), nothing) <= 0.0
-
-sol = solve(probB, Bisection(; exact_left = true, exact_right = true); immutable = false)
-@test f(sol.left, nothing) < 0.0
-@test f(nextfloat(sol.left), nothing) >= 0.0
-@test f(sol.right, nothing) >= 0.0
-@test f(prevfloat(sol.right), nothing) <= 0.0
-
-# Test that `SimpleTrustRegion` passes a test that `SimpleNewtonRaphson` fails on.
-u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
-global g, f
-f = (u, p) -> 0.010000000000000002 .+
-              10.000000000000002 ./ (1 .+
-               (0.21640425613334457 .+
-                216.40425613334457 ./ (1 .+
-                 (0.21640425613334457 .+
-                  216.40425613334457 ./
-                  (1 .+ 0.0006250000000000001(u .^ 2.0))) .^ 2.0)) .^ 2.0) .-
-              0.0011552453009332421u .- p
-g = function (p)
-    probN = NonlinearProblem{false}(f, u0, p)
-    sol = solve(probN, SimpleTrustRegion())
-    return sol.u
-end
-p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-u = g(p)
-f(u, p)
-@test all(abs.(f(u, p)) .< 1e-10)
-
-# Test kwars in `SimpleTrustRegion`
-max_trust_radius = [10.0, 100.0, 1000.0]
-initial_trust_radius = [10.0, 1.0, 0.1]
-step_threshold = [0.0, 0.01, 0.25]
-shrink_threshold = [0.25, 0.3, 0.5]
-expand_threshold = [0.5, 0.8, 0.9]
-shrink_factor = [0.1, 0.3, 0.5]
-expand_factor = [1.5, 2.0, 3.0]
-max_shrink_times = [10, 20, 30]
-
-list_of_options = zip(max_trust_radius, initial_trust_radius, step_threshold,
-    shrink_threshold, expand_threshold, shrink_factor,
-    expand_factor, max_shrink_times)
-for options in list_of_options
-    local probN, sol, alg
-    alg = SimpleTrustRegion(max_trust_radius = options[1],
-        initial_trust_radius = options[2],
-        step_threshold = options[3],
-        shrink_threshold = options[4],
-        expand_threshold = options[5],
-        shrink_factor = options[6],
-        expand_factor = options[7],
-        max_shrink_times = options[8])
-
-    probN = NonlinearProblem(f, u0, p)
-    sol = solve(probN, alg)
-    @test all(abs.(f(u, p)) .< 1e-10)
-end
 
-# Test that `SimpleDFSane` passes a test that `SimpleNewtonRaphson` fails on.
-u0 = [-10.0, -1.0, 1.0, 2.0, 3.0, 4.0, 10.0]
-global g, f
-f = (u, p) -> 0.010000000000000002 .+
-              10.000000000000002 ./ (1 .+
-               (0.21640425613334457 .+
-                216.40425613334457 ./ (1 .+
-                 (0.21640425613334457 .+
-                  216.40425613334457 ./
-                  (1 .+ 0.0006250000000000001(u .^ 2.0))) .^ 2.0)) .^ 2.0) .-
-              0.0011552453009332421u .- p
-g = function (p)
-    probN = NonlinearProblem{false}(f, u0, p)
-    sol = solve(probN, SimpleDFSane())
-    return sol.u
+    # ForwardDiff allocates for hessian since we don't propagate the chunksize
+    # SimpleLimitedMemoryBroyden needs to do views on the low rank matrices so the sizes
+    # are dynamic. This can be fixed but no without maintaining the simplicity of the code
+    try
+        nlsolve(nlprob_sa, alg)
+        @test true
+    catch e
+        @error e
+        @test false broken=(alg isa SimpleHalley || alg isa SimpleLimitedMemoryBroyden)
+    end
 end
-p = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
-u = g(p)
-f(u, p)
-@test all(abs.(f(u, p)) .< 1e-10)
-
-# Test kwars in `SimpleDFSane`
-σ_min = [1e-10, 1e-5, 1e-4]
-σ_max = [1e10, 1e5, 1e4]
-σ_1 = [1.0, 0.5, 2.0]
-M = [10, 1, 100]
-γ = [1e-4, 1e-3, 1e-5]
-τ_min = [0.1, 0.2, 0.3]
-τ_max = [0.5, 0.8, 0.9]
-nexp = [2, 1, 2]
-η_strategy = [
-    (f_1, k, x, F) -> f_1 / k^2,
-    (f_1, k, x, F) -> f_1 / k^3,
-    (f_1, k, x, F) -> f_1 / k^4,
-]
 
-list_of_options = zip(σ_min, σ_max, σ_1, M, γ, τ_min, τ_max, nexp,
-    η_strategy)
-for options in list_of_options
-    local probN, sol, alg
-    alg = SimpleDFSane(σ_min = options[1],
-        σ_max = options[2],
-        σ_1 = options[3],
-        M = options[4],
-        γ = options[5],
-        τ_min = options[6],
-        τ_max = options[7],
-        nexp = options[8],
-        η_strategy = options[9])
-
-    probN = NonlinearProblem(f, u0, p)
-    sol = solve(probN, alg)
-    @test all(abs.(f(u, p)) .< 1e-10)
-end
+# --- Interval Nonlinear Problems ---
 
-f, u0 = (u, p) -> u .* u .- p, randn(1, 3)
+@testset "Interval Nonlinear Problem: $(alg)" for alg in (Bisection(), Falsi(), Ridder(),
+    Brent(), ITP(), Alefeld())
+    tspan = (1.0, 20.0)
 
-p = [2.0 1.0 5.0];
-probN = NonlinearProblem{false}(f, u0, p);
+    function g(p)
+        probN = IntervalNonlinearProblem{false}(quadratic_f, typeof(p).(tspan), p)
+        sol = solve(probN, alg; abstol = 1e-9)
+        return sol.left
+    end
 
-sol = solve(probN, Broyden(batched = true))
+    for p in 1.1:0.1:100.0
+        @test g(p)≈sqrt(p) atol=1e-3 rtol=1e-3
+        @test ForwardDiff.derivative(g, p)≈1 / (2 * sqrt(p)) atol=1e-3 rtol=1e-3
+    end
 
-@test abs.(sol.u) ≈ sqrt.(p)
+    t = (p) -> [sqrt(p[2] / p[1])]
+    p = [0.9, 50.0]
 
-@testset "Batched Solver: $(nameof(typeof(alg)))" for alg in (BATCHED_BROYDEN_SOLVERS...,
-    BATCHED_LBROYDEN_SOLVERS...,
-    BATCHED_DFSANE_SOLVERS...,
-    BATCHED_RAPHSON_SOLVERS...)
-    sol = solve(probN, alg; abstol = 1e-3, reltol = 1e-3)
+    function g2(p)
+        probN = IntervalNonlinearProblem{false}((u, p) -> p[1] * u * u - p[2], tspan, p)
+        sol = solve(probN, alg; abstol = 1e-9)
+        return [sol.u]
+    end
 
-    @test sol.retcode == ReturnCode.Success
-    @test abs.(sol.u)≈sqrt.(p) atol=1e-3 rtol=1e-3
-end
+    @test g2(p)≈[sqrt(p[2] / p[1])] atol=1e-3 rtol=1e-3
+    @test ForwardDiff.jacobian(g2, p)≈ForwardDiff.jacobian(t, p) atol=1e-3 rtol=1e-3
 
-## User specified Jacobian
+    probB = IntervalNonlinearProblem{false}(quadratic_f, (1.0, 2.0), 2.0)
+    sol = solve(probB, alg; abstol = 1e-9)
+    @test sol.left≈sqrt(2.0) atol=1e-3 rtol=1e-3
 
-f, u0 = (u, p) -> u .* u .- p, randn(3)
+    if !(alg isa Bisection || alg isa Falsi)
+        probB = IntervalNonlinearProblem{false}(quadratic_f, (sqrt(2.0), 10.0), 2.0)
+        sol = solve(probB, alg; abstol = 1e-9)
+        @test sol.left≈sqrt(2.0) atol=1e-3 rtol=1e-3
 
-f_jac(u, p) = begin
-    diagm(2 * u)
+        probB = IntervalNonlinearProblem{false}(quadratic_f, (0.0, sqrt(2.0)), 2.0)
+        sol = solve(probB, alg; abstol = 1e-9)
+        @test sol.left≈sqrt(2.0) atol=1e-3 rtol=1e-3
+    end
 end
 
-p = [2.0, 1.0, 5.0];
+@testset "Tolerance Tests Interval Methods: $(alg)" for alg in (Bisection(), Falsi(), ITP())
+    tspan = (1.0, 20.0)
+    probB = IntervalNonlinearProblem(quadratic_f, tspan, 2.0)
+    tols = [0.1, 0.01, 0.001, 0.0001, 1e-5, 1e-6, 1e-7]
+    ϵ = eps(1.0) #least possible tol for all methods
 
-probN = NonlinearProblem(NonlinearFunction(f, jac = f_jac), u0, p)
+    for atol in tols
+        sol = solve(probB, alg; abstol = atol)
+        @test abs(sol.u - sqrt(2)) < atol
+        @test abs(sol.u - sqrt(2)) > ϵ #test that the solution is not calculated upto max precision
+    end
+end
+
+@testset "Tolerance Tests Interval Methods: $(alg)" for alg in (Ridder(), Brent())
+    tspan = (1.0, 20.0)
+    probB = IntervalNonlinearProblem(quadratic_f, tspan, 2.0)
+    tols = [0.1] # Ridder and Brent converge rapidly so as we lower tolerance below 0.01, it converges with max precision to the solution
+    ϵ = eps(1.0) #least possible tol for all methods
 
-for alg in (SimpleNewtonRaphson(), SimpleTrustRegion())
-    sol = solve(probN, alg)
-    @test abs.(sol.u) ≈ sqrt.(p)
+    for atol in tols
+        sol = solve(probB, alg; abstol = atol)
+        @test abs(sol.u - sqrt(2)) < atol
+        @test abs(sol.u - sqrt(2)) > ϵ #test that the solution is not calculated upto max precision
+    end
 end
 
-# Flipped signs & reversed tspan test for bracketing algorithms
-f1(u, p) = u * u - p
-f2(u, p) = p - u * u
+@testset "Flipped Signs and Reversed Tspan: $(alg)" for alg in (Alefeld(), Bisection(),
+    Falsi(), Brent(), ITP(), Ridder())
+    f1(u, p) = u * u - p
+    f2(u, p) = p - u * u
 
-for alg in (Alefeld(), Bisection(), Falsi(), Brent(), ITP(), Ridder())
     for p in 1:4
         inp1 = IntervalNonlinearProblem(f1, (1.0, 2.0), p)
         inp2 = IntervalNonlinearProblem(f2, (1.0, 2.0), p)
@@ -584,3 +372,31 @@ for alg in (Alefeld(), Bisection(), Falsi(), Brent(), ITP(), Ridder())
         @test abs.(solve(inp4, alg).u) ≈ sqrt.(p)
     end
 end
+
+# The following tests were included in the previos versions but these kwargs never did
+# anything!
+# # Garuntee Tests for Bisection
+# f = function (u, p)
+#     if u < 2.0
+#         return u - 2.0
+#     elseif u > 3.0
+#         return u - 3.0
+#     else
+#         return 0.0
+#     end
+# end
+# probB = IntervalNonlinearProblem(f, (0.0, 4.0))
+
+# sol = solve(probB, Bisection(; exact_left = true))
+# @test f(sol.left, nothing) < 0.0
+# @test f(nextfloat(sol.left), nothing) >= 0.0
+
+# sol = solve(probB, Bisection(; exact_right = true))
+# @test f(sol.right, nothing) >= 0.0
+# @test f(prevfloat(sol.right), nothing) <= 0.0
+
+# sol = solve(probB, Bisection(; exact_left = true, exact_right = true); immutable = false)
+# @test f(sol.left, nothing) < 0.0
+# @test f(nextfloat(sol.left), nothing) >= 0.0
+# @test f(sol.right, nothing) >= 0.0
+# @test f(prevfloat(sol.right), nothing) <= 0.0
diff --git a/test/inplace.jl b/test/inplace.jl
deleted file mode 100644
index 2e9d033..0000000
--- a/test/inplace.jl
+++ /dev/null
@@ -1,52 +0,0 @@
-using SimpleNonlinearSolve,
-    StaticArrays, BenchmarkTools, DiffEqBase, LinearAlgebra, Test,
-    NNlib
-
-# Supported Solvers: BatchedBroyden, BatchedSimpleDFSane, BatchedSimpleNewtonRaphson
-function f!(du::AbstractArray{<:Number, N},
-        u::AbstractArray{<:Number, N},
-        p::AbstractVector) where {N}
-    u_ = reshape(u, :, size(u, N))
-    du .= reshape(sum(abs2, u_; dims = 1) .- u_ .- reshape(p, 1, :), size(u))
-    return du
-end
-
-function f!(du::AbstractMatrix, u::AbstractMatrix, p::AbstractVector)
-    du .= sum(abs2, u; dims = 1) .- u .- reshape(p, 1, :)
-    return du
-end
-
-function f!(du::AbstractVector, u::AbstractVector, p::AbstractVector)
-    du .= sum(abs2, u) .- u .- p
-    return du
-end
-
-@testset "Solver: $(nameof(typeof(solver)))" for solver in (Broyden(; batched = true),
-    SimpleDFSane(; batched = true),
-    SimpleNewtonRaphson(; batched = true))
-    @testset "T: $T" for T in (Float32, Float64)
-        p = rand(T, 5)
-        @testset "size(u0): $sz" for sz in ((2, 5), (1, 5), (2, 3, 5))
-            u0 = ones(T, sz)
-            prob = NonlinearProblem{true}(f!, u0, p)
-
-            sol = solve(prob, solver)
-
-            @test SciMLBase.successful_retcode(sol.retcode)
-
-            @test sol.resid≈zero(sol.resid) atol=5e-3
-        end
-
-        p = rand(T, 1)
-        @testset "size(u0): $sz" for sz in ((3,), (5,), (10,))
-            u0 = ones(T, sz)
-            prob = NonlinearProblem{true}(f!, u0, p)
-
-            sol = solve(prob, solver)
-
-            @test SciMLBase.successful_retcode(sol.retcode)
-
-            @test sol.resid≈zero(sol.resid) atol=5e-3
-        end
-    end
-end
diff --git a/test/least_squares.jl b/test/least_squares.jl
index a7003f6..bc80142 100644
--- a/test/least_squares.jl
+++ b/test/least_squares.jl
@@ -13,7 +13,9 @@ end
 
 θ_init = θ_true .+ 0.1
 prob_oop = NonlinearLeastSquaresProblem{false}(loss_function, θ_init, x)
-sol = solve(prob_oop, SimpleNewtonRaphson())
-sol = solve(prob_oop, SimpleGaussNewton())
 
-@test norm(sol.resid) < 1e-12
+for solver in [SimpleNewtonRaphson(AutoForwardDiff()), SimpleGaussNewton(AutoForwardDiff()),
+    SimpleNewtonRaphson(AutoFiniteDiff()), SimpleGaussNewton(AutoFiniteDiff())]
+    sol = solve(prob_oop, solver)
+    @test norm(sol.resid) < 1e-12
+end
diff --git a/test/matrix_resizing_tests.jl b/test/matrix_resizing_tests.jl
index 9a1989b..66f6a3d 100644
--- a/test/matrix_resizing_tests.jl
+++ b/test/matrix_resizing_tests.jl
@@ -6,6 +6,7 @@ p = 2.0
 vecprob = NonlinearProblem(ff, vec(u0), p)
 prob = NonlinearProblem(ff, u0, p)
 
-for alg in (Klement(), Broyden(), SimpleNewtonRaphson())
-    @test vec(solve(prob, alg).u) == solve(vecprob, alg).u
+@testset "$(alg)" for alg in (SimpleKlement(), SimpleBroyden(), SimpleNewtonRaphson(),
+    SimpleDFSane(), SimpleLimitedMemoryBroyden(), SimpleTrustRegion())
+    @test vec(solve(prob, alg).u) ≈ solve(vecprob, alg).u
 end
diff --git a/test/runtests.jl b/test/runtests.jl
index d0fd1ff..cc4cd70 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,12 +1,12 @@
-using SafeTestsets
+using SafeTestsets, Test
 
 const GROUP = get(ENV, "GROUP", "All")
 
-@time begin
+@time @testset "SimpleNonlinearSolve.jl" begin
     if GROUP == "All" || GROUP == "Core"
         @time @safetestset "Basic Tests + Some AD" include("basictests.jl")
-        @time @safetestset "Inplace Tests" include("inplace.jl")
         @time @safetestset "Matrix Resizing Tests" include("matrix_resizing_tests.jl")
         @time @safetestset "Least Squares Tests" include("least_squares.jl")
+        @time @safetestset "23 Test Problems" include("23_test_problems.jl")
     end
 end