diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
index 9dc416799..af44cd769 100644
--- a/.github/workflows/Documentation.yml
+++ b/.github/workflows/Documentation.yml
@@ -22,7 +22,7 @@ jobs:
           Pkg.Registry.update()
           # Install packages present in subdirectories
           dev_pks = Pkg.PackageSpec[]
-          for path in ("lib/SciMLJacobianOperators", ".", "lib/SimpleNonlinearSolve", "lib/NonlinearSolveBase", "lib/BracketingNonlinearSolve")
+          for path in ("lib/SciMLJacobianOperators", ".", "lib/SimpleNonlinearSolve", "lib/NonlinearSolveBase", "lib/BracketingNonlinearSolve", "lib/NonlinearSolveFirstOrder", "lib/NonlinearSolveQuasiNewton", "lib/NonlinearSolveSpectralMethods")
               push!(dev_pks, Pkg.PackageSpec(; path))
           end
           Pkg.develop(dev_pks)
diff --git a/docs/Project.toml b/docs/Project.toml
index 9ed14da4e..09392f6f4 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -14,6 +14,9 @@ InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
+NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
+NonlinearSolveQuasiNewton = "9a2c21bd-3a47-402d-9113-8faf9a0ee114"
+NonlinearSolveSpectralMethods = "26075421-4e9a-44e1-8bd1-420ed7ad02b2"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 PETSc = "ace2c81b-2b5f-4b1e-a30d-d662738edfe0"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
@@ -42,8 +45,11 @@ InteractiveUtils = "<0.0.1, 1"
 LinearSolve = "2"
 NonlinearSolve = "4"
 NonlinearSolveBase = "1"
+NonlinearSolveFirstOrder = "1"
+NonlinearSolveQuasiNewton = "1"
+NonlinearSolveSpectralMethods = "1"
 OrdinaryDiffEqTsit5 = "1.1.0"
-PETSc = "0.2"
+PETSc = "0.3"
 Plots = "1"
 Random = "1.10"
 SciMLBase = "2.4"
diff --git a/docs/make.jl b/docs/make.jl
index bdc2b7fe9..29b1535dd 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -1,13 +1,23 @@
 using Documenter, DocumenterCitations, DocumenterInterLinks
-using NonlinearSolve, SimpleNonlinearSolve, Sundials, SteadyStateDiffEq, SciMLBase,
-      BracketingNonlinearSolve, NonlinearSolveBase
-using SciMLJacobianOperators
 import DiffEqBase
 
-cp(joinpath(@__DIR__, "Manifest.toml"),
-    joinpath(@__DIR__, "src/assets/Manifest.toml"), force = true)
-cp(joinpath(@__DIR__, "Project.toml"),
-    joinpath(@__DIR__, "src/assets/Project.toml"), force = true)
+using Sundials
+using NonlinearSolveBase, SciMLBase, DiffEqBase
+using SimpleNonlinearSolve, BracketingNonlinearSolve
+using NonlinearSolveFirstOrder, NonlinearSolveQuasiNewton, NonlinearSolveSpectralMethods
+using SciMLJacobianOperators
+using NonlinearSolve, SteadyStateDiffEq
+
+cp(
+    joinpath(@__DIR__, "Manifest.toml"),
+    joinpath(@__DIR__, "src/assets/Manifest.toml");
+    force = true
+)
+cp(
+    joinpath(@__DIR__, "Project.toml"),
+    joinpath(@__DIR__, "src/assets/Project.toml");
+    force = true
+)
 
 include("pages.jl")
 
@@ -20,20 +30,29 @@ interlinks = InterLinks(
 
 makedocs(;
     sitename = "NonlinearSolve.jl",
-    authors = "Chris Rackauckas",
-    modules = [NonlinearSolve, SimpleNonlinearSolve, SteadyStateDiffEq, DiffEqBase,
-        Sundials, NonlinearSolveBase, SciMLBase, SciMLJacobianOperators,
-        BracketingNonlinearSolve],
+    authors = "SciML",
+    modules = [
+        NonlinearSolveBase, SciMLBase, DiffEqBase,
+        SimpleNonlinearSolve, BracketingNonlinearSolve,
+        NonlinearSolveFirstOrder, NonlinearSolveQuasiNewton, NonlinearSolveSpectralMethods,
+        Sundials,
+        SciMLJacobianOperators,
+        NonlinearSolve, SteadyStateDiffEq
+    ],
     clean = true,
     doctest = false,
     linkcheck = true,
-    linkcheck_ignore = ["https://twitter.com/ChrisRackauckas/status/1544743542094020615",
-        "https://link.springer.com/article/10.1007/s40096-020-00339-4"],
+    linkcheck_ignore = [
+        "https://twitter.com/ChrisRackauckas/status/1544743542094020615",
+        "https://link.springer.com/article/10.1007/s40096-020-00339-4"
+    ],
     checkdocs = :exports,
     warnonly = [:missing_docs],
     plugins = [bib, interlinks],
-    format = Documenter.HTML(assets = ["assets/favicon.ico", "assets/citations.css"],
-        canonical = "https://docs.sciml.ai/NonlinearSolve/stable/"),
+    format = Documenter.HTML(
+        assets = ["assets/favicon.ico", "assets/citations.css"],
+        canonical = "https://docs.sciml.ai/NonlinearSolve/stable/"
+    ),
     pages
 )
 
diff --git a/docs/src/basics/autodiff.md b/docs/src/basics/autodiff.md
index d2d01e00b..5cdffa75b 100644
--- a/docs/src/basics/autodiff.md
+++ b/docs/src/basics/autodiff.md
@@ -3,7 +3,7 @@
 !!! note
     
     We support all backends supported by DifferentiationInterface.jl. Please refer to
-    the [backends page](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/stable/explanation/backends/)
+    the [backends page](https://juliadiff.org/DifferentiationInterface.jl/DifferentiationInterface/stable/explanation/backends/)
     for more information.
 
 ## Summary of Finite Differencing Backends
diff --git a/docs/src/basics/faq.md b/docs/src/basics/faq.md
index 265cee264..4d428250a 100644
--- a/docs/src/basics/faq.md
+++ b/docs/src/basics/faq.md
@@ -102,7 +102,8 @@ It is hard to say why your code is not fast. Take a look at the
 there is type instability.
 
 If you are using the defaults for the autodiff and your problem is not a scalar or using
-static arrays, ForwardDiff will create type unstable code. See this simple example:
+static arrays, ForwardDiff will create type unstable code and lead to dynamic dispatch
+internally. See this simple example:
 
 ```@example type_unstable
 using NonlinearSolve, InteractiveUtils
@@ -136,14 +137,17 @@ prob = NonlinearProblem(f, [1.0, 2.0], 2.0)
 nothing # hide
 ```
 
-Oh no! This is type unstable. This is because ForwardDiff.jl will chunk the jacobian
-computation and the type of this chunksize can't be statically inferred. To fix this, we
-directly specify the chunksize:
+Ah it is still type stable. But internally since the chunksize is not statically inferred,
+it will be dynamic and lead to dynamic dispatch. To fix this, we directly specify the
+chunksize:
 
 ```@example type_unstable
-@code_warntype solve(prob,
+@code_warntype solve(
+    prob,
     NewtonRaphson(;
-        autodiff = AutoForwardDiff(; chunksize = NonlinearSolve.pickchunksize(prob.u0))))
+        autodiff = AutoForwardDiff(; chunksize = NonlinearSolve.pickchunksize(prob.u0))
+    )
+)
 nothing # hide
 ```
 
diff --git a/docs/src/devdocs/algorithm_helpers.md b/docs/src/devdocs/algorithm_helpers.md
index c945b6003..2bd038792 100644
--- a/docs/src/devdocs/algorithm_helpers.md
+++ b/docs/src/devdocs/algorithm_helpers.md
@@ -3,8 +3,8 @@
 ## Pseudo Transient Method
 
 ```@docs
-NonlinearSolve.SwitchedEvolutionRelaxation
-NonlinearSolve.SwitchedEvolutionRelaxationCache
+NonlinearSolveFirstOrder.SwitchedEvolutionRelaxation
+NonlinearSolveFirstOrder.SwitchedEvolutionRelaxationCache
 ```
 
 ## Approximate Jacobian Methods
@@ -12,54 +12,54 @@ NonlinearSolve.SwitchedEvolutionRelaxationCache
 ### Initialization
 
 ```@docs
-NonlinearSolve.IdentityInitialization
-NonlinearSolve.TrueJacobianInitialization
-NonlinearSolve.BroydenLowRankInitialization
+NonlinearSolveQuasiNewton.IdentityInitialization
+NonlinearSolveQuasiNewton.TrueJacobianInitialization
+NonlinearSolveQuasiNewton.BroydenLowRankInitialization
 ```
 
 ### Jacobian Structure
 
 ```@docs
-NonlinearSolve.FullStructure
-NonlinearSolve.DiagonalStructure
+NonlinearSolveQuasiNewton.FullStructure
+NonlinearSolveQuasiNewton.DiagonalStructure
 ```
 
 ### Jacobian Caches
 
 ```@docs
-NonlinearSolve.InitializedApproximateJacobianCache
+NonlinearSolveQuasiNewton.InitializedApproximateJacobianCache
 ```
 
 ### Reset Methods
 
 ```@docs
-NonlinearSolve.NoChangeInStateReset
-NonlinearSolve.IllConditionedJacobianReset
+NonlinearSolveQuasiNewton.NoChangeInStateReset
+NonlinearSolveQuasiNewton.IllConditionedJacobianReset
 ```
 
 ### Update Rules
 
 ```@docs
-NonlinearSolve.GoodBroydenUpdateRule
-NonlinearSolve.BadBroydenUpdateRule
-NonlinearSolve.KlementUpdateRule
+NonlinearSolveQuasiNewton.GoodBroydenUpdateRule
+NonlinearSolveQuasiNewton.BadBroydenUpdateRule
+NonlinearSolveQuasiNewton.KlementUpdateRule
 ```
 
 ## Levenberg Marquardt Method
 
 ```@docs
-NonlinearSolve.LevenbergMarquardtTrustRegion
+NonlinearSolveFirstOrder.LevenbergMarquardtTrustRegion
 ```
 
 ## Trust Region Method
 
 ```@docs
-NonlinearSolve.GenericTrustRegionScheme
+NonlinearSolveFirstOrder.GenericTrustRegionScheme
 ```
 
 ## Miscellaneous
 
 ```@docs
-NonlinearSolve.callback_into_cache!
-NonlinearSolve.concrete_jac
+NonlinearSolveBase.callback_into_cache!
+NonlinearSolveBase.concrete_jac
 ```
diff --git a/docs/src/devdocs/internal_interfaces.md b/docs/src/devdocs/internal_interfaces.md
index 53e20d754..f4474e1aa 100644
--- a/docs/src/devdocs/internal_interfaces.md
+++ b/docs/src/devdocs/internal_interfaces.md
@@ -3,50 +3,43 @@
 ## Solvers
 
 ```@docs
-NonlinearSolve.AbstractNonlinearSolveAlgorithm
-NonlinearSolve.AbstractNonlinearSolveExtensionAlgorithm
-NonlinearSolve.AbstractNonlinearSolveCache
+NonlinearSolveBase.AbstractNonlinearSolveAlgorithm
+NonlinearSolveBase.AbstractNonlinearSolveCache
 ```
 
-## Descent Algorithms
+## Descent Directions
 
 ```@docs
-NonlinearSolve.AbstractDescentAlgorithm
-NonlinearSolve.AbstractDescentCache
+NonlinearSolveBase.AbstractDescentDirection
+NonlinearSolveBase.AbstractDescentCache
 ```
 
-## Descent Results
+### Descent Results
 
 ```@docs
-NonlinearSolve.DescentResult
+NonlinearSolveBase.DescentResult
 ```
 
 ## Approximate Jacobian
 
 ```@docs
-NonlinearSolve.AbstractApproximateJacobianStructure
-NonlinearSolve.AbstractJacobianInitialization
-NonlinearSolve.AbstractApproximateJacobianUpdateRule
-NonlinearSolve.AbstractApproximateJacobianUpdateRuleCache
-NonlinearSolve.AbstractResetCondition
+NonlinearSolveBase.AbstractApproximateJacobianStructure
+NonlinearSolveBase.AbstractJacobianInitialization
+NonlinearSolveBase.AbstractApproximateJacobianUpdateRule
+NonlinearSolveBase.AbstractApproximateJacobianUpdateRuleCache
+NonlinearSolveBase.AbstractResetCondition
 ```
 
 ## Damping Algorithms
 
 ```@docs
-NonlinearSolve.AbstractDampingFunction
-NonlinearSolve.AbstractDampingFunctionCache
+NonlinearSolveBase.AbstractDampingFunction
+NonlinearSolveBase.AbstractDampingFunctionCache
 ```
 
 ## Trust Region
 
 ```@docs
-NonlinearSolve.AbstractTrustRegionMethod
-NonlinearSolve.AbstractTrustRegionMethodCache
-```
-
-## Tracing
-
-```@docs
-NonlinearSolve.AbstractNonlinearSolveTraceLevel
+NonlinearSolveBase.AbstractTrustRegionMethod
+NonlinearSolveBase.AbstractTrustRegionMethodCache
 ```
diff --git a/docs/src/devdocs/jacobian.md b/docs/src/devdocs/jacobian.md
index 082478f7c..3e54966f4 100644
--- a/docs/src/devdocs/jacobian.md
+++ b/docs/src/devdocs/jacobian.md
@@ -1,6 +1,5 @@
 # Jacobian Wrappers
 
 ```@docs
-NonlinearSolve.AbstractNonlinearSolveJacobianCache
-NonlinearSolve.JacobianCache
+NonlinearSolveBase.construct_jacobian_cache
 ```
diff --git a/docs/src/devdocs/linear_solve.md b/docs/src/devdocs/linear_solve.md
index 88fa87440..dadd7dea1 100644
--- a/docs/src/devdocs/linear_solve.md
+++ b/docs/src/devdocs/linear_solve.md
@@ -1,6 +1,6 @@
 # Linear Solve
 
 ```@docs
-NonlinearSolve.AbstractLinearSolverCache
-NonlinearSolve.LinearSolverCache
+NonlinearSolveBase.AbstractLinearSolverCache
+NonlinearSolveBase.construct_linear_solver
 ```
diff --git a/docs/src/devdocs/operators.md b/docs/src/devdocs/operators.md
index fcbadc778..249eda29d 100644
--- a/docs/src/devdocs/operators.md
+++ b/docs/src/devdocs/operators.md
@@ -3,5 +3,5 @@
 ## Low-Rank Jacobian Operators
 
 ```@docs
-NonlinearSolve.BroydenLowRankJacobian
+NonlinearSolveQuasiNewton.BroydenLowRankJacobian
 ```
diff --git a/docs/src/native/diagnostics.md b/docs/src/native/diagnostics.md
index 35f11552f..1c0571c87 100644
--- a/docs/src/native/diagnostics.md
+++ b/docs/src/native/diagnostics.md
@@ -5,9 +5,9 @@
 These functions are not exported since the names have a potential for conflict.
 
 ```@docs
-NonlinearSolve.enable_timer_outputs
-NonlinearSolve.disable_timer_outputs
-NonlinearSolve.@static_timeit
+NonlinearSolveBase.enable_timer_outputs
+NonlinearSolveBase.disable_timer_outputs
+NonlinearSolveBase.@static_timeit
 ```
 
 ## Tracing API
@@ -17,6 +17,3 @@ TraceAll
 TraceWithJacobianConditionNumber
 TraceMinimal
 ```
-
-For details about the arguments refer to the documentation of
-[`NonlinearSolve.AbstractNonlinearSolveTraceLevel`](@ref).
diff --git a/docs/src/native/solvers.md b/docs/src/native/solvers.md
index 5653f1fab..a3aa900e0 100644
--- a/docs/src/native/solvers.md
+++ b/docs/src/native/solvers.md
@@ -81,7 +81,7 @@ All of the previously mentioned solvers are wrappers around the following solver
 are meant for advanced users and allow building custom solvers.
 
 ```@docs
-ApproximateJacobianSolveAlgorithm
+QuasiNewtonAlgorithm
 GeneralizedFirstOrderAlgorithm
 GeneralizedDFSane
 ```
diff --git a/lib/NonlinearSolveBase/src/abstract_types.jl b/lib/NonlinearSolveBase/src/abstract_types.jl
index 91afc2901..6105f5125 100644
--- a/lib/NonlinearSolveBase/src/abstract_types.jl
+++ b/lib/NonlinearSolveBase/src/abstract_types.jl
@@ -542,7 +542,7 @@ Define custom operations on `internalcache` tightly coupled with the calling `ca
 `args...` contain the sequence of caches calling into `internalcache`.
 
 This unfortunately makes code very tightly coupled and not modular. It is recommended to not
-use this functionality unless it can't be avoided (like in [`LevenbergMarquardt`](@ref)).
+use this functionality unless it can't be avoided (like in `LevenbergMarquardt`).
 """
 callback_into_cache!(cache, internalcache, args...) = nothing  # By default do nothing
 
diff --git a/lib/NonlinearSolveBase/src/descent/geodesic_acceleration.jl b/lib/NonlinearSolveBase/src/descent/geodesic_acceleration.jl
index 409ed9ce9..f7758b1e7 100644
--- a/lib/NonlinearSolveBase/src/descent/geodesic_acceleration.jl
+++ b/lib/NonlinearSolveBase/src/descent/geodesic_acceleration.jl
@@ -5,7 +5,7 @@ Uses the `descent` algorithm to compute the velocity and acceleration terms for
 geodesic acceleration method. The velocity and acceleration terms are then combined to
 compute the descent direction.
 
-This method in its current form was developed for [`LevenbergMarquardt`](@ref). Performance
+This method in its current form was developed for `LevenbergMarquardt`. Performance
 for other methods are not theorectically or experimentally verified.
 
 ### Keyword Arguments
diff --git a/lib/NonlinearSolveBase/src/timer_outputs.jl b/lib/NonlinearSolveBase/src/timer_outputs.jl
index 0b6d23496..623e7ea15 100644
--- a/lib/NonlinearSolveBase/src/timer_outputs.jl
+++ b/lib/NonlinearSolveBase/src/timer_outputs.jl
@@ -18,7 +18,7 @@ end
     @static_timeit to name expr
 
 Like `TimerOutputs.@timeit_debug` but has zero overhead if `TimerOutputs` is disabled via
-[`NonlinearSolve.disable_timer_outputs()`](@ref).
+[`NonlinearSolveBase.disable_timer_outputs()`](@ref).
 """
 macro static_timeit(to, name, expr)
     @static if TIMER_OUTPUTS_ENABLED
diff --git a/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl b/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl
index a06378dcb..f9c11ae28 100644
--- a/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl
+++ b/lib/NonlinearSolveFirstOrder/src/levenberg_marquardt.jl
@@ -30,7 +30,7 @@ nonlinear systems.
     should not be disabled.
 
 For the remaining arguments, see [`GeodesicAcceleration`](@ref) and
-[`NonlinearSolve.LevenbergMarquardtTrustRegion`](@ref) documentations.
+[`NonlinearSolveFirstOrder.LevenbergMarquardtTrustRegion`](@ref) documentations.
 """
 function LevenbergMarquardt(;
         linsolve = nothing, precs = nothing,
diff --git a/lib/NonlinearSolveFirstOrder/src/solve.jl b/lib/NonlinearSolveFirstOrder/src/solve.jl
index b2711ca76..c9c8c77a8 100644
--- a/lib/NonlinearSolveFirstOrder/src/solve.jl
+++ b/lib/NonlinearSolveFirstOrder/src/solve.jl
@@ -15,9 +15,9 @@ order of convergence.
 ### Keyword Arguments
 
   - `trustregion`: Globalization using a Trust Region Method. This needs to follow the
-    [`NonlinearSolve.AbstractTrustRegionMethod`](@ref) interface.
+    [`NonlinearSolveBase.AbstractTrustRegionMethod`](@ref) interface.
   - `descent`: The descent method to use to compute the step. This needs to follow the
-    [`NonlinearSolve.AbstractDescentAlgorithm`](@ref) interface.
+    [`NonlinearSolveBase.AbstractDescentDirection`](@ref) interface.
   - `max_shrink_times`: The maximum number of times the trust region radius can be shrunk
     before the algorithm terminates.
 """
diff --git a/lib/NonlinearSolveFirstOrder/src/trust_region.jl b/lib/NonlinearSolveFirstOrder/src/trust_region.jl
index f76b77601..92d2655ac 100644
--- a/lib/NonlinearSolveFirstOrder/src/trust_region.jl
+++ b/lib/NonlinearSolveFirstOrder/src/trust_region.jl
@@ -19,7 +19,7 @@ for large-scale and numerically-difficult nonlinear systems.
     `RadiusUpdateSchemes.Simple`. See [`RadiusUpdateSchemes`](@ref) for more details. For a
     review on trust region radius update schemes, see [yuan2015recent](@citet).
 
-For the remaining arguments, see [`NonlinearSolve.GenericTrustRegionScheme`](@ref)
+For the remaining arguments, see [`NonlinearSolveFirstOrder.GenericTrustRegionScheme`](@ref)
 documentation.
 """
 function TrustRegion(;
diff --git a/lib/NonlinearSolveQuasiNewton/src/initialization.jl b/lib/NonlinearSolveQuasiNewton/src/initialization.jl
index 4af1fb921..f730b541b 100644
--- a/lib/NonlinearSolveQuasiNewton/src/initialization.jl
+++ b/lib/NonlinearSolveQuasiNewton/src/initialization.jl
@@ -10,7 +10,8 @@ A cache for Approximate Jacobian.
   - `J`: The current Jacobian.
   - `structure`: The structure of the Jacobian.
   - `alg`: The initialization algorithm.
-  - `cache`: The Jacobian cache [`NonlinearSolve.JacobianCache`](@ref) (if needed).
+  - `cache`: The Jacobian cache [`NonlinearSolveBase.construct_jacobian_cache`](@ref)
+    (if needed).
   - `initialized`: A boolean indicating whether the Jacobian has been initialized.
   - `internalnorm`: The norm to be used.
 
diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl
index 2ee0ac575..7d6ac7b09 100644
--- a/src/NonlinearSolve.jl
+++ b/src/NonlinearSolve.jl
@@ -46,8 +46,6 @@ using SimpleNonlinearSolve: SimpleNonlinearSolve
 
 const SII = SymbolicIndexingInterface
 
-include("helpers.jl")
-
 include("polyalg.jl")
 include("extension_algs.jl")
 
diff --git a/src/forward_diff.jl b/src/forward_diff.jl
index 28534c59e..410e818c3 100644
--- a/src/forward_diff.jl
+++ b/src/forward_diff.jl
@@ -87,3 +87,12 @@ end
 nodual_value(x) = x
 nodual_value(x::Dual) = ForwardDiff.value(x)
 nodual_value(x::AbstractArray{<:Dual}) = map(ForwardDiff.value, x)
+
+"""
+    pickchunksize(x) = pickchunksize(length(x))
+    pickchunksize(x::Int)
+
+Determine the chunk size for ForwardDiff and PolyesterForwardDiff based on the input length.
+"""
+@inline pickchunksize(x) = pickchunksize(length(x))
+@inline pickchunksize(x::Int) = ForwardDiff.pickchunksize(x)
diff --git a/src/helpers.jl b/src/helpers.jl
deleted file mode 100644
index 8b1378917..000000000
--- a/src/helpers.jl
+++ /dev/null
@@ -1 +0,0 @@
-