diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 9ca45c89f..3993fff35 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -18,8 +18,8 @@ jobs:
       fail-fast: false
       matrix:
         group:
-          - Core
-          - NLLS
+          - RootFinding
+          - NLLSSolvers
           - 23TestProblems
           - Wrappers
           - Miscellaneous
diff --git a/Project.toml b/Project.toml
index ebc0676d9..3668dfb99 100644
--- a/Project.toml
+++ b/Project.toml
@@ -1,7 +1,7 @@
 name = "NonlinearSolve"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 authors = ["SciML"]
-version = "3.3.0"
+version = "3.4.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
@@ -10,6 +10,7 @@ ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
+FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
@@ -64,6 +65,7 @@ DiffEqBase = "6.144"
 EnumX = "1"
 Enzyme = "0.11.11"
 FastBroadcast = "0.2.8"
+FastClosures = "0.3"
 FastLevenbergMarquardt = "0.1"
 FiniteDiff = "2.21"
 FixedPointAcceleration = "0.3"
@@ -85,16 +87,17 @@ Printf = "1.9"
 Random = "1.91"
 RecursiveArrayTools = "3.2"
 Reexport = "1.2"
+SIAMFANLEquations = "1.0.1"
 SafeTestsets = "0.1"
 SciMLBase = "2.11"
 SciMLOperators = "0.3.7"
-SIAMFANLEquations = "1.0.1"
 SimpleNonlinearSolve = "1.0.2"
 SparseArrays = "1.9"
 SparseDiffTools = "2.14"
 SpeedMapping = "0.3"
 StableRNGs = "1"
 StaticArrays = "1.7"
+Sundials = "4.23.1"
 Symbolics = "5.13"
 Test = "1"
 UnPack = "1.0"
@@ -120,15 +123,16 @@ NonlinearProblemLibrary = "b7050fa9-e91f-4b37-bcee-a89a063da141"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
-SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SIAMFANLEquations = "084e46ad-d928-497d-ad5e-07fa361a48c4"
+SafeTestsets = "1bc83da4-3b8d-516f-aca4-4fe02f6d838f"
 SparseDiffTools = "47a9eef4-7e08-11e9-0b38-333d64bd3804"
 SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "StableRNGs", "MINPACK", "NLsolve", "OrdinaryDiffEq", "SpeedMapping", "FixedPointAcceleration", "SIAMFANLEquations"]
+test = ["Aqua", "Enzyme", "BenchmarkTools", "SafeTestsets", "Pkg", "Test", "ForwardDiff", "StaticArrays", "Symbolics", "LinearSolve", "Random", "LinearAlgebra", "Zygote", "SparseDiffTools", "NonlinearProblemLibrary", "LeastSquaresOptim", "FastLevenbergMarquardt", "NaNMath", "BandedMatrices", "DiffEqBase", "StableRNGs", "MINPACK", "NLsolve", "OrdinaryDiffEq", "SpeedMapping", "FixedPointAcceleration", "SIAMFANLEquations", "Sundials"]
diff --git a/docs/make.jl b/docs/make.jl
index e096f58d7..f494f711c 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -14,7 +14,7 @@ makedocs(; sitename = "NonlinearSolve.jl",
         DiffEqBase, SciMLBase],
     clean = true, doctest = false, linkcheck = true,
     linkcheck_ignore = ["https://twitter.com/ChrisRackauckas/status/1544743542094020615"],
-    warnonly = [:cross_references], checkdocs = :export,
+    checkdocs = :export,
     format = Documenter.HTML(assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/NonlinearSolve/stable/"),
     pages)
diff --git a/docs/src/api/siamfanlequations.md b/docs/src/api/siamfanlequations.md
index 2848a18ed..5ee36ba12 100644
--- a/docs/src/api/siamfanlequations.md
+++ b/docs/src/api/siamfanlequations.md
@@ -1,6 +1,7 @@
 # SIAMFANLEquations.jl
 
-This is an extension for importing solvers from [SIAMFANLEquations.jl](https://github.com/ctkelley/SIAMFANLEquations.jl) into the SciML
+This is an extension for importing solvers from
+[SIAMFANLEquations.jl](https://github.com/ctkelley/SIAMFANLEquations.jl) into the SciML
 interface. Note that these solvers do not come by default, and thus one needs to install
 the package before using these solvers:
 
diff --git a/docs/src/basics/solve.md b/docs/src/basics/solve.md
index 0f205aba1..cf78e1212 100644
--- a/docs/src/basics/solve.md
+++ b/docs/src/basics/solve.md
@@ -14,8 +14,8 @@ solve(prob::SciMLBase.NonlinearProblem, args...; kwargs...)
 ## Iteration Controls
 
   - `maxiters::Int`: The maximum number of iterations to perform. Defaults to `1000`.
-  - `abstol::Number`: The absolute tolerance.
-  - `reltol::Number`: The relative tolerance.
+  - `abstol::Number`: The absolute tolerance. Defaults to `real(oneunit(T)) * (eps(real(one(T))))^(4 // 5)`.
+  - `reltol::Number`: The relative tolerance. Defaults to `real(oneunit(T)) * (eps(real(one(T))))^(4 // 5)`.
   - `termination_condition`: Termination Condition from DiffEqBase. Defaults to
     `AbsSafeBestTerminationMode()` for `NonlinearSolve.jl` and `AbsTerminateMode()` for
     `SimpleNonlinearSolve.jl`.
diff --git a/docs/src/solvers/FixedPointSolvers.md b/docs/src/solvers/FixedPointSolvers.md
index 4eb3fdf0e..ed507098c 100644
--- a/docs/src/solvers/FixedPointSolvers.md
+++ b/docs/src/solvers/FixedPointSolvers.md
@@ -40,3 +40,10 @@ We are only listing the methods that natively solve fixed point problems.
 
   - `FixedPointAccelerationJL()`: accelerates the convergence of a mapping to a fixed point
     by the Anderson acceleration algorithm and a few other methods.
+
+### NLsolve.jl
+
+In our tests, we have found the anderson method implemented here to NOT be the most
+robust.
+
+  - `NLsolveJL(; method = :anderson)`: Anderson acceleration for fixed point problems.
diff --git a/docs/src/solvers/NonlinearSystemSolvers.md b/docs/src/solvers/NonlinearSystemSolvers.md
index bba941d86..c15948814 100644
--- a/docs/src/solvers/NonlinearSystemSolvers.md
+++ b/docs/src/solvers/NonlinearSystemSolvers.md
@@ -143,3 +143,10 @@ Newton-Krylov form. However, KINSOL is known to be less stable than some other
 implementations, as it has no line search or globalizer (trust region).
 
   - `KINSOL()`: The KINSOL method of the SUNDIALS C library
+
+### SIAMFANLEquations.jl
+
+SIAMFANLEquations.jl is a wrapper for the methods in the SIAMFANLEquations.jl library.
+
+  - `SIAMFANLEquationsJL()`: A wrapper for using the methods in
+    [SIAMFANLEquations.jl](https://github.com/ctkelley/SIAMFANLEquations.jl)
diff --git a/ext/NonlinearSolveFastLevenbergMarquardtExt.jl b/ext/NonlinearSolveFastLevenbergMarquardtExt.jl
index 6d33c23ba..fcda6e34d 100644
--- a/ext/NonlinearSolveFastLevenbergMarquardtExt.jl
+++ b/ext/NonlinearSolveFastLevenbergMarquardtExt.jl
@@ -5,7 +5,7 @@ import ConcreteStructs: @concrete
 import FastLevenbergMarquardt as FastLM
 import FiniteDiff, ForwardDiff
 
-function _fast_lm_solver(::FastLevenbergMarquardtJL{linsolve}, x) where {linsolve}
+@inline function _fast_lm_solver(::FastLevenbergMarquardtJL{linsolve}, x) where {linsolve}
     if linsolve === :cholesky
         return FastLM.CholeskySolver(ArrayInterface.undefmatrix(x))
     elseif linsolve === :qr
@@ -15,6 +15,7 @@ function _fast_lm_solver(::FastLevenbergMarquardtJL{linsolve}, x) where {linsolv
     end
 end
 
+# TODO: Implement reinit
 @concrete struct FastLevenbergMarquardtJLCache
     f!
     J!
@@ -25,68 +26,27 @@ end
     kwargs
 end
 
-@concrete struct InplaceFunction{iip} <: Function
-    f
-end
-
-(f::InplaceFunction{true})(fx, x, p) = f.f(fx, x, p)
-(f::InplaceFunction{false})(fx, x, p) = (fx .= f.f(x, p))
-
 function SciMLBase.__init(prob::NonlinearLeastSquaresProblem,
-        alg::FastLevenbergMarquardtJL, args...; alias_u0 = false, abstol = 1e-8,
-        reltol = 1e-8, maxiters = 1000, kwargs...)
+        alg::FastLevenbergMarquardtJL, args...; alias_u0 = false, abstol = nothing,
+        reltol = nothing, maxiters = 1000, kwargs...)
+    # FIXME: Support scalar u0
+    prob.u0 isa Number &&
+        throw(ArgumentError("FastLevenbergMarquardtJL does not support scalar `u0`"))
     iip = SciMLBase.isinplace(prob)
     u = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
     fu = NonlinearSolve.evaluate_f(prob, u)
 
-    f! = InplaceFunction{iip}(prob.f)
+    f! = NonlinearSolve.__make_inplace{iip}(prob.f, nothing)
+
+    abstol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u))
+    reltol = NonlinearSolve.DEFAULT_TOLERANCE(reltol, eltype(u))
 
     if prob.f.jac === nothing
-        use_forward_diff = if alg.autodiff === nothing
-            ForwardDiff.can_dual(eltype(u))
-        else
-            alg.autodiff isa AutoForwardDiff
-        end
-        uf = SciMLBase.JacobianWrapper{iip}(prob.f, prob.p)
-        if use_forward_diff
-            cache = iip ? ForwardDiff.JacobianConfig(uf, fu, u) :
-                    ForwardDiff.JacobianConfig(uf, u)
-        else
-            cache = FiniteDiff.JacobianCache(u, fu)
-        end
-        J! = if iip
-            if use_forward_diff
-                fu_cache = similar(fu)
-                function (J, x, p)
-                    uf.p = p
-                    ForwardDiff.jacobian!(J, uf, fu_cache, x, cache)
-                    return J
-                end
-            else
-                function (J, x, p)
-                    uf.p = p
-                    FiniteDiff.finite_difference_jacobian!(J, uf, x, cache)
-                    return J
-                end
-            end
-        else
-            if use_forward_diff
-                function (J, x, p)
-                    uf.p = p
-                    ForwardDiff.jacobian!(J, uf, x, cache)
-                    return J
-                end
-            else
-                function (J, x, p)
-                    uf.p = p
-                    J_ = FiniteDiff.finite_difference_jacobian(uf, x, cache)
-                    copyto!(J, J_)
-                    return J
-                end
-            end
-        end
+        alg = NonlinearSolve.get_concrete_algorithm(alg, prob)
+        J! = NonlinearSolve.__construct_jac(prob, alg, u;
+            can_handle_arbitrary_dims = Val(true))
     else
-        J! = InplaceFunction{iip}(prob.f.jac)
+        J! = NonlinearSolve.__make_inplace{iip}(prob.f.jac, nothing)
     end
 
     J = similar(u, length(fu), length(u))
@@ -95,17 +55,16 @@ function SciMLBase.__init(prob::NonlinearLeastSquaresProblem,
     LM = FastLM.LMWorkspace(u, fu, J)
 
     return FastLevenbergMarquardtJLCache(f!, J!, prob, alg, LM, solver,
-        (; xtol = abstol, ftol = reltol, maxit = maxiters, alg.factor, alg.factoraccept,
-            alg.factorreject, alg.minscale, alg.maxscale, alg.factorupdate, alg.minfactor,
-            alg.maxfactor, kwargs...))
+        (; xtol = reltol, ftol = reltol, gtol = abstol, maxit = maxiters, alg.factor,
+            alg.factoraccept, alg.factorreject, alg.minscale, alg.maxscale,
+            alg.factorupdate, alg.minfactor, alg.maxfactor))
 end
 
 function SciMLBase.solve!(cache::FastLevenbergMarquardtJLCache)
     res, fx, info, iter, nfev, njev, LM, solver = FastLM.lmsolve!(cache.f!, cache.J!,
         cache.lmworkspace, cache.prob.p; cache.solver, cache.kwargs...)
     stats = SciMLBase.NLStats(nfev, njev, -1, -1, iter)
-    retcode = info == 1 ? ReturnCode.Success :
-              (info == -1 ? ReturnCode.MaxIters : ReturnCode.Default)
+    retcode = info == -1 ? ReturnCode.MaxIters : ReturnCode.Success
     return SciMLBase.build_solution(cache.prob, cache.alg, res, fx;
         retcode, original = (res, fx, info, iter, nfev, njev, LM, solver), stats)
 end
diff --git a/ext/NonlinearSolveFixedPointAccelerationExt.jl b/ext/NonlinearSolveFixedPointAccelerationExt.jl
index e39946652..2c7ed376e 100644
--- a/ext/NonlinearSolveFixedPointAccelerationExt.jl
+++ b/ext/NonlinearSolveFixedPointAccelerationExt.jl
@@ -7,43 +7,27 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::FixedPointAccelerationJL
         show_trace::Val{PrintReports} = Val(false), termination_condition = nothing,
         kwargs...) where {PrintReports}
     @assert (termination_condition ===
-             nothing)||(termination_condition isa AbsNormTerminationMode) "SpeedMappingJL does not support termination conditions!"
+             nothing)||(termination_condition isa AbsNormTerminationMode) "FixedPointAccelerationJL does not support termination conditions!"
 
-    u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
-    u_size = size(u0)
-    T = eltype(u0)
-    iip = isinplace(prob)
-    p = prob.p
+    f, u0 = NonlinearSolve.__construct_f(prob; alias_u0, make_fixed_point = Val(true),
+        force_oop = Val(true))
 
-    if !iip && prob.u0 isa Number
-        # FixedPointAcceleration makes the scalar problem into a vector problem
-        f = (u) -> [prob.f(u[1], p) .+ u[1]]
-    elseif !iip && prob.u0 isa AbstractVector
-        f = (u) -> (prob.f(u, p) .+ u)
-    elseif !iip && prob.u0 isa AbstractArray
-        f = (u) -> vec(prob.f(reshape(u, u_size), p) .+ u)
-    elseif iip && prob.u0 isa AbstractVector
-        du = similar(u0)
-        f = (u) -> (prob.f(du, u, p); du .+ u)
-    else
-        du = similar(u0)
-        f = (u) -> (prob.f(du, reshape(u, u_size), p); vec(du) .+ u)
-    end
-
-    tol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol
+    tol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u0))
 
-    sol = fixed_point(f, NonlinearSolve._vec(u0); Algorithm = alg.algorithm,
+    sol = fixed_point(f, u0; Algorithm = alg.algorithm,
         ConvergenceMetricThreshold = tol, MaxIter = maxiters, MaxM = alg.m,
         ExtrapolationPeriod = alg.extrapolation_period, Dampening = alg.dampening,
         PrintReports, ReplaceInvalids = alg.replace_invalids,
         ConditionNumberThreshold = alg.condition_number_threshold, quiet_errors = true)
 
-    res = prob.u0 isa Number ? first(sol.FixedPoint_) : sol.FixedPoint_
-    if res === missing
+    if sol.FixedPoint_ === missing
+        u0 = prob.u0 isa Number ? u0[1] : u0
         resid = NonlinearSolve.evaluate_f(prob, u0)
         res = u0
         converged = false
     else
+        res = prob.u0 isa Number ? first(sol.FixedPoint_) :
+              reshape(sol.FixedPoint_, size(prob.u0))
         resid = NonlinearSolve.evaluate_f(prob, res)
         converged = maximum(abs, resid) ≤ tol
     end
diff --git a/ext/NonlinearSolveLeastSquaresOptimExt.jl b/ext/NonlinearSolveLeastSquaresOptimExt.jl
index 004c26670..e50469cec 100644
--- a/ext/NonlinearSolveLeastSquaresOptimExt.jl
+++ b/ext/NonlinearSolveLeastSquaresOptimExt.jl
@@ -4,7 +4,7 @@ using NonlinearSolve, SciMLBase
 import ConcreteStructs: @concrete
 import LeastSquaresOptim as LSO
 
-function _lso_solver(::LeastSquaresOptimJL{alg, linsolve}) where {alg, linsolve}
+@inline function _lso_solver(::LeastSquaresOptimJL{alg, linsolve}) where {alg, linsolve}
     ls = linsolve === :qr ? LSO.QR() :
          (linsolve === :cholesky ? LSO.Cholesky() :
           (linsolve === :lsmr ? LSO.LSMR() : nothing))
@@ -17,6 +17,7 @@ function _lso_solver(::LeastSquaresOptimJL{alg, linsolve}) where {alg, linsolve}
     end
 end
 
+# TODO: Implement reinit
 @concrete struct LeastSquaresOptimJLCache
     prob
     alg
@@ -24,34 +25,29 @@ end
     kwargs
 end
 
-@concrete struct FunctionWrapper{iip}
-    f
-    p
-end
-
-(f::FunctionWrapper{true})(du, u) = f.f(du, u, f.p)
-(f::FunctionWrapper{false})(du, u) = (du .= f.f(u, f.p))
-
 function SciMLBase.__init(prob::NonlinearLeastSquaresProblem, alg::LeastSquaresOptimJL,
-        args...; alias_u0 = false, abstol = 1e-8, reltol = 1e-8, verbose = false,
-        maxiters = 1000, kwargs...)
+        args...; alias_u0 = false, abstol = nothing, show_trace::Val{ShT} = Val(false),
+        trace_level = TraceMinimal(), store_trace::Val{StT} = Val(false), maxiters = 1000,
+        reltol = nothing, kwargs...) where {ShT, StT}
     iip = SciMLBase.isinplace(prob)
-    u = alias_u0 ? prob.u0 : deepcopy(prob.u0)
+    u = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
+
+    abstol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u))
+    reltol = NonlinearSolve.DEFAULT_TOLERANCE(reltol, eltype(u))
 
-    f! = FunctionWrapper{iip}(prob.f, prob.p)
-    g! = prob.f.jac === nothing ? nothing : FunctionWrapper{iip}(prob.f.jac, prob.p)
+    f! = NonlinearSolve.__make_inplace{iip}(prob.f, prob.p)
+    g! = NonlinearSolve.__make_inplace{iip}(prob.f.jac, prob.p)
 
     resid_prototype = prob.f.resid_prototype === nothing ?
-                      (!iip ? prob.f(u, prob.p) : zeros(u)) :
-                      prob.f.resid_prototype
+                      (!iip ? prob.f(u, prob.p) : zeros(u)) : prob.f.resid_prototype
 
     lsoprob = LSO.LeastSquaresProblem(; x = u, f!, y = resid_prototype, g!,
         J = prob.f.jac_prototype, alg.autodiff, output_length = length(resid_prototype))
     allocated_prob = LSO.LeastSquaresProblemAllocated(lsoprob, _lso_solver(alg))
 
     return LeastSquaresOptimJLCache(prob, alg, allocated_prob,
-        (; x_tol = abstol, f_tol = reltol, iterations = maxiters, show_trace = verbose,
-            kwargs...))
+        (; x_tol = reltol, f_tol = abstol, g_tol = abstol, iterations = maxiters,
+            show_trace = ShT, store_trace = StT, show_every = trace_level.print_frequency))
 end
 
 function SciMLBase.solve!(cache::LeastSquaresOptimJLCache)
diff --git a/ext/NonlinearSolveMINPACKExt.jl b/ext/NonlinearSolveMINPACKExt.jl
index d7d73de6b..0d3b8fc42 100644
--- a/ext/NonlinearSolveMINPACKExt.jl
+++ b/ext/NonlinearSolveMINPACKExt.jl
@@ -2,89 +2,52 @@ module NonlinearSolveMINPACKExt
 
 using NonlinearSolve, DiffEqBase, SciMLBase
 using MINPACK
+import FastClosures: @closure
 
 function SciMLBase.__solve(prob::Union{NonlinearProblem{uType, iip},
             NonlinearLeastSquaresProblem{uType, iip}}, alg::CMINPACK, args...;
-        abstol = nothing, maxiters = 100000, alias_u0::Bool = false,
-        termination_condition = nothing, kwargs...) where {uType, iip}
+        abstol = nothing, maxiters = 1000, alias_u0::Bool = false,
+        show_trace::Val{ShT} = Val(false), store_trace::Val{StT} = Val(false),
+        termination_condition = nothing, kwargs...) where {uType, iip, ShT, StT}
     @assert (termination_condition ===
              nothing)||(termination_condition isa AbsNormTerminationMode) "CMINPACK does not support termination conditions!"
 
-    if prob.u0 isa Number
-        u0 = [prob.u0]
-    else
-        u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
-    end
-
-    T = eltype(u0)
-    sizeu = size(prob.u0)
-    p = prob.p
+    f!_, u0 = NonlinearSolve.__construct_f(prob; alias_u0)
+    f! = @closure (du, u) -> (f!_(du, u); Cint(0))
 
-    # unwrapping alg params
-    show_trace = alg.show_trace
-    tracing = alg.tracing
-
-    if !iip && prob.u0 isa Number
-        f! = (du, u) -> (du .= prob.f(first(u), p); Cint(0))
-    elseif !iip && prob.u0 isa AbstractVector
-        f! = (du, u) -> (du .= prob.f(u, p); Cint(0))
-    elseif !iip && prob.u0 isa AbstractArray
-        f! = (du, u) -> (du .= vec(prob.f(reshape(u, sizeu), p)); Cint(0))
-    elseif prob.u0 isa AbstractVector
-        f! = (du, u) -> prob.f(du, u, p)
-    else # Then it's an in-place function on an abstract array
-        f! = (du, u) -> (prob.f(reshape(du, sizeu), reshape(u, sizeu), p); du = vec(du); 0)
-    end
-
-    u = zero(u0)
-    resid = NonlinearSolve.evaluate_f(prob, u)
+    resid = NonlinearSolve.evaluate_f(prob, prob.u0)
     m = length(resid)
-    size_jac = (length(resid), length(u))
 
     method = ifelse(alg.method === :auto,
         ifelse(prob isa NonlinearLeastSquaresProblem, :lm, :hybr), alg.method)
 
-    abstol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol
+    show_trace = alg.show_trace || ShT
+    tracing = alg.tracing || StT
+    tol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u0))
+
+    jac!_ = NonlinearSolve.__construct_jac(prob, alg, u0)
 
-    if SciMLBase.has_jac(prob.f)
-        if !iip && prob.u0 isa Number
-            g! = (du, u) -> (du .= prob.f.jac(first(u), p); Cint(0))
-        elseif !iip && prob.u0 isa AbstractVector
-            g! = (du, u) -> (du .= prob.f.jac(u, p); Cint(0))
-        elseif !iip && prob.u0 isa AbstractArray
-            g! = (du, u) -> (du .= vec(prob.f.jac(reshape(u, sizeu), p)); Cint(0))
-        elseif prob.u0 isa AbstractVector
-            g! = (du, u) -> prob.f.jac(du, u, p)
-        else # Then it's an in-place function on an abstract array
-            g! = function (du, u)
-                prob.f.jac(reshape(du, size_jac), reshape(u, sizeu), p)
-                return Cint(0)
-            end
-        end
-        original = MINPACK.fsolve(f!, g!, vec(u0), m; tol = abstol, show_trace, tracing,
-            method, iterations = maxiters, kwargs...)
+    if jac!_ === nothing
+        original = MINPACK.fsolve(f!, u0, m; tol, show_trace, tracing, method,
+            iterations = maxiters)
     else
-        original = MINPACK.fsolve(f!, vec(u0), m; tol = abstol, show_trace, tracing,
-            method, iterations = maxiters, kwargs...)
+        jac! = @closure((J, u)->(jac!_(J, u); Cint(0)))
+        original = MINPACK.fsolve(f!, jac!, u0, m; tol, show_trace, tracing, method,
+            iterations = maxiters)
     end
 
-    u = reshape(original.x, size(u))
-    resid = original.f
-    # retcode = original.converged ? ReturnCode.Success : ReturnCode.Failure
-    # MINPACK lies about convergence? or maybe uses some other criteria?
-    # We just check for absolute tolerance on the residual
-    objective = maximum(abs, resid)
-    retcode = ifelse(objective ≤ abstol, ReturnCode.Success, ReturnCode.Failure)
+    u = original.x
+    resid_ = original.f
+    objective = maximum(abs, resid_)
+    retcode = ifelse(objective ≤ tol, ReturnCode.Success, ReturnCode.Failure)
 
-    # These are only meaningful if `tracing = true`
+    # These are only meaningful if `store_trace = Val(true)`
     stats = SciMLBase.NLStats(original.trace.f_calls, original.trace.g_calls,
         original.trace.g_calls, original.trace.g_calls, -1)
 
-    if prob.u0 isa Number
-        return SciMLBase.build_solution(prob, alg, u[1], resid[1]; stats, retcode, original)
-    else
-        return SciMLBase.build_solution(prob, alg, u, resid; stats, retcode, original)
-    end
+    u_ = prob.u0 isa Number ? original.x[1] : reshape(original.x, size(prob.u0))
+    resid_ = prob.u0 isa Number ? resid_[1] : reshape(resid_, size(resid))
+    return SciMLBase.build_solution(prob, alg, u_, resid_; retcode, original, stats)
 end
 
 end
diff --git a/ext/NonlinearSolveNLsolveExt.jl b/ext/NonlinearSolveNLsolveExt.jl
index d1a3e557f..7d1eff02d 100644
--- a/ext/NonlinearSolveNLsolveExt.jl
+++ b/ext/NonlinearSolveNLsolveExt.jl
@@ -1,7 +1,6 @@
 module NonlinearSolveNLsolveExt
 
 using NonlinearSolve, NLsolve, DiffEqBase, SciMLBase
-import UnPack: @unpack
 
 function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...;
         abstol = nothing, maxiters = 1000, alias_u0::Bool = false,
@@ -9,84 +8,47 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::NLsolveJL, args...;
     @assert (termination_condition ===
              nothing)||(termination_condition isa AbsNormTerminationMode) "NLsolveJL does not support termination conditions!"
 
-    if prob.u0 isa Number
-        u0 = [prob.u0]
-    else
-        u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
-    end
-
-    T = eltype(u0)
-    iip = isinplace(prob)
-
-    sizeu = size(prob.u0)
-    p = prob.p
+    f!, u0 = NonlinearSolve.__construct_f(prob; alias_u0)
 
     # unwrapping alg params
-    @unpack method, autodiff, store_trace, extended_trace, linesearch, linsolve = alg
-    @unpack factor, autoscale, m, beta, show_trace = alg
-
-    if !iip && prob.u0 isa Number
-        f! = (du, u) -> (du .= prob.f(first(u), p); Cint(0))
-    elseif !iip && prob.u0 isa AbstractVector
-        f! = (du, u) -> (du .= prob.f(u, p); Cint(0))
-    elseif !iip && prob.u0 isa AbstractArray
-        f! = (du, u) -> (du .= vec(prob.f(reshape(u, sizeu), p)); Cint(0))
-    elseif prob.u0 isa AbstractVector
-        f! = (du, u) -> prob.f(du, u, p)
-    else # Then it's an in-place function on an abstract array
-        f! = (du, u) -> (prob.f(reshape(du, sizeu), reshape(u, sizeu), p); du = vec(du); 0)
-    end
+    (; method, autodiff, store_trace, extended_trace, linesearch, linsolve, factor,
+    autoscale, m, beta, show_trace) = alg
 
     if prob.u0 isa Number
         resid = [NonlinearSolve.evaluate_f(prob, first(u0))]
     else
-        resid = NonlinearSolve.evaluate_f(prob, u0)
+        resid = NonlinearSolve.evaluate_f(prob, prob.u0)
     end
 
-    size_jac = (length(resid), length(u0))
+    jac! = NonlinearSolve.__construct_jac(prob, alg, u0)
 
-    if SciMLBase.has_jac(prob.f)
-        if !iip && prob.u0 isa Number
-            g! = (du, u) -> (du .= prob.f.jac(first(u), p); Cint(0))
-        elseif !iip && prob.u0 isa AbstractVector
-            g! = (du, u) -> (du .= prob.f.jac(u, p); Cint(0))
-        elseif !iip && prob.u0 isa AbstractArray
-            g! = (du, u) -> (du .= vec(prob.f.jac(reshape(u, sizeu), p)); Cint(0))
-        elseif prob.u0 isa AbstractVector
-            g! = (du, u) -> prob.f.jac(du, u, p)
-        else # Then it's an in-place function on an abstract array
-            g! = function (du, u)
-                prob.f.jac(reshape(du, size_jac), reshape(u, sizeu), p)
-                return Cint(0)
-            end
-        end
+    if jac! === nothing
+        df = OnceDifferentiable(f!, vec(u0), vec(resid); autodiff)
+    else
         if prob.f.jac_prototype !== nothing
             J = zero(prob.f.jac_prototype)
-            df = OnceDifferentiable(f!, g!, vec(u0), vec(resid), J)
+            df = OnceDifferentiable(f!, jac!, vec(u0), vec(resid), J)
         else
-            df = OnceDifferentiable(f!, g!, vec(u0), vec(resid))
+            df = OnceDifferentiable(f!, jac!, vec(u0), vec(resid))
         end
-    else
-        df = OnceDifferentiable(f!, vec(u0), vec(resid); autodiff)
     end
 
-    abstol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol
+    abstol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u0))
 
     original = nlsolve(df, vec(u0); ftol = abstol, iterations = maxiters, method,
         store_trace, extended_trace, linesearch, linsolve, factor, autoscale, m, beta,
         show_trace)
 
-    u = reshape(original.zero, size(u0))
-    f!(vec(resid), vec(u))
+    f!(vec(resid), original.zero)
+    u = prob.u0 isa Number ? original.zero[1] : reshape(original.zero, size(prob.u0))
+    resid = prob.u0 isa Number ? resid[1] : resid
+
     retcode = original.x_converged || original.f_converged ? ReturnCode.Success :
               ReturnCode.Failure
     stats = SciMLBase.NLStats(original.f_calls, original.g_calls, original.g_calls,
         original.g_calls, original.iterations)
-    if prob.u0 isa Number
-        return SciMLBase.build_solution(prob, alg, u[1], resid[1]; retcode, original, stats)
-    else
-        return SciMLBase.build_solution(prob, alg, u, resid; retcode, original, stats)
-    end
+
+    return SciMLBase.build_solution(prob, alg, u, resid; retcode, original, stats)
 end
 
 end
diff --git a/ext/NonlinearSolveSIAMFANLEquationsExt.jl b/ext/NonlinearSolveSIAMFANLEquationsExt.jl
index 47ecc96c9..a5bbd5e1b 100644
--- a/ext/NonlinearSolveSIAMFANLEquationsExt.jl
+++ b/ext/NonlinearSolveSIAMFANLEquationsExt.jl
@@ -2,136 +2,125 @@ module NonlinearSolveSIAMFANLEquationsExt
 
 using NonlinearSolve, SciMLBase
 using SIAMFANLEquations
-import ConcreteStructs: @concrete
-import UnPack: @unpack
 
-function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, args...; abstol = nothing,
-        reltol = nothing, alias_u0::Bool = false, maxiters = 1000, termination_condition = nothing, kwargs...)
-    @assert (termination_condition === nothing) || (termination_condition isa AbsNormTerminationMode) "SIAMFANLEquationsJL does not support termination conditions!"
+@inline function __siam_fanl_equations_retcode_mapping(sol)
+    if sol.errcode == 0
+        return ReturnCode.Success
+    elseif sol.errcode == 10
+        return ReturnCode.MaxIters
+    elseif sol.errcode == 1
+        return ReturnCode.Failure
+    elseif sol.errcode == -1
+        return ReturnCode.Default
+    end
+end
 
-    @unpack method, show_trace, delta, linsolve = alg
+@inline function __zeros_like(x, args...)
+    z = similar(x, args...)
+    fill!(z, zero(eltype(x)))
+    return z
+end
+
+# pseudo transient continuation has a fixed cost per iteration, iteration statistics are
+# not interesting here.
+@inline function __siam_fanl_equations_stats_mapping(method, sol)
+    method === :pseudotransient && return nothing
+    return SciMLBase.NLStats(sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0,
+        sum(sol.stats.iarm))
+end
+
+function SciMLBase.__solve(prob::NonlinearProblem, alg::SIAMFANLEquationsJL, args...;
+        abstol = nothing, reltol = nothing, alias_u0::Bool = false, maxiters = 1000,
+        termination_condition = nothing, show_trace::Val{ShT} = Val(false),
+        kwargs...) where {ShT}
+    @assert (termination_condition ===
+             nothing)||(termination_condition isa AbsNormTerminationMode) "SIAMFANLEquationsJL does not support termination conditions!"
+
+    (; method, delta, linsolve) = alg
 
     iip = SciMLBase.isinplace(prob)
-    T = eltype(prob.u0)
 
-    atol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol
-    rtol = reltol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : reltol
+    T = eltype(prob.u0)
+    atol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, T)
+    rtol = NonlinearSolve.DEFAULT_TOLERANCE(reltol, T)
 
     if prob.u0 isa Number
-        f! = if iip
-            function (u)
-                du = similar(u)
-                prob.f(du, u, prob.p)
-                return du
-            end
-        else
-            u -> prob.f(u, prob.p)
-        end
+        f = (u) -> prob.f(u, prob.p)
 
         if method == :newton
-            sol = nsolsc(f!, prob.u0; maxit = maxiters, atol = atol, rtol = rtol, printerr = show_trace)
+            sol = nsolsc(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT)
         elseif method == :pseudotransient
-            sol = ptcsolsc(f!, prob.u0; delta0 = delta, maxit = maxiters, atol = atol, rtol=rtol, printerr = show_trace)
+            sol = ptcsolsc(f, prob.u0; delta0 = delta, maxit = maxiters, atol, rtol,
+                printerr = ShT)
         elseif method == :secant
-            sol = secant(f!, prob.u0; maxit = maxiters, atol = atol, rtol = rtol, printerr = show_trace)
+            sol = secant(f, prob.u0; maxit = maxiters, atol, rtol, printerr = ShT)
         end
 
-        if sol.errcode == 0
-            retcode = ReturnCode.Success
-        elseif sol.errcode == 10
-            retcode = ReturnCode.MaxIters
-        elseif sol.errcode == 1
-            retcode = ReturnCode.Failure
-        elseif sol.errcode == -1
-            retcode = ReturnCode.Default
-        end
-        stats = method == :pseudotransient ? nothing : (SciMLBase.NLStats(sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0, sum(sol.stats.iarm)))
-        return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode, stats, original = sol)
-    else
-        u = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
+        retcode = __siam_fanl_equations_retcode_mapping(sol)
+        stats = __siam_fanl_equations_stats_mapping(method, sol)
+        resid = NonlinearSolve.evaluate_f(prob, sol.solution)
+        return SciMLBase.build_solution(prob, alg, sol.solution, resid; retcode,
+            stats, original = sol)
     end
 
-    if iip
-        f! = function (du, u)
-            prob.f(du, u, prob.p)
-            return du
-        end
-    else
-        f! = function (du, u)
-            du .= prob.f(u, prob.p)
-            return du
-        end
-    end
+    f!, u = NonlinearSolve.__construct_f(prob; alias_u0,
+        can_handle_arbitrary_dims = Val(true))
 
     # Allocate ahead for function
     N = length(u)
-    FS = zeros(T, N)
+    FS = __zeros_like(u, N)
 
     # Jacobian free Newton Krylov
     if linsolve !== nothing
         # Allocate ahead for Krylov basis
-        JVS = linsolve == :gmres ? zeros(T, N, 3) : zeros(T, N)
-        # `linsolve` as a Symbol to keep unified interface with other EXTs, SIAMFANLEquations directly use String to choose between different linear solvers
+        JVS = linsolve == :gmres ? __zeros_like(u, N, 3) : __zeros_like(u, N)
+        # `linsolve` as a Symbol to keep unified interface with other EXTs,
+        # SIAMFANLEquations directly use String to choose between different linear solvers
         linsolve_alg = String(linsolve)
 
         if method == :newton
-            sol = nsoli(f!, u, FS, JVS; lsolver = linsolve_alg, maxit = maxiters, atol = atol, rtol = rtol, printerr = show_trace)
+            sol = nsoli(f!, u, FS, JVS; lsolver = linsolve_alg, maxit = maxiters, atol,
+                rtol, printerr = ShT)
         elseif method == :pseudotransient
-            sol = ptcsoli(f!, u, FS, JVS; lsolver = linsolve_alg, maxit = maxiters, atol = atol, rtol = rtol, printerr = show_trace)
-        end
-        
-        if sol.errcode == 0
-            retcode = ReturnCode.Success
-        elseif sol.errcode == 10
-            retcode = ReturnCode.MaxIters
-        elseif sol.errcode == 1
-            retcode = ReturnCode.Failure
-        elseif sol.errcode == -1
-            retcode = ReturnCode.Default
+            sol = ptcsoli(f!, u, FS, JVS; lsolver = linsolve_alg, maxit = maxiters, atol,
+                rtol, printerr = ShT)
         end
-        stats = method == :pseudotransient ? nothing : (SciMLBase.NLStats(sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0, sum(sol.stats.iarm)))
-        return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode, stats, original = sol)
+
+        retcode = __siam_fanl_equations_retcode_mapping(sol)
+        stats = __siam_fanl_equations_stats_mapping(method, sol)
+        resid = NonlinearSolve.evaluate_f(prob, sol.solution)
+        return SciMLBase.build_solution(prob, alg, sol.solution, resid; retcode,
+            stats, original = sol)
     end
 
     # Allocate ahead for Jacobian
-    FPS = zeros(T, N, N)
+    FPS = __zeros_like(u, N, N)
+
     if prob.f.jac === nothing
         # Use the built-in Jacobian machinery
         if method == :newton
-            sol = nsol(f!, u, FS, FPS;
-                        sham=1, atol = atol, rtol = rtol, maxit = maxiters,
-                        printerr = show_trace)
+            sol = nsol(f!, u, FS, FPS; sham = 1, atol, rtol, maxit = maxiters,
+                printerr = ShT)
         elseif method == :pseudotransient
-            sol = ptcsol(f!, u, FS, FPS;
-                        atol = atol, rtol = rtol, maxit = maxiters,
-                        delta0 = delta, printerr = show_trace)
+            sol = ptcsol(f!, u, FS, FPS; atol, rtol, maxit = maxiters,
+                delta0 = delta, printerr = ShT)
         end
     else
         AJ!(J, u, x) = prob.f.jac(J, x, prob.p)
         if method == :newton
-            sol = nsol(f!, u, FS, FPS, AJ!;
-                        sham=1, atol = atol, rtol = rtol, maxit = maxiters,
-                        printerr = show_trace)
+            sol = nsol(f!, u, FS, FPS, AJ!; sham = 1, atol, rtol, maxit = maxiters,
+                printerr = ShT)
         elseif method == :pseudotransient
-            sol = ptcsol(f!, u, FS, FPS, AJ!;
-                        atol = atol, rtol = rtol, maxit = maxiters,
-                        delta0 = delta, printerr = show_trace)
+            sol = ptcsol(f!, u, FS, FPS, AJ!; atol, rtol, maxit = maxiters,
+                delta0 = delta, printerr = ShT)
         end
     end
 
-    if sol.errcode == 0
-        retcode = ReturnCode.Success
-    elseif sol.errcode == 10
-        retcode = ReturnCode.MaxIters
-    elseif sol.errcode == 1
-        retcode = ReturnCode.Failure
-    elseif sol.errcode == -1
-        retcode = ReturnCode.Default
-    end
-
-    # pseudo transient continuation has a fixed cost per iteration, iteration statistics are not interesting here.
-    stats = method == :pseudotransient ? nothing : (SciMLBase.NLStats(sum(sol.stats.ifun), sum(sol.stats.ijac), 0, 0, sum(sol.stats.iarm)))
-    return SciMLBase.build_solution(prob, alg, sol.solution, sol.history; retcode, stats, original = sol)
+    retcode = __siam_fanl_equations_retcode_mapping(sol)
+    stats = __siam_fanl_equations_stats_mapping(method, sol)
+    resid = NonlinearSolve.evaluate_f(prob, sol.solution)
+    return SciMLBase.build_solution(prob, alg, sol.solution, resid; retcode, stats,
+        original = sol)
 end
 
-end
\ No newline at end of file
+end
diff --git a/ext/NonlinearSolveSpeedMappingExt.jl b/ext/NonlinearSolveSpeedMappingExt.jl
index be637b124..9f15ab97b 100644
--- a/ext/NonlinearSolveSpeedMappingExt.jl
+++ b/ext/NonlinearSolveSpeedMappingExt.jl
@@ -9,30 +9,15 @@ function SciMLBase.__solve(prob::NonlinearProblem, alg::SpeedMappingJL, args...;
     @assert (termination_condition ===
              nothing)||(termination_condition isa AbsNormTerminationMode) "SpeedMappingJL does not support termination conditions!"
 
-    if typeof(prob.u0) <: Number
-        u0 = [prob.u0]
-    else
-        u0 = NonlinearSolve.__maybe_unaliased(prob.u0, alias_u0)
-    end
+    m!, u0 = NonlinearSolve.__construct_f(prob; alias_u0, make_fixed_point = Val(true),
+        can_handle_arbitrary_dims = Val(true))
 
-    T = eltype(u0)
-    iip = isinplace(prob)
-    p = prob.p
-
-    if !iip && prob.u0 isa Number
-        m! = (du, u) -> (du .= prob.f(first(u), p) .+ first(u))
-    elseif !iip
-        m! = (du, u) -> (du .= prob.f(u, p) .+ u)
-    else
-        m! = (du, u) -> (prob.f(du, u, p); du .+= u)
-    end
-
-    tol = abstol === nothing ? real(oneunit(T)) * (eps(real(one(T))))^(4 // 5) : abstol
+    tol = NonlinearSolve.DEFAULT_TOLERANCE(abstol, eltype(u0))
 
     sol = speedmapping(u0; m!, tol, Lp = Inf, maps_limit = maxiters, alg.orders,
         alg.check_obj, store_info, alg.σ_min, alg.stabilize)
     res = prob.u0 isa Number ? first(sol.minimizer) : sol.minimizer
-    resid = NonlinearSolve.evaluate_f(prob, sol.minimizer)
+    resid = NonlinearSolve.evaluate_f(prob, res)
 
     return SciMLBase.build_solution(prob, alg, res, resid;
         retcode = sol.converged ? ReturnCode.Success : ReturnCode.Failure,
diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl
index bf77135c8..9b8786380 100644
--- a/src/NonlinearSolve.jl
+++ b/src/NonlinearSolve.jl
@@ -17,6 +17,7 @@ import PrecompileTools: @recompile_invalidations, @compile_workload, @setup_work
     import ConcreteStructs: @concrete
     import EnumX: @enumx
     import FastBroadcast: @..
+    import FastClosures: @closure
     import FiniteDiff
     import ForwardDiff
     import ForwardDiff: Dual
@@ -169,6 +170,7 @@ function SciMLBase.solve!(cache::AbstractNonlinearSolveCache)
 end
 
 include("utils.jl")
+include("function_wrappers.jl")
 include("trace.jl")
 include("extension_algs.jl")
 include("linesearch.jl")
@@ -237,7 +239,8 @@ export RadiusUpdateSchemes
 export NewtonRaphson, TrustRegion, LevenbergMarquardt, DFSane, GaussNewton, PseudoTransient,
     Broyden, Klement, LimitedMemoryBroyden
 export LeastSquaresOptimJL,
-    FastLevenbergMarquardtJL, CMINPACK, NLsolveJL, FixedPointAccelerationJL, SpeedMappingJL, SIAMFANLEquationsJL
+    FastLevenbergMarquardtJL, CMINPACK, NLsolveJL, FixedPointAccelerationJL, SpeedMappingJL,
+    SIAMFANLEquationsJL
 export NonlinearSolvePolyAlgorithm,
     RobustMultiNewton, FastShortcutNonlinearPolyalg, FastShortcutNLLSPolyalg
 
diff --git a/src/extension_algs.jl b/src/extension_algs.jl
index ea92454b0..195cced8d 100644
--- a/src/extension_algs.jl
+++ b/src/extension_algs.jl
@@ -59,7 +59,7 @@ for solving `NonlinearLeastSquaresProblem`.
     This algorithm is only available if `FastLevenbergMarquardt.jl` is installed.
 """
 @concrete struct FastLevenbergMarquardtJL{linsolve} <: AbstractNonlinearSolveAlgorithm
-    autodiff
+    ad
     factor
     factoraccept
     factorreject
@@ -70,6 +70,12 @@ for solving `NonlinearLeastSquaresProblem`.
     maxfactor
 end
 
+function set_ad(alg::FastLevenbergMarquardtJL{linsolve}, ad) where {linsolve}
+    return FastLevenbergMarquardtJL{linsolve}(ad, alg.factor, alg.factoraccept,
+        alg.factorreject, alg.factorupdate, alg.minscale, alg.maxscale, alg.minfactor,
+        alg.maxfactor)
+end
+
 function FastLevenbergMarquardtJL(linsolve::Symbol = :cholesky; factor = 1e-6,
         factoraccept = 13.0, factorreject = 3.0, factorupdate = :marquardt,
         minscale = 1e-12, maxscale = 1e16, minfactor = 1e-28, maxfactor = 1e32,
@@ -88,12 +94,10 @@ function FastLevenbergMarquardtJL(linsolve::Symbol = :cholesky; factor = 1e-6,
 end
 
 """
-    CMINPACK(; show_trace::Bool=false, tracing::Bool=false, method::Symbol=:auto)
+    CMINPACK(; method::Symbol = :auto)
 
 ### Keyword Arguments
 
-  - `show_trace`: whether to show the trace.
-  - `tracing`: who the hell knows what this does. If you find out, please open an issue/PR.
   - `method`: the choice of method for the solver.
 
 ### Method Choices
@@ -134,27 +138,42 @@ struct CMINPACK <: AbstractNonlinearSolveAlgorithm
     method::Symbol
 end
 
-function CMINPACK(; show_trace::Bool = false, tracing::Bool = false, method::Symbol = :auto)
+function CMINPACK(; show_trace = missing, tracing = missing, method::Symbol = :auto)
     if Base.get_extension(@__MODULE__, :NonlinearSolveMINPACKExt) === nothing
         error("CMINPACK requires MINPACK.jl to be loaded")
     end
 
+    if show_trace !== missing
+        Base.depwarn("`show_trace` for CMINPACK has been deprecated and will be removed \
+                      in v4. Use the `show_trace` keyword argument via the logging API \
+                      https://docs.sciml.ai/NonlinearSolve/stable/basics/Logging/ \
+                      instead.", :CMINPACK)
+    else
+        show_trace = false
+    end
+
+    if tracing !== missing
+        Base.depwarn("`tracing` for CMINPACK has been deprecated and will be removed \
+                      in v4. Use the `store_trace` keyword argument via the logging API \
+                      https://docs.sciml.ai/NonlinearSolve/stable/basics/Logging/ \
+                      instead.", :CMINPACK)
+    else
+        tracing = false
+    end
+
     return CMINPACK(show_trace, tracing, method)
 end
 
 """
-    NLsolveJL(; method=:trust_region, autodiff=:central, store_trace=false,
-        extended_trace=false, linesearch=LineSearches.Static(),
-        linsolve=(x, A, b) -> copyto!(x, A\\b), factor = one(Float64), autoscale=true,
-        m=10, beta=one(Float64), show_trace=false)
+    NLsolveJL(; method = :trust_region, autodiff = :central, linesearch = Static(),
+        linsolve = (x, A, b) -> copyto!(x, A\\b), factor = one(Float64), autoscale = true,
+        m = 10, beta = one(Float64))
 
 ### Keyword Arguments
 
   - `method`: the choice of method for solving the nonlinear system.
   - `autodiff`: the choice of method for generating the Jacobian. Defaults to `:central` or
     central differencing via FiniteDiff.jl. The other choices are `:forward`
-  - `show_trace`: should a trace of the optimization algorithm's state be shown on `STDOUT`?
-  - `extended_trace`: should additional algorithm internals be added to the state trace?
   - `linesearch`: the line search method to be used within the solver method. The choices
     are line search types from
     [LineSearches.jl](https://github.com/JuliaNLSolvers/LineSearches.jl).
@@ -168,7 +187,6 @@ end
     constants are close to 1. If convergence fails, though, you may consider lowering it.
   - `beta`: It is also known as DIIS or Pulay mixing, this method is based on the acceleration
     of the fixed-point iteration xₙ₊₁ = xₙ + beta*f(xₙ), where by default beta = 1.
-  - `store_trace``: should a trace of the optimization algorithm's state be stored?
 
 ### Submethod Choice
 
@@ -195,20 +213,47 @@ Choices for methods in `NLsolveJL`:
     show_trace::Bool
 end
 
-function NLsolveJL(; method = :trust_region, autodiff = :central, store_trace = false,
-        extended_trace = false, linesearch = LineSearches.Static(),
+function NLsolveJL(; method = :trust_region, autodiff = :central, store_trace = missing,
+        extended_trace = missing, linesearch = LineSearches.Static(),
         linsolve = (x, A, b) -> copyto!(x, A \ b), factor = 1.0, autoscale = true, m = 10,
-        beta = one(Float64), show_trace = false)
+        beta = one(Float64), show_trace = missing)
     if Base.get_extension(@__MODULE__, :NonlinearSolveNLsolveExt) === nothing
         error("NLsolveJL requires NLsolve.jl to be loaded")
     end
 
+    if show_trace !== missing
+        Base.depwarn("`show_trace` for NLsolveJL has been deprecated and will be removed \
+                      in v4. Use the `show_trace` keyword argument via the logging API \
+                      https://docs.sciml.ai/NonlinearSolve/stable/basics/Logging/ \
+                      instead.", :NLsolveJL)
+    else
+        show_trace = false
+    end
+
+    if store_trace !== missing
+        Base.depwarn("`store_trace` for NLsolveJL has been deprecated and will be removed \
+                      in v4. Use the `store_trace` keyword argument via the logging API \
+                      https://docs.sciml.ai/NonlinearSolve/stable/basics/Logging/ \
+                      instead.", :NLsolveJL)
+    else
+        store_trace = false
+    end
+
+    if extended_trace !== missing
+        Base.depwarn("`extended_trace` for NLsolveJL has been deprecated and will be \
+                      removed in v4. Use the `trace_level = TraceAll()` keyword argument \
+                      via the logging API \
+                      https://docs.sciml.ai/NonlinearSolve/stable/basics/Logging/ instead.",
+            :NLsolveJL)
+    else
+        extended_trace = false
+    end
+
     return NLsolveJL(method, autodiff, store_trace, extended_trace, linesearch, linsolve,
         factor, autoscale, m, beta, show_trace)
 end
 
 """
-
     SpeedMappingJL(; σ_min = 0.0, stabilize::Bool = false, check_obj::Bool = false,
         orders::Vector{Int} = [3, 3, 2], time_limit::Real = 1000)
 
@@ -325,13 +370,11 @@ function FixedPointAccelerationJL(; algorithm = :Anderson, m = missing,
 end
 
 """
-
-    SIAMFANLEquationsJL(; method = :newton, autodiff = :central, show_trace = false, delta = 1e-3, linsolve = nothing)
+    SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing)
 
 ### Keyword Arguments
 
   - `method`: the choice of method for solving the nonlinear system.
-  - `show_trace`: whether to show the trace.
   - `delta`: initial pseudo time step, default is 1e-3.
   - `linsolve` : JFNK linear solvers, choices are `gmres` and `bicgstab`.
 
@@ -341,16 +384,16 @@ end
   - `:pseudotransient`: Pseudo transient method.
   - `:secant`: Secant method for scalar equations.
 """
-@concrete struct SIAMFANLEquationsJL <: AbstractNonlinearAlgorithm
+@concrete struct SIAMFANLEquationsJL{L <: Union{Symbol, Nothing}} <:
+                 AbstractNonlinearSolveAlgorithm
     method::Symbol
-    show_trace::Bool
     delta
-    linsolve::Union{Symbol, Nothing}
+    linsolve::L
 end
 
-function SIAMFANLEquationsJL(; method = :newton, show_trace = false, delta = 1e-3, linsolve = nothing)
+function SIAMFANLEquationsJL(; method = :newton, delta = 1e-3, linsolve = nothing)
     if Base.get_extension(@__MODULE__, :NonlinearSolveSIAMFANLEquationsExt) === nothing
         error("SIAMFANLEquationsJL requires SIAMFANLEquations.jl to be loaded")
     end
-  return SIAMFANLEquationsJL(method, show_trace, delta, linsolve)
+    return SIAMFANLEquationsJL(method, delta, linsolve)
 end
diff --git a/src/function_wrappers.jl b/src/function_wrappers.jl
new file mode 100644
index 000000000..599127f39
--- /dev/null
+++ b/src/function_wrappers.jl
@@ -0,0 +1,188 @@
+# NonlinearSolve can handle all NonlinearFunction specifications but that is not true for
+# downstream packages. Make conversion to those easier.
+function __construct_f(prob; alias_u0::Bool = false, can_handle_oop::Val{OOP} = Val(false),
+        can_handle_scalar::Val{SCALAR} = Val(false), make_fixed_point::Val{FP} = Val(false),
+        can_handle_arbitrary_dims::Val{DIMS} = Val(false),
+        force_oop::Val{FOOP} = Val(false)) where {SCALAR, OOP, DIMS, FP, FOOP}
+    if !OOP && SCALAR
+        error("Incorrect Specification: OOP not supported but scalar supported.")
+    end
+
+    resid = evaluate_f(prob, prob.u0)
+
+    if SCALAR || !(prob.u0 isa Number)
+        u0 = __maybe_unaliased(prob.u0, alias_u0)
+    else
+        u0 = [prob.u0]
+    end
+
+    f = if FP
+        if isinplace(prob)
+            @closure (du, u, p) -> begin
+                prob.f(du, u, p)
+                @. du += u
+            end
+        else
+            @closure (u, p) -> prob.f(u, p) .+ u
+        end
+    else
+        prob.f
+    end
+
+    ff = if isinplace(prob)
+        ninputs = 2
+        if DIMS || u0 isa AbstractVector
+            @closure (du, u) -> (f(du, u, prob.p); du)
+        else
+            u0_size = size(u0)
+            du_size = size(resid)
+            @closure (du, u) -> (f(reshape(du, du_size), reshape(u, u0_size), prob.p); du)
+        end
+    else
+        if prob.u0 isa Number
+            if SCALAR
+                ninputs = 1
+                @closure (u) -> f(u, prob.p)
+            elseif OOP
+                ninputs = 1
+                @closure (u) -> [f(first(u), prob.p)]
+            else
+                ninputs = 2
+                resid = [resid]
+                @closure (du, u) -> (du[1] = f(first(u), prob.p); du)
+            end
+        else
+            if OOP
+                ninputs = 1
+                if DIMS
+                    @closure (u) -> f(u, prob.p)
+                else
+                    u0_size = size(u0)
+                    @closure (u) -> _vec(f(reshape(u, u0_size), prob.p))
+                end
+            else
+                ninputs = 2
+                if DIMS
+                    @closure (du, u) -> (copyto!(du, f(u, prob.p)); du)
+                else
+                    u0_size = size(u0)
+                    @closure (du, u) -> begin
+                        copyto!(vec(du), vec(f(reshape(u, u0_size), prob.p)))
+                        return du
+                    end
+                end
+            end
+        end
+    end
+
+    f_final = if FOOP
+        if ninputs == 1
+            ff
+        else
+            du_ = DIMS ? similar(resid) : _vec(similar(resid))
+            @closure (u) -> (ff(du_, u); du_)
+        end
+    else
+        ff
+    end
+
+    return f_final, ifelse(DIMS, u0, _vec(u0))
+end
+
+function __construct_jac(prob, alg, u0; can_handle_oop::Val{OOP} = Val(false),
+        can_handle_scalar::Val{SCALAR} = Val(false),
+        can_handle_arbitrary_dims::Val{DIMS} = Val(false)) where {SCALAR, OOP, DIMS}
+    if SciMLBase.has_jac(prob.f)
+        jac = prob.f.jac
+
+        jac_final = if isinplace(prob)
+            if DIMS || u0 isa AbstractVector
+                @closure (J, u) -> (jac(reshape(J, :, length(u)), u, prob.p); J)
+            else
+                u0_size = size(u0)
+                @closure (J, u) -> (jac(reshape(J, :, length(u)), reshape(u, u0_size),
+                    prob.p);
+                J)
+            end
+        else
+            if prob.u0 isa Number
+                if SCALAR
+                    @closure (u) -> jac(u, prob.p)
+                elseif OOP
+                    @closure (u) -> [jac(first(u), prob.p)]
+                else
+                    @closure (J, u) -> (J[1] = jac(first(u), prob.p); J)
+                end
+            else
+                if OOP
+                    if DIMS
+                        @closure (u) -> jac(u, prob.p)
+                    else
+                        u0_size = size(u0)
+                        @closure (u) -> jac(reshape(u, u0_size), prob.p)
+                    end
+                else
+                    if DIMS
+                        @closure (J, u) -> (copyto!(J, jac(u, prob.p)); J)
+                    else
+                        u0_size = size(u0)
+                        @closure (J, u) -> begin
+                            copyto!(J, jac(reshape(u, u0_size), prob.p))
+                            return J
+                        end
+                    end
+                end
+            end
+        end
+
+        return jac_final
+    end
+
+    hasfield(typeof(alg), :ad) || return nothing
+
+    uf, _, J, fu, jac_cache, _, _, _ = jacobian_caches(alg, prob.f, u0, prob.p,
+        Val{isinplace(prob)}(); lininit = Val(false), linsolve_with_JᵀJ = Val(false))
+    stats = SciMLBase.NLStats(0, 0, 0, 0, 0)
+    return JacobianFunctionCache{isinplace(prob)}(J, prob.f, uf, u0, prob.p, jac_cache,
+        alg, fu, stats)
+end
+
+# Currently used only in some of the extensions. Plan is to eventually use it in all the
+# native algorithms and other extensions to provide better jacobian support
+@concrete struct JacobianFunctionCache{iip, U, P} <: Function
+    J
+    f
+    uf
+    u::U
+    p::P
+    jac_cache
+    alg
+    fu_cache
+    stats
+end
+
+SciMLBase.isinplace(::JacobianFunctionCache{iip}) where {iip} = iip
+
+function (jac_cache::JacobianFunctionCache{iip, U, P})(J::AbstractMatrix, u::U,
+        p::P = jac_cache.p) where {iip, U, P}
+    jacobian!!(J, jac_cache; u, p)
+    return J
+end
+function (jac_cache::JacobianFunctionCache{iip, U, P})(u::U, p::P) where {iip, U, P}
+    return jacobian!!(cache.J, jac_cache; u, p)
+end
+
+@concrete struct InplaceFunction{iip} <: Function
+    f
+    p
+end
+
+(f::InplaceFunction{true})(du, u) = f.f(du, u, f.p)
+(f::InplaceFunction{true})(du, u, p) = f.f(du, u, p)
+(f::InplaceFunction{false})(du, u) = (du .= f.f(u, f.p))
+(f::InplaceFunction{false})(du, u, p) = (du .= f.f(u, p))
+
+struct __make_inplace{iip} end
+
+@inline __make_inplace{iip}(f::F, p) where {iip, F} = InplaceFunction{iip}(f, p)
+@inline __make_inplace{iip}(::Nothing, p) where {iip} = nothing
diff --git a/src/jacobian.jl b/src/jacobian.jl
index e52021d55..20825ebda 100644
--- a/src/jacobian.jl
+++ b/src/jacobian.jl
@@ -9,7 +9,7 @@ isinplace(JᵀJ::KrylovJᵀJ) = isinplace(JᵀJ.Jᵀ)
 
 # Select if we are going to use sparse differentiation or not
 sparsity_detection_alg(_, _) = NoSparsityDetection()
-function sparsity_detection_alg(f, ad::AbstractSparseADType)
+function sparsity_detection_alg(f::NonlinearFunction, ad::AbstractSparseADType)
     if f.sparsity === nothing
         if f.jac_prototype === nothing
             if is_extension_loaded(Val(:Symbolics))
@@ -47,11 +47,12 @@ function sparsity_detection_alg(f, ad::AbstractSparseADType)
 end
 
 # NoOp for Jacobian if it is not a Abstract Array -- For eg, JacVec Operator
-jacobian!!(J, _) = J
+jacobian!!(J, cache; u = nothing, p = nothing) = J
 # `!!` notation is from BangBang.jl since J might be jacobian in case of oop `f.jac`
 # and we don't want wasteful `copyto!`
-function jacobian!!(J::Union{AbstractMatrix{<:Number}, Nothing}, cache)
-    @unpack f, uf, u, p, jac_cache, alg, fu_cache = cache
+function jacobian!!(J::Union{AbstractMatrix{<:Number}, Nothing}, cache; u = cache.u,
+        p = cache.p)
+    @unpack f, uf, jac_cache, alg, fu_cache = cache
     cache.stats.njacs += 1
     iip = isinplace(cache)
     if iip
@@ -72,9 +73,9 @@ function jacobian!!(J::Union{AbstractMatrix{<:Number}, Nothing}, cache)
     end
 end
 # Scalar case
-function jacobian!!(::Number, cache)
+function jacobian!!(::Number, cache; u = cache.u, p = cache.p)
     cache.stats.njacs += 1
-    return last(value_derivative(cache.uf, cache.u))
+    return last(value_derivative(cache.uf, u))
 end
 
 # Build Jacobian Caches
diff --git a/src/utils.jl b/src/utils.jl
index 094ddb040..9bf4f6987 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -1,5 +1,7 @@
 const DEFAULT_NORM = DiffEqBase.NONLINEARSOLVE_DEFAULT_NORM
 
+@inline DEFAULT_TOLERANCE(args...) = DiffEqBase._get_tolerance(args...)
+
 @concrete mutable struct FakeLinearSolveJLCache
     A
     b
diff --git a/test/23_test_problems.jl b/test/core/23_test_problems.jl
similarity index 100%
rename from test/23_test_problems.jl
rename to test/core/23_test_problems.jl
diff --git a/test/forward_ad.jl b/test/core/forward_ad.jl
similarity index 87%
rename from test/forward_ad.jl
rename to test/core/forward_ad.jl
index 65d4e3c15..41fe5015e 100644
--- a/test/forward_ad.jl
+++ b/test/core/forward_ad.jl
@@ -1,4 +1,5 @@
-using ForwardDiff, NonlinearSolve, MINPACK, NLsolve, StaticArrays, Test, LinearAlgebra
+using ForwardDiff,
+    NonlinearSolve, MINPACK, NLsolve, StaticArrays, Sundials, Test, LinearAlgebra
 
 test_f!(du, u, p) = (@. du = u^2 - p)
 test_f(u, p) = (@. u^2 - p)
@@ -41,10 +42,12 @@ __compatible(::Number, ::AbstractArray) = false
 __compatible(u::AbstractArray, p::AbstractArray) = size(u) == size(p)
 
 __compatible(u::Number, ::SciMLBase.AbstractNonlinearAlgorithm) = true
-__compatible(u::Number, ::Union{CMINPACK, NLsolveJL}) = true
+__compatible(u::Number, ::Union{CMINPACK, NLsolveJL, KINSOL}) = true
 __compatible(u::AbstractArray, ::SciMLBase.AbstractNonlinearAlgorithm) = true
+__compatible(u::AbstractArray{T, N}, ::KINSOL) where {T, N} = N == 1  # Need to be fixed upstream
+__compatible(u::StaticArray{S, T, N}, ::KINSOL) where {S <: Tuple, T, N} = false
 __compatible(u::StaticArray, ::SciMLBase.AbstractNonlinearAlgorithm) = true
-__compatible(u::StaticArray, ::Union{CMINPACK, NLsolveJL}) = false
+__compatible(u::StaticArray, ::Union{CMINPACK, NLsolveJL, KINSOL}) = false
 __compatible(u, ::Nothing) = true
 
 __compatible(::Any, ::Any) = true
@@ -52,10 +55,12 @@ __compatible(::CMINPACK, ::Val{:iip_cache}) = false
 __compatible(::CMINPACK, ::Val{:oop_cache}) = false
 __compatible(::NLsolveJL, ::Val{:iip_cache}) = false
 __compatible(::NLsolveJL, ::Val{:oop_cache}) = false
+__compatible(::KINSOL, ::Val{:iip_cache}) = false
+__compatible(::KINSOL, ::Val{:oop_cache}) = false
 
 @testset "ForwardDiff.jl Integration: $(alg)" for alg in (NewtonRaphson(), TrustRegion(),
     LevenbergMarquardt(), PseudoTransient(; alpha_initial = 10.0), Broyden(), Klement(),
-    DFSane(), nothing, NLsolveJL(), CMINPACK())
+    DFSane(), nothing, NLsolveJL(), CMINPACK(), KINSOL())
     us = (2.0, @SVector[1.0, 1.0], [1.0, 1.0], ones(2, 2), @SArray ones(2, 2))
 
     @testset "Scalar AD" begin
diff --git a/test/core/nlls.jl b/test/core/nlls.jl
new file mode 100644
index 000000000..331f84faa
--- /dev/null
+++ b/test/core/nlls.jl
@@ -0,0 +1,87 @@
+using NonlinearSolve,
+    LinearSolve, LinearAlgebra, Test, StableRNGs, Random, ForwardDiff, Zygote
+
+true_function(x, θ) = @. θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4])
+true_function(y, x, θ) = (@. y = θ[1] * exp(θ[2] * x) * cos(θ[3] * x + θ[4]))
+
+θ_true = [1.0, 0.1, 2.0, 0.5]
+
+x = [-1.0, -0.5, 0.0, 0.5, 1.0]
+
+y_target = true_function(x, θ_true)
+
+function loss_function(θ, p)
+    ŷ = true_function(p, θ)
+    return ŷ .- y_target
+end
+
+function loss_function(resid, θ, p)
+    true_function(resid, p, θ)
+    resid .= resid .- y_target
+    return resid
+end
+
+θ_init = θ_true .+ randn!(StableRNG(0), similar(θ_true)) * 0.1
+prob_oop = NonlinearLeastSquaresProblem{false}(loss_function, θ_init, x)
+prob_iip = NonlinearLeastSquaresProblem(NonlinearFunction(loss_function;
+        resid_prototype = zero(y_target)), θ_init, x)
+
+nlls_problems = [prob_oop, prob_iip]
+
+solvers = []
+for linsolve in [nothing, LUFactorization(), KrylovJL_GMRES()]
+    vjp_autodiffs = linsolve isa KrylovJL ? [nothing, AutoZygote(), AutoFiniteDiff()] :
+                    [nothing]
+    for linesearch in [Static(), BackTracking(), HagerZhang(), StrongWolfe(), MoreThuente()],
+        vjp_autodiff in vjp_autodiffs
+
+        push!(solvers, GaussNewton(; linsolve, linesearch, vjp_autodiff))
+    end
+end
+append!(solvers,
+    [
+        LevenbergMarquardt(),
+        LevenbergMarquardt(; linsolve = LUFactorization()),
+        nothing,
+    ])
+for radius_update_scheme in [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.NocedalWright,
+    RadiusUpdateSchemes.NLsolve, RadiusUpdateSchemes.Hei, RadiusUpdateSchemes.Yuan,
+    RadiusUpdateSchemes.Fan, RadiusUpdateSchemes.Bastin]
+    push!(solvers, TrustRegion(; radius_update_scheme))
+end
+
+for prob in nlls_problems, solver in solvers
+    @time sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
+    @test SciMLBase.successful_retcode(sol)
+    @test maximum(abs, sol.resid) < 1e-6
+end
+
+# This is just for testing that we can use vjp provided by the user
+function vjp(v, θ, p)
+    resid = zeros(length(p))
+    J = ForwardDiff.jacobian((resid, θ) -> loss_function(resid, θ, p), resid, θ)
+    return vec(v' * J)
+end
+
+function vjp!(Jv, v, θ, p)
+    resid = zeros(length(p))
+    J = ForwardDiff.jacobian((resid, θ) -> loss_function(resid, θ, p), resid, θ)
+    mul!(vec(Jv), v', J)
+    return nothing
+end
+
+probs = [
+    NonlinearLeastSquaresProblem(NonlinearFunction{true}(loss_function;
+            resid_prototype = zero(y_target), vjp = vjp!), θ_init, x),
+    NonlinearLeastSquaresProblem(NonlinearFunction{false}(loss_function;
+            resid_prototype = zero(y_target), vjp = vjp), θ_init, x),
+]
+
+for prob in probs, solver in solvers
+    !(solver isa GaussNewton) && continue
+    !(solver.linsolve isa KrylovJL) && continue
+    @test_warn "Currently we don't make use of user provided `jvp`. This is planned to be \
+    fixed in the near future." sol=solve(prob, solver; maxiters = 10000, abstol = 1e-8)
+    sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
+    @test maximum(abs, sol.resid) < 1e-6
+end
diff --git a/test/basictests.jl b/test/core/rootfind.jl
similarity index 100%
rename from test/basictests.jl
rename to test/core/rootfind.jl
diff --git a/test/GPU/Project.toml b/test/gpu/Project.toml
similarity index 100%
rename from test/GPU/Project.toml
rename to test/gpu/Project.toml
diff --git a/test/gpu.jl b/test/gpu/core.jl
similarity index 100%
rename from test/gpu.jl
rename to test/gpu/core.jl
diff --git a/test/minpack.jl b/test/minpack.jl
deleted file mode 100644
index 46cc6be6e..000000000
--- a/test/minpack.jl
+++ /dev/null
@@ -1,35 +0,0 @@
-using NonlinearSolve, MINPACK, Test
-
-function f_iip(du, u, p)
-    du[1] = 2 - 2u[1]
-    du[2] = u[1] - 4u[2]
-end
-u0 = zeros(2)
-prob_iip = NonlinearProblem{true}(f_iip, u0)
-abstol = 1e-8
-
-for alg in [CMINPACK()]
-    local sol
-    sol = solve(prob_iip, alg)
-    @test sol.retcode == ReturnCode.Success
-    p = nothing
-
-    du = zeros(2)
-    f_iip(du, sol.u, nothing)
-    @test maximum(du) < 1e-6
-end
-
-# OOP Tests
-f_oop(u, p) = [2 - 2u[1], u[1] - 4u[2]]
-u0 = zeros(2)
-prob_oop = NonlinearProblem{false}(f_oop, u0)
-
-for alg in [CMINPACK()]
-    local sol
-    sol = solve(prob_oop, alg)
-    @test sol.retcode == ReturnCode.Success
-
-    du = zeros(2)
-    du = f_oop(sol.u, nothing)
-    @test maximum(du) < 1e-6
-end
diff --git a/test/misc.jl b/test/misc/banded_matrices.jl
similarity index 100%
rename from test/misc.jl
rename to test/misc/banded_matrices.jl
diff --git a/test/sparse.jl b/test/misc/bruss.jl
similarity index 100%
rename from test/sparse.jl
rename to test/misc/bruss.jl
diff --git a/test/infeasible.jl b/test/misc/infeasible.jl
similarity index 100%
rename from test/infeasible.jl
rename to test/misc/infeasible.jl
diff --git a/test/matrix_resizing.jl b/test/misc/matrix_resizing.jl
similarity index 100%
rename from test/matrix_resizing.jl
rename to test/misc/matrix_resizing.jl
diff --git a/test/polyalgs.jl b/test/misc/polyalgs.jl
similarity index 100%
rename from test/polyalgs.jl
rename to test/misc/polyalgs.jl
diff --git a/test/qa.jl b/test/misc/qa.jl
similarity index 87%
rename from test/qa.jl
rename to test/misc/qa.jl
index 960f7fc8b..9d123470d 100644
--- a/test/qa.jl
+++ b/test/misc/qa.jl
@@ -2,7 +2,7 @@ using NonlinearSolve, Aqua
 
 @testset "Aqua" begin
     Aqua.find_persistent_tasks_deps(NonlinearSolve)
-    Aqua.test_ambiguities(NonlinearSolve, recursive = false)
+    Aqua.test_ambiguities(NonlinearSolve; recursive = false)
     Aqua.test_deps_compat(NonlinearSolve)
     Aqua.test_piracies(NonlinearSolve,
         treat_as_own = [NonlinearProblem, NonlinearLeastSquaresProblem])
diff --git a/test/nlsolve.jl b/test/nlsolve.jl
deleted file mode 100644
index ebda5d272..000000000
--- a/test/nlsolve.jl
+++ /dev/null
@@ -1,138 +0,0 @@
-using NonlinearSolve, NLsolve, LinearAlgebra, Test
-
-# IIP Tests
-function f_iip(du, u, p, t)
-    du[1] = 2 - 2u[1]
-    du[2] = u[1] - 4u[2]
-end
-u0 = zeros(2)
-prob_iip = SteadyStateProblem(f_iip, u0)
-abstol = 1e-8
-
-for alg in [NLsolveJL()]
-    sol = solve(prob_iip, alg)
-    @test sol.retcode == ReturnCode.Success
-    p = nothing
-
-    du = zeros(2)
-    f_iip(du, sol.u, nothing, 0)
-    @test maximum(du) < 1e-6
-end
-
-# OOP Tests
-f_oop(u, p, t) = [2 - 2u[1], u[1] - 4u[2]]
-u0 = zeros(2)
-prob_oop = SteadyStateProblem(f_oop, u0)
-
-for alg in [NLsolveJL()]
-    sol = solve(prob_oop, alg)
-    @test sol.retcode == ReturnCode.Success
-    # test the solver is doing reasonable things for linear solve
-    # and that the stats are working properly
-    @test 1 <= sol.stats.nf < 10
-
-    du = zeros(2)
-    du = f_oop(sol.u, nothing, 0)
-    @test maximum(du) < 1e-6
-end
-
-# NonlinearProblem Tests
-
-function f_iip(du, u, p)
-    du[1] = 2 - 2u[1]
-    du[2] = u[1] - 4u[2]
-end
-u0 = zeros(2)
-prob_iip = NonlinearProblem{true}(f_iip, u0)
-abstol = 1e-8
-for alg in [NLsolveJL()]
-    local sol
-    sol = solve(prob_iip, alg)
-    @test sol.retcode == ReturnCode.Success
-    p = nothing
-
-    du = zeros(2)
-    f_iip(du, sol.u, nothing)
-    @test maximum(du) < 1e-6
-end
-
-# OOP Tests
-f_oop(u, p) = [2 - 2u[1], u[1] - 4u[2]]
-u0 = zeros(2)
-prob_oop = NonlinearProblem{false}(f_oop, u0)
-for alg in [NLsolveJL()]
-    local sol
-    sol = solve(prob_oop, alg)
-    @test sol.retcode == ReturnCode.Success
-
-    du = zeros(2)
-    du = f_oop(sol.u, nothing)
-    @test maximum(du) < 1e-6
-end
-
-# tolerance tests
-f_tol(u, p) = u^2 - 2
-prob_tol = NonlinearProblem(f_tol, 1.0)
-for tol in [1e-1, 1e-3, 1e-6, 1e-10, 1e-15]
-    sol = solve(prob_tol, NLsolveJL(), abstol = tol)
-    @test abs(sol.u[1] - sqrt(2)) < tol
-end
-
-# Test the finite differencing technique
-function f!(fvec, x, p)
-    fvec[1] = (x[1] + 3) * (x[2]^3 - 7) + 18
-    fvec[2] = sin(x[2] * exp(x[1]) - 1)
-end
-
-prob = NonlinearProblem{true}(f!, [0.1; 1.2])
-sol = solve(prob, NLsolveJL(autodiff = :central))
-
-du = zeros(2)
-f!(du, sol.u, nothing)
-@test maximum(du) < 1e-6
-
-# Test the autodiff technique
-function f!(fvec, x, p)
-    fvec[1] = (x[1] + 3) * (x[2]^3 - 7) + 18
-    fvec[2] = sin(x[2] * exp(x[1]) - 1)
-end
-
-prob = NonlinearProblem{true}(f!, [0.1; 1.2])
-sol = solve(prob, NLsolveJL(autodiff = :forward))
-
-du = zeros(2)
-f!(du, sol.u, nothing)
-@test maximum(du) < 1e-6
-
-function problem(x, A)
-    return x .^ 2 - A
-end
-
-function problemJacobian(x, A)
-    return diagm(2 .* x)
-end
-
-function f!(F, u, p)
-    F[1:152] = problem(u, p)
-end
-
-function j!(J, u, p)
-    J[1:152, 1:152] = problemJacobian(u, p)
-end
-
-f = NonlinearFunction(f!)
-
-init = ones(152);
-A = ones(152);
-A[6] = 0.8
-
-f = NonlinearFunction(f!, jac = j!)
-
-p = A
-
-ProbN = NonlinearProblem(f, init, p)
-sol = solve(ProbN, NLsolveJL(), reltol = 1e-8, abstol = 1e-8)
-
-init = ones(Complex{Float64}, 152);
-ProbN = NonlinearProblem(f, init, p)
-sol = solve(ProbN, NLsolveJL(), reltol = 1e-8, abstol = 1e-8)
diff --git a/test/runtests.jl b/test/runtests.jl
index c5932ac6d..2e1a241ed 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -1,48 +1,44 @@
-using Pkg
-using SafeTestsets
-const LONGER_TESTS = false
+using Pkg, SafeTestsets
 
 const GROUP = get(ENV, "GROUP", "All")
-const is_APPVEYOR = Sys.iswindows() && haskey(ENV, "APPVEYOR")
 
-function activate_downstream_env()
-    Pkg.activate("GPU")
+function activate_env(env)
+    Pkg.activate(env)
     Pkg.develop(PackageSpec(path = dirname(@__DIR__)))
     Pkg.instantiate()
 end
 
 @time begin
-    if GROUP == "All" || GROUP == "Core"
-        @time @safetestset "Basic Tests" include("basictests.jl")
-        @time @safetestset "Forward AD" include("forward_ad.jl")
+    if GROUP == "All" || GROUP == "RootFinding"
+        @time @safetestset "Basic Root Finding Tests" include("core/rootfind.jl")
+        @time @safetestset "Forward AD" include("core/forward_ad.jl")
     end
 
-    if GROUP == "All" || GROUP == "NLLS"
-        @time @safetestset "Nonlinear Least Squares" include("nonlinear_least_squares.jl")
+    if GROUP == "All" || GROUP == "NLLSSolvers"
+        @time @safetestset "Basic NLLS Solvers" include("core/nlls.jl")
     end
 
     if GROUP == "All" || GROUP == "Wrappers"
-        @time @safetestset "MINPACK" include("minpack.jl")
-        @time @safetestset "NLsolve" include("nlsolve.jl")
-        @time @safetestset "SIAMFANLEquations" include("siamfanlequations.jl")
-        @time @safetestset "SpeedMapping" include("speedmapping.jl")
-        @time @safetestset "FixedPointAcceleration" include("fixed_point_acceleration.jl")
+        @time @safetestset "Fixed Point Solvers" include("wrappers/fixedpoint.jl")
+        @time @safetestset "Root Finding Solvers" include("wrappers/rootfind.jl")
+        @time @safetestset "Nonlinear Least Squares Solvers" include("wrappers/nlls.jl")
     end
 
     if GROUP == "All" || GROUP == "23TestProblems"
-        @time @safetestset "23 Test Problems" include("23_test_problems.jl")
+        @time @safetestset "23 Test Problems" include("core/23_test_problems.jl")
     end
 
     if GROUP == "All" || GROUP == "Miscellaneous"
-        @time @safetestset "Quality Assurance" include("qa.jl")
-        @time @safetestset "Sparsity Tests" include("sparse.jl")
-        @time @safetestset "Polyalgs" include("polyalgs.jl")
-        @time @safetestset "Matrix Resizing" include("matrix_resizing.jl")
-        @time @safetestset "Infeasible Problems" include("infeasible.jl")
+        @time @safetestset "Quality Assurance" include("misc/qa.jl")
+        @time @safetestset "Sparsity Tests: Bruss Steady State" include("misc/bruss.jl")
+        @time @safetestset "Polyalgs" include("misc/polyalgs.jl")
+        @time @safetestset "Matrix Resizing" include("misc/matrix_resizing.jl")
+        @time @safetestset "Infeasible Problems" include("misc/infeasible.jl")
+        @time @safetestset "Banded Matrices" include("misc/banded_matrices.jl")
     end
 
     if GROUP == "GPU"
-        activate_downstream_env()
-        @time @safetestset "GPU Tests" include("gpu.jl")
+        activate_env("gpu")
+        @time @safetestset "GPU Tests" include("gpu/core.jl")
     end
 end
diff --git a/test/siamfanlequations.jl b/test/siamfanlequations.jl
deleted file mode 100644
index ed35485b2..000000000
--- a/test/siamfanlequations.jl
+++ /dev/null
@@ -1,152 +0,0 @@
-using NonlinearSolve, SIAMFANLEquations, LinearAlgebra, Test
-
-# IIP Tests
-function f_iip(du, u, p, t)
-    du[1] = 2 - 2u[1]
-    du[2] = u[1] - 4u[2]
-end
-u0 = zeros(2)
-prob_iip = SteadyStateProblem(f_iip, u0)
-abstol = 1e-8
-
-for alg in [SIAMFANLEquationsJL()]
-    sol = solve(prob_iip, alg)
-    @test sol.retcode == ReturnCode.Success
-    p = nothing
-
-    du = zeros(2)
-    f_iip(du, sol.u, nothing, 0)
-    @test maximum(du) < 1e-6
-end
-
-# OOP Tests
-f_oop(u, p, t) = [2 - 2u[1], u[1] - 4u[2]]
-u0 = zeros(2)
-prob_oop = SteadyStateProblem(f_oop, u0)
-
-for alg in [SIAMFANLEquationsJL()]
-    sol = solve(prob_oop, alg)
-    @test sol.retcode == ReturnCode.Success
-    # test the solver is doing reasonable things for linear solve
-    # and that the stats are working properly
-    @test 1 <= sol.stats.nf < 10
-
-    du = zeros(2)
-    du = f_oop(sol.u, nothing, 0)
-    @test maximum(du) < 1e-6
-end
-
-# NonlinearProblem Tests
-
-function f_iip(du, u, p)
-    du[1] = 2 - 2u[1]
-    du[2] = u[1] - 4u[2]
-end
-u0 = zeros(2)
-prob_iip = NonlinearProblem{true}(f_iip, u0)
-abstol = 1e-8
-for alg in [SIAMFANLEquationsJL()]
-    local sol
-    sol = solve(prob_iip, alg)
-    @test sol.retcode == ReturnCode.Success
-    p = nothing
-
-    du = zeros(2)
-    f_iip(du, sol.u, nothing)
-    @test maximum(du) < 1e-6
-end
-
-# OOP Tests
-f_oop(u, p) = [2 - 2u[1], u[1] - 4u[2]]
-u0 = zeros(2)
-prob_oop = NonlinearProblem{false}(f_oop, u0)
-for alg in [SIAMFANLEquationsJL()]
-    local sol
-    sol = solve(prob_oop, alg)
-    @test sol.retcode == ReturnCode.Success
-
-    du = zeros(2)
-    du = f_oop(sol.u, nothing)
-    @test maximum(du) < 1e-6
-end
-
-# tolerance tests for scalar equation solvers
-f_tol(u, p) = u^2 - 2
-prob_tol = NonlinearProblem(f_tol, 1.0)
-for tol in [1e-1, 1e-3, 1e-6, 1e-10, 1e-11]
-    for method = [:newton, :pseudotransient, :secant]
-        sol = solve(prob_tol, SIAMFANLEquationsJL(method = method), abstol = tol)
-        @test abs(sol.u[1] - sqrt(2)) < tol
-    end
-end
-
-# Test the JFNK technique
-f_jfnk(u, p) = u^2 - 2
-prob_jfnk = NonlinearProblem(f_jfnk, 1.0)
-for tol in [1e-1, 1e-3, 1e-6, 1e-10, 1e-11]
-    sol = solve(prob_jfnk, SIAMFANLEquationsJL(linsolve = :gmres), abstol = tol)
-    @test abs(sol.u[1] - sqrt(2)) < tol
-end
-
-# Test the finite differencing technique
-function f!(fvec, x, p)
-    fvec[1] = (x[1] + 3) * (x[2]^3 - 7) + 18
-    fvec[2] = sin(x[2] * exp(x[1]) - 1)
-end
-
-prob = NonlinearProblem{true}(f!, [0.1; 1.2])
-sol = solve(prob, SIAMFANLEquationsJL())
-
-du = zeros(2)
-f!(du, sol.u, nothing)
-@test maximum(du) < 1e-6
-
-# Test the autodiff technique
-function f!(fvec, x, p)
-    fvec[1] = (x[1] + 3) * (x[2]^3 - 7) + 18
-    fvec[2] = sin(x[2] * exp(x[1]) - 1)
-end
-
-prob = NonlinearProblem{true}(f!, [0.1; 1.2])
-sol = solve(prob, SIAMFANLEquationsJL())
-
-du = zeros(2)
-f!(du, sol.u, nothing)
-@test maximum(du) < 1e-6
-
-function problem(x, A)
-    return x .^ 2 - A
-end
-
-function problemJacobian(x, A)
-    return diagm(2 .* x)
-end
-
-function f!(F, u, p)
-    F[1:152] = problem(u, p)
-end
-
-function j!(J, u, p)
-    J[1:152, 1:152] = problemJacobian(u, p)
-end
-
-f = NonlinearFunction(f!)
-
-init = ones(152);
-A = ones(152);
-A[6] = 0.8
-
-f = NonlinearFunction(f!, jac = j!)
-
-p = A
-
-ProbN = NonlinearProblem(f, init, p)
-for method = [:newton, :pseudotransient]
-    sol = solve(ProbN, SIAMFANLEquationsJL(method = method), reltol = 1e-8, abstol = 1e-8)
-end
-
-#= doesn't support complex numbers handling
-init = ones(Complex{Float64}, 152);
-ProbN = NonlinearProblem(f, init, p)
-sol = solve(ProbN, SIAMFANLEquationsJL(), reltol = 1e-8, abstol = 1e-8)
-=#
\ No newline at end of file
diff --git a/test/speedmapping.jl b/test/speedmapping.jl
deleted file mode 100644
index 97bfed4ef..000000000
--- a/test/speedmapping.jl
+++ /dev/null
@@ -1,46 +0,0 @@
-using NonlinearSolve, SpeedMapping, LinearAlgebra, Test
-
-# Fixed Point for Power Method
-# Taken from https://github.com/NicolasL-S/SpeedMapping.jl/blob/95951db8f8a4457093090e18802ad382db1c76da/test/runtests.jl
-@testset "Power Method" begin
-    C = [1 2 3; 4 5 6; 7 8 9]
-    A = C + C'
-    B = Hermitian(ones(10) * ones(10)' .* im + Diagonal(1:10))
-
-    function power_method!(du, u, A)
-        mul!(du, A, u)
-        du ./= norm(du, Inf)
-        du .-= u  # Convert to a root finding problem
-        return nothing
-    end
-
-    prob = NonlinearProblem(power_method!, ones(3), A)
-
-    sol = solve(prob, SpeedMappingJL())
-    @test sol.u' * A[:, 3] ≈ 32.916472867168096
-
-    sol = solve(prob, SpeedMappingJL(; orders = [3, 2]))
-    @test sol.u' * A[:, 3] ≈ 32.916472867168096
-
-    sol = solve(prob, SpeedMappingJL(; stabilize = true))
-    @test sol.u' * A[:, 3] ≈ 32.91647286145264
-
-    # Non vector inputs
-    function power_method_nonvec!(du, u, A)
-        mul!(vec(du), A, vec(u))
-        du ./= norm(du, Inf)
-        du .-= u  # Convert to a root finding problem
-        return nothing
-    end
-
-    prob = NonlinearProblem(power_method_nonvec!, ones(1, 3, 1), A)
-
-    sol = solve(prob, SpeedMappingJL())
-    @test vec(sol.u)' * A[:, 3] ≈ 32.916472867168096
-
-    sol = solve(prob, SpeedMappingJL(; orders = [3, 2]))
-    @test vec(sol.u)' * A[:, 3] ≈ 32.916472867168096
-
-    sol = solve(prob, SpeedMappingJL(; stabilize = true))
-    @test vec(sol.u)' * A[:, 3] ≈ 32.91647286145264
-end
diff --git a/test/fixed_point_acceleration.jl b/test/wrappers/fixedpoint.jl
similarity index 56%
rename from test/fixed_point_acceleration.jl
rename to test/wrappers/fixedpoint.jl
index 099ccff9e..56b646378 100644
--- a/test/fixed_point_acceleration.jl
+++ b/test/wrappers/fixedpoint.jl
@@ -1,4 +1,4 @@
-using NonlinearSolve, FixedPointAcceleration, LinearAlgebra, Test
+using NonlinearSolve, FixedPointAcceleration, SpeedMapping, NLsolve, LinearAlgebra, Test
 
 # Simple Scalar Problem
 @testset "Simple Scalar Problem" begin
@@ -6,8 +6,14 @@ using NonlinearSolve, FixedPointAcceleration, LinearAlgebra, Test
     prob = NonlinearProblem(f1, 1.1)
 
     for alg in (:Anderson, :MPE, :RRE, :VEA, :SEA, :Simple, :Aitken, :Newton)
-        @test abs(solve(prob, FixedPointAccelerationJL()).resid) ≤ 1e-10
+        @test abs(solve(prob, FixedPointAccelerationJL(; algorithm = alg)).resid) ≤ 1e-10
     end
+
+    @test abs(solve(prob, SpeedMappingJL()).resid) ≤ 1e-10
+    @test abs(solve(prob, SpeedMappingJL(; orders = [3, 2])).resid) ≤ 1e-10
+    @test abs(solve(prob, SpeedMappingJL(; stabilize = true)).resid) ≤ 1e-10
+
+    @test abs(solve(prob, NLsolveJL(; method = :anderson)).resid) ≤ 1e-10
 end
 
 # Simple Vector Problem
@@ -18,6 +24,12 @@ end
     for alg in (:Anderson, :MPE, :RRE, :VEA, :SEA, :Simple, :Aitken, :Newton)
         @test maximum(abs.(solve(prob, FixedPointAccelerationJL()).resid)) ≤ 1e-10
     end
+
+    @test maximum(abs.(solve(prob, SpeedMappingJL()).resid)) ≤ 1e-10
+    @test maximum(abs.(solve(prob, SpeedMappingJL(; orders = [3, 2])).resid)) ≤ 1e-10
+    @test maximum(abs.(solve(prob, SpeedMappingJL(; stabilize = true)).resid)) ≤ 1e-10
+
+    @test_broken maximum(abs.(solve(prob, NLsolveJL(; method = :anderson)).resid)) ≤ 1e-10
 end
 
 # Fixed Point for Power Method
@@ -46,6 +58,15 @@ end
         end
     end
 
+    for kwargs in ((;), (; orders = [3, 2]), (; stabilize = true))
+        alg = SpeedMappingJL(; kwargs...)
+        sol = solve(prob, alg)
+        @test sol.u' * A[:, 3] ≈ 32.916472867168096
+    end
+
+    sol = solve(prob, NLsolveJL(; method = :anderson))
+    @test_broken sol.u' * A[:, 3] ≈ 32.916472867168096
+
     # Non vector inputs
     function power_method_nonvec!(du, u, A)
         mul!(vec(du), A, vec(u))
@@ -59,10 +80,19 @@ end
     for alg in (:Anderson, :MPE, :RRE, :VEA, :SEA, :Simple, :Aitken, :Newton)
         sol = solve(prob, FixedPointAccelerationJL(; algorithm = alg))
         if SciMLBase.successful_retcode(sol)
-            @test sol.u' * A[:, 3] ≈ 32.916472867168096
+            @test vec(sol.u)' * A[:, 3] ≈ 32.916472867168096
         else
             @warn "Power Method failed for FixedPointAccelerationJL(; algorithm = $alg)"
-            @test_broken sol.u' * A[:, 3] ≈ 32.916472867168096
+            @test_broken vec(sol.u)' * A[:, 3] ≈ 32.916472867168096
         end
     end
+
+    for kwargs in ((;), (; orders = [3, 2]), (; stabilize = true))
+        alg = SpeedMappingJL(; kwargs...)
+        sol = solve(prob, alg)
+        @test vec(sol.u)' * A[:, 3] ≈ 32.916472867168096
+    end
+
+    sol = solve(prob, NLsolveJL(; method = :anderson))
+    @test_broken vec(sol.u)' * A[:, 3] ≈ 32.916472867168096
 end
diff --git a/test/nonlinear_least_squares.jl b/test/wrappers/nlls.jl
similarity index 80%
rename from test/nonlinear_least_squares.jl
rename to test/wrappers/nlls.jl
index e97da43eb..3e31b47de 100644
--- a/test/nonlinear_least_squares.jl
+++ b/test/wrappers/nlls.jl
@@ -1,3 +1,4 @@
+
 using NonlinearSolve,
     LinearSolve, LinearAlgebra, Test, StableRNGs, Random, ForwardDiff, Zygote
 import FastLevenbergMarquardt, LeastSquaresOptim, MINPACK
@@ -29,29 +30,10 @@ prob_iip = NonlinearLeastSquaresProblem(NonlinearFunction(loss_function;
 
 nlls_problems = [prob_oop, prob_iip]
 
-solvers = []
-for linsolve in [nothing, LUFactorization(), KrylovJL_GMRES()]
-    vjp_autodiffs = linsolve isa KrylovJL ? [nothing, AutoZygote(), AutoFiniteDiff()] :
-                    [nothing]
-    for linesearch in [Static(), BackTracking(), HagerZhang(), StrongWolfe(), MoreThuente()],
-        vjp_autodiff in vjp_autodiffs
-
-        push!(solvers, GaussNewton(; linsolve, linesearch, vjp_autodiff))
-    end
-end
-append!(solvers,
-    [
-        LevenbergMarquardt(),
-        LevenbergMarquardt(; linsolve = LUFactorization()),
-        LeastSquaresOptimJL(:lm),
-        LeastSquaresOptimJL(:dogleg),
-        nothing,
-    ])
-for radius_update_scheme in [RadiusUpdateSchemes.Simple, RadiusUpdateSchemes.NocedalWright,
-    RadiusUpdateSchemes.NLsolve, RadiusUpdateSchemes.Hei, RadiusUpdateSchemes.Yuan,
-    RadiusUpdateSchemes.Fan, RadiusUpdateSchemes.Bastin]
-    push!(solvers, TrustRegion(; radius_update_scheme))
-end
+solvers = [
+    LeastSquaresOptimJL(:lm),
+    LeastSquaresOptimJL(:dogleg),
+]
 
 for prob in nlls_problems, solver in solvers
     @time sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
@@ -110,7 +92,7 @@ push!(solvers, CMINPACK())
 
 for solver in solvers, prob in probs
     @time sol = solve(prob, solver; maxiters = 10000, abstol = 1e-8)
-    @test norm(sol.resid) < 1e-6
+    @test maximum(abs, sol.resid) < 1e-6
 end
 
 probs = [
diff --git a/test/wrappers/rootfind.jl b/test/wrappers/rootfind.jl
new file mode 100644
index 000000000..47021fb77
--- /dev/null
+++ b/test/wrappers/rootfind.jl
@@ -0,0 +1,115 @@
+using NonlinearSolve, LinearAlgebra, SciMLBase, Test
+import NLsolve, SIAMFANLEquations, MINPACK
+
+@testset "Steady State Problems" begin
+    # IIP Tests
+    function f_iip(du, u, p, t)
+        du[1] = 2 - 2u[1]
+        du[2] = u[1] - 4u[2]
+    end
+    u0 = zeros(2)
+    prob_iip = SteadyStateProblem(f_iip, u0)
+
+    for alg in [NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
+        sol = solve(prob_iip, alg)
+        @test SciMLBase.successful_retcode(sol.retcode)
+        @test maximum(abs, sol.resid) < 1e-6
+    end
+
+    # OOP Tests
+    f_oop(u, p, t) = [2 - 2u[1], u[1] - 4u[2]]
+    u0 = zeros(2)
+    prob_oop = SteadyStateProblem(f_oop, u0)
+
+    for alg in [NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
+        sol = solve(prob_oop, alg)
+        @test SciMLBase.successful_retcode(sol.retcode)
+        @test maximum(abs, sol.resid) < 1e-6
+    end
+end
+
+@testset "Nonlinear Root Finding Problems" begin
+    # IIP Tests
+    function f_iip(du, u, p)
+        du[1] = 2 - 2u[1]
+        du[2] = u[1] - 4u[2]
+    end
+    u0 = zeros(2)
+    prob_iip = NonlinearProblem{true}(f_iip, u0)
+
+    for alg in [NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
+        local sol
+        sol = solve(prob_iip, alg)
+        @test SciMLBase.successful_retcode(sol.retcode)
+        @test maximum(abs, sol.resid) < 1e-6
+    end
+
+    # OOP Tests
+    f_oop(u, p) = [2 - 2u[1], u[1] - 4u[2]]
+    u0 = zeros(2)
+    prob_oop = NonlinearProblem{false}(f_oop, u0)
+    for alg in [NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL()]
+        local sol
+        sol = solve(prob_oop, alg)
+        @test SciMLBase.successful_retcode(sol.retcode)
+        @test maximum(abs, sol.resid) < 1e-6
+    end
+
+    # Tolerance Tests
+    f_tol(u, p) = u^2 - 2
+    prob_tol = NonlinearProblem(f_tol, 1.0)
+    for tol in [1e-1, 1e-3, 1e-6, 1e-10, 1e-15],
+        alg in [NLsolveJL(), CMINPACK(), SIAMFANLEquationsJL(; method = :newton),
+            SIAMFANLEquationsJL(; method = :pseudotransient),
+            SIAMFANLEquationsJL(; method = :secant)]
+
+        sol = solve(prob_tol, alg, abstol = tol)
+        @test abs(sol.u[1] - sqrt(2)) < tol
+    end
+
+    f_jfnk(u, p) = u^2 - 2
+    prob_jfnk = NonlinearProblem(f_jfnk, 1.0)
+    for tol in [1e-1, 1e-3, 1e-6, 1e-10, 1e-11]
+        sol = solve(prob_jfnk, SIAMFANLEquationsJL(linsolve = :gmres), abstol = tol)
+        @test abs(sol.u[1] - sqrt(2)) < tol
+    end
+
+    # Test the finite differencing technique
+    function f!(fvec, x, p)
+        fvec[1] = (x[1] + 3) * (x[2]^3 - 7) + 18
+        fvec[2] = sin(x[2] * exp(x[1]) - 1)
+    end
+
+    prob = NonlinearProblem{true}(f!, [0.1; 1.2])
+    sol = solve(prob, NLsolveJL(autodiff = :central))
+    @test maximum(abs, sol.resid) < 1e-6
+    sol = solve(prob, SIAMFANLEquationsJL())
+    @test maximum(abs, sol.resid) < 1e-6
+
+    # Test the autodiff technique
+    sol = solve(prob, NLsolveJL(autodiff = :forward))
+    @test maximum(abs, sol.resid) < 1e-6
+
+    # Custom Jacobian
+    problem(x, A) = x .^ 2 - A
+    problemJacobian(x, A) = diagm(2 .* x)
+
+    f!(F, u, p) = (F[1:152] = problem(u, p))
+    j!(J, u, p) = (J[1:152, 1:152] = problemJacobian(u, p))
+
+    init = ones(152)
+    A = ones(152)
+    A[6] = 0.8
+
+    f = NonlinearFunction(f!; jac = j!)
+    p = A
+
+    ProbN = NonlinearProblem(f, init, p)
+
+    sol = solve(ProbN, NLsolveJL(); abstol = 1e-8)
+    @test maximum(abs, sol.resid) < 1e-6
+    sol = solve(ProbN, SIAMFANLEquationsJL(; method = :newton); abstol = 1e-8)
+    @test maximum(abs, sol.resid) < 1e-6
+    sol = solve(ProbN, SIAMFANLEquationsJL(; method = :pseudotransient); abstol = 1e-8)
+    @test maximum(abs, sol.resid) < 1e-6
+end