diff --git a/src/algorithms/klement.jl b/src/algorithms/klement.jl
index 428dc382c..66daec1b7 100644
--- a/src/algorithms/klement.jl
+++ b/src/algorithms/klement.jl
@@ -25,8 +25,8 @@ over this.
         differentiable problems.
 """
 function Klement(; max_resets::Int = 100, linsolve = nothing, alpha = nothing,
-        linesearch = NoLineSearch(), precs = DEFAULT_PRECS, autodiff = nothing,
-        init_jacobian::Val{IJ} = Val(:identity)) where {IJ}
+        linesearch = NoLineSearch(), precs = DEFAULT_PRECS,
+        autodiff = nothing, init_jacobian::Val{IJ} = Val(:identity)) where {IJ}
     if !(linesearch isa AbstractNonlinearSolveLineSearchAlgorithm)
         Base.depwarn(
             "Passing in a `LineSearches.jl` algorithm directly is deprecated. \
diff --git a/src/core/approximate_jacobian.jl b/src/core/approximate_jacobian.jl
index d0471e99b..4204b2db7 100644
--- a/src/core/approximate_jacobian.jl
+++ b/src/core/approximate_jacobian.jl
@@ -66,8 +66,8 @@ function ApproximateJacobianSolveAlgorithm{concrete_jac, name}(;
         linesearch = LineSearchesJL(; method = linesearch)
     end
     return ApproximateJacobianSolveAlgorithm{concrete_jac, name}(
-        linesearch, trustregion, descent, update_rule, reinit_rule,
-        max_resets, max_shrink_times, initialization)
+        linesearch, trustregion, descent, update_rule,
+        reinit_rule, max_resets, max_shrink_times, initialization)
 end
 
 @inline concrete_jac(::ApproximateJacobianSolveAlgorithm{CJ}) where {CJ} = CJ
diff --git a/src/default.jl b/src/default.jl
index 4d40fb995..931b5bf1f 100644
--- a/src/default.jl
+++ b/src/default.jl
@@ -1,6 +1,7 @@
 # Poly Algorithms
 """
-    NonlinearSolvePolyAlgorithm(algs, ::Val{pType} = Val(:NLS)) where {pType}
+    NonlinearSolvePolyAlgorithm(algs, ::Val{pType} = Val(:NLS);
+        start_index = 1) where {pType}
 
 A general way to define PolyAlgorithms for `NonlinearProblem` and
 `NonlinearLeastSquaresProblem`. This is a container for a tuple of algorithms that will be
@@ -15,6 +16,10 @@ residual is returned.
     `NonlinearLeastSquaresProblem`. This is used to determine the correct problem type to
     dispatch on.
 
+### Keyword Arguments
+
+  - `start_index`: the index to start at. Defaults to `1`.
+
 ### Example
 
 ```julia
@@ -25,11 +30,14 @@ alg = NonlinearSolvePolyAlgorithm((NewtonRaphson(), Broyden()))
 """
 struct NonlinearSolvePolyAlgorithm{pType, N, A} <: AbstractNonlinearSolveAlgorithm{:PolyAlg}
     algs::A
+    start_index::Int
 
-    function NonlinearSolvePolyAlgorithm(algs, ::Val{pType} = Val(:NLS)) where {pType}
+    function NonlinearSolvePolyAlgorithm(
+            algs, ::Val{pType} = Val(:NLS); start_index::Int = 1) where {pType}
         @assert pType ∈ (:NLS, :NLLS)
+        @assert 0 < start_index ≤ length(algs)
         algs = Tuple(algs)
-        return new{pType, length(algs), typeof(algs)}(algs)
+        return new{pType, length(algs), typeof(algs)}(algs, start_index)
     end
 end
 
@@ -73,7 +81,7 @@ end
 
 function reinit_cache!(cache::NonlinearSolvePolyAlgorithmCache, args...; kwargs...)
     foreach(c -> reinit_cache!(c, args...; kwargs...), cache.caches)
-    cache.current = 1
+    cache.current = cache.alg.start_index
     cache.nsteps = 0
     cache.total_time = 0.0
 end
@@ -91,7 +99,7 @@ for (probType, pType) in ((:NonlinearProblem, :NLS), (:NonlinearLeastSquaresProb
                     alg.algs),
                 alg,
                 -1,
-                1,
+                alg.start_index,
                 0,
                 0.0,
                 maxtime,
@@ -263,6 +271,7 @@ function RobustMultiNewton(::Type{T} = Float64; concrete_jac = nothing, linsolve
         algs = (TrustRegion(; concrete_jac, linsolve, precs, autodiff),
             TrustRegion(; concrete_jac, linsolve, precs, autodiff,
                 radius_update_scheme = RadiusUpdateSchemes.Bastin),
+            NewtonRaphson(; concrete_jac, linsolve, precs, autodiff),
             NewtonRaphson(; concrete_jac, linsolve, precs,
                 linesearch = LineSearchesJL(; method = BackTracking()), autodiff),
             TrustRegion(; concrete_jac, linsolve, precs,
@@ -276,7 +285,8 @@ end
 """
     FastShortcutNonlinearPolyalg(::Type{T} = Float64; concrete_jac = nothing,
         linsolve = nothing, precs = DEFAULT_PRECS, must_use_jacobian::Val = Val(false),
-        prefer_simplenonlinearsolve::Val{SA} = Val(false), autodiff = nothing) where {T}
+        prefer_simplenonlinearsolve::Val{SA} = Val(false), autodiff = nothing,
+        u0_len::Union{Int, Nothing} = nothing) where {T}
 
 A polyalgorithm focused on balancing speed and robustness. It first tries less robust methods
 for more performance and then tries more robust techniques if the faster ones fail.
@@ -285,12 +295,19 @@ for more performance and then tries more robust techniques if the faster ones fa
 
   - `T`: The eltype of the initial guess. It is only used to check if some of the algorithms
     are compatible with the problem type. Defaults to `Float64`.
+
+### Keyword Arguments
+
+  - `u0_len`: The length of the initial guess. If this is `nothing`, then the length of the
+    initial guess is not checked. If this is an integer and it is less than `25`, we use
+    jacobian based methods.
 """
 function FastShortcutNonlinearPolyalg(
         ::Type{T} = Float64; concrete_jac = nothing, linsolve = nothing,
         precs = DEFAULT_PRECS, must_use_jacobian::Val{JAC} = Val(false),
         prefer_simplenonlinearsolve::Val{SA} = Val(false),
-        autodiff = nothing) where {T, JAC, SA}
+        u0_len::Union{Int, Nothing} = nothing, autodiff = nothing) where {T, JAC, SA}
+    start_index = 1
     if JAC
         if __is_complex(T)
             algs = (NewtonRaphson(; concrete_jac, linsolve, precs, autodiff),)
@@ -312,6 +329,7 @@ function FastShortcutNonlinearPolyalg(
                     SimpleKlement(),
                     NewtonRaphson(; concrete_jac, linsolve, precs, autodiff))
             else
+                start_index = u0_len !== nothing ? (u0_len ≤ 25 ? 4 : 1) : 1
                 algs = (SimpleBroyden(),
                     Broyden(; init_jacobian = Val(:true_jacobian), autodiff),
                     SimpleKlement(),
@@ -327,6 +345,8 @@ function FastShortcutNonlinearPolyalg(
                     Klement(; linsolve, precs, autodiff),
                     NewtonRaphson(; concrete_jac, linsolve, precs, autodiff))
             else
+                # TODO: This number requires a bit rigorous testing
+                start_index = u0_len !== nothing ? (u0_len ≤ 25 ? 4 : 1) : 1
                 algs = (Broyden(; autodiff),
                     Broyden(; init_jacobian = Val(:true_jacobian), autodiff),
                     Klement(; linsolve, precs, autodiff),
@@ -339,7 +359,7 @@ function FastShortcutNonlinearPolyalg(
             end
         end
     end
-    return NonlinearSolvePolyAlgorithm(algs, Val(:NLS))
+    return NonlinearSolvePolyAlgorithm(algs, Val(:NLS); start_index)
 end
 
 """
diff --git a/test/core/forward_ad_tests.jl b/test/core/forward_ad_tests.jl
index a13259939..24258e157 100644
--- a/test/core/forward_ad_tests.jl
+++ b/test/core/forward_ad_tests.jl
@@ -64,8 +64,8 @@ end
 
 @testitem "ForwardDiff.jl Integration" setup=[ForwardADTesting] begin
     for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(),
-        PseudoTransient(; alpha_initial = 10.0), Broyden(), Klement(), DFSane(), nothing,
-        NLsolveJL(), CMINPACK(), KINSOL(; globalization_strategy = :LineSearch))
+        PseudoTransient(; alpha_initial = 10.0), Broyden(), Klement(), DFSane(),
+        nothing, NLsolveJL(), CMINPACK(), KINSOL(; globalization_strategy = :LineSearch))
         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/rootfind_tests.jl b/test/core/rootfind_tests.jl
index 86dc5ff2c..87ba7cd35 100644
--- a/test/core/rootfind_tests.jl
+++ b/test/core/rootfind_tests.jl
@@ -476,8 +476,8 @@ end
 
 @testitem "Broyden" setup=[CoreRootfindTesting] begin
     @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad)) Init Jacobian: $(init_jacobian) Update Rule: $(update_rule)" for lsmethod in (
-            Static(), StrongWolfe(), BackTracking(), HagerZhang(),
-            MoreThuente(), LiFukushimaLineSearch()),
+            Static(), StrongWolfe(), BackTracking(),
+            HagerZhang(), MoreThuente(), LiFukushimaLineSearch()),
         ad in (AutoFiniteDiff(), AutoZygote()),
         init_jacobian in (Val(:identity), Val(:true_jacobian)),
         update_rule in (Val(:good_broyden), Val(:bad_broyden), Val(:diagonal))
@@ -575,8 +575,8 @@ end
 
 @testitem "LimitedMemoryBroyden" setup=[CoreRootfindTesting] begin
     @testset "LineSearch: $(_nameof(lsmethod)) LineSearch AD: $(_nameof(ad))" for lsmethod in (
-            Static(), StrongWolfe(), BackTracking(), HagerZhang(),
-            MoreThuente(), LiFukushimaLineSearch()),
+            Static(), StrongWolfe(), BackTracking(),
+            HagerZhang(), MoreThuente(), LiFukushimaLineSearch()),
         ad in (AutoFiniteDiff(), AutoZygote())
 
         linesearch = LineSearchesJL(; method = lsmethod, autodiff = ad)
diff --git a/test/misc/matrix_resizing_tests.jl b/test/misc/matrix_resizing_tests.jl
index 1d17a696a..9210e3e6b 100644
--- a/test/misc/matrix_resizing_tests.jl
+++ b/test/misc/matrix_resizing_tests.jl
@@ -7,8 +7,8 @@
     vecprob = NonlinearProblem(ff, vec(u0), p)
     prob = NonlinearProblem(ff, u0, p)
 
-    for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(),
-        RobustMultiNewton(), FastShortcutNonlinearPolyalg(),
+    for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(),
+        PseudoTransient(), RobustMultiNewton(), FastShortcutNonlinearPolyalg(),
         Broyden(), Klement(), LimitedMemoryBroyden(; threshold = 2))
         @test vec(solve(prob, alg).u) == solve(vecprob, alg).u
     end
@@ -23,8 +23,8 @@ end
     vecprob = NonlinearProblem(fiip, vec(u0), p)
     prob = NonlinearProblem(fiip, u0, p)
 
-    for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(), PseudoTransient(),
-        RobustMultiNewton(), FastShortcutNonlinearPolyalg(),
+    for alg in (NewtonRaphson(), TrustRegion(), LevenbergMarquardt(),
+        PseudoTransient(), RobustMultiNewton(), FastShortcutNonlinearPolyalg(),
         Broyden(), Klement(), LimitedMemoryBroyden(; threshold = 2))
         @test vec(solve(prob, alg).u) == solve(vecprob, alg).u
     end