diff --git a/Project.toml b/Project.toml
index 184aed113..572dbb6fe 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,7 +9,6 @@ ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
-FastBroadcast = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 FiniteDiff = "6a86dc24-6348-571c-b903-95158fe2bd41"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -18,6 +17,7 @@ LineSearch = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 MaybeInplace = "bb5d69b7-63fc-4a16-80bd-7e42200c7bdb"
+NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -28,7 +28,6 @@ SciMLJacobianOperators = "19f34311-ddf3-4b8b-af20-060888a46c0e"
 SciMLOperators = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
 SimpleNonlinearSolve = "727e6d20-b764-4bd8-a329-72de5adea6c7"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
-SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"
 StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
 SymbolicIndexingInterface = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
@@ -69,7 +68,6 @@ DiffEqBase = "6.155.3"
 DifferentiationInterface = "0.6.1"
 Enzyme = "0.13.2"
 ExplicitImports = "1.5"
-FastBroadcast = "0.3.5"
 FastClosures = "0.3.2"
 FastLevenbergMarquardt = "0.1"
 FiniteDiff = "2.24"
@@ -85,11 +83,11 @@ LinearAlgebra = "1.10"
 LinearSolve = "2.35"
 MINPACK = "1.2"
 MaybeInplace = "0.1.4"
-ModelingToolkit = "9.41.0"
 NLSolvers = "0.5"
 NLsolve = "4.5"
 NaNMath = "1"
 NonlinearProblemLibrary = "0.1.2"
+NonlinearSolveBase = "1"
 OrdinaryDiffEqTsit5 = "1.1.0"
 Pkg = "1.10"
 PrecompileTools = "1.2"
@@ -141,6 +139,7 @@ Pkg = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
 SIAMFANLEquations = "084e46ad-d928-497d-ad5e-07fa361a48c4"
+SparseConnectivityTracer = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
 SpeedMapping = "f1835b91-879b-4a3f-a438-e4baacf14412"
 StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
@@ -149,4 +148,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f"
 
 [targets]
-test = ["Aqua", "BandedMatrices", "BenchmarkTools", "CUDA", "Enzyme", "ExplicitImports", "FastLevenbergMarquardt", "FixedPointAcceleration", "Hwloc", "InteractiveUtils", "LeastSquaresOptim", "LineSearches", "MINPACK", "NLSolvers", "NLsolve", "NaNMath", "NonlinearProblemLibrary", "OrdinaryDiffEqTsit5", "Pkg", "Random", "ReTestItems", "SIAMFANLEquations", "SpeedMapping", "StableRNGs", "StaticArrays", "Sundials", "Test", "Zygote"]
+test = ["Aqua", "BandedMatrices", "BenchmarkTools", "CUDA", "Enzyme", "ExplicitImports", "FastLevenbergMarquardt", "FixedPointAcceleration", "Hwloc", "InteractiveUtils", "LeastSquaresOptim", "LineSearches", "MINPACK", "NLSolvers", "NLsolve", "NaNMath", "NonlinearProblemLibrary", "OrdinaryDiffEqTsit5", "Pkg", "Random", "ReTestItems", "SIAMFANLEquations", "SparseConnectivityTracer", "SpeedMapping", "StableRNGs", "StaticArrays", "Sundials", "Test", "Zygote"]
diff --git a/src/NonlinearSolve.jl b/src/NonlinearSolve.jl
index a41027531..15e1ee3a9 100644
--- a/src/NonlinearSolve.jl
+++ b/src/NonlinearSolve.jl
@@ -12,14 +12,14 @@ using DiffEqBase: DiffEqBase, AbstractNonlinearTerminationMode,
                   NormTerminationMode, RelNormTerminationMode, RelSafeBestTerminationMode,
                   RelSafeTerminationMode, RelTerminationMode,
                   SimpleNonlinearSolveTerminationMode, SteadyStateDiffEqTerminationMode
-using FastBroadcast: @..
 using FastClosures: @closure
 using LazyArrays: LazyArrays, ApplyArray, cache
 using LinearAlgebra: LinearAlgebra, ColumnNorm, Diagonal, I, LowerTriangular, Symmetric,
                      UpperTriangular, axpy!, cond, diag, diagind, dot, issuccess, istril,
                      istriu, lu, mul!, norm, pinv, tril!, triu!
 using LineSearch: LineSearch, AbstractLineSearchAlgorithm, AbstractLineSearchCache,
-                  NoLineSearch, RobustNonMonotoneLineSearch, BackTracking, LineSearchesJL
+                  NoLineSearch, RobustNonMonotoneLineSearch, BackTracking, LineSearchesJL,
+                  LiFukushimaLineSearch
 using LinearSolve: LinearSolve, LUFactorization, QRFactorization,
                    needs_concrete_A, AbstractFactorization,
                    DefaultAlgorithmChoice, DefaultLinearSolver
@@ -47,7 +47,6 @@ using SciMLJacobianOperators: AbstractJacobianOperator, JacobianOperator, VecJac
 
 ## Sparse AD Support
 using SparseArrays: AbstractSparseMatrix, SparseMatrixCSC
-using SparseConnectivityTracer: TracerSparsityDetector # This can be dropped in the next release
 using SparseMatrixColorings: ConstantColoringAlgorithm, GreedyColoringAlgorithm,
                              LargestFirst
 
@@ -182,7 +181,8 @@ export NewtonDescent, SteepestDescent, Dogleg, DampedNewtonDescent, GeodesicAcce
 
 # Globalization
 ## Line Search Algorithms
-export LineSearch, BackTracking, NoLineSearch, RobustNonMonotoneLineSearch, LineSearchesJL
+export LineSearch, BackTracking, NoLineSearch, RobustNonMonotoneLineSearch,
+       LiFukushimaLineSearch, LineSearchesJL
 ## Trust Region Algorithms
 export RadiusUpdateSchemes
 
diff --git a/src/algorithms/levenberg_marquardt.jl b/src/algorithms/levenberg_marquardt.jl
index e8a2fd0d7..501a5dd29 100644
--- a/src/algorithms/levenberg_marquardt.jl
+++ b/src/algorithms/levenberg_marquardt.jl
@@ -156,18 +156,16 @@ end
 @inline function __update_LM_diagonal!!(y::Diagonal, x::AbstractMatrix)
     if __can_setindex(y.diag)
         if fast_scalar_indexing(y.diag)
-            @inbounds for i in axes(x, 1)
-                y.diag[i] = max(y.diag[i], x[i, i])
+            @simd for i in axes(x, 1)
+                @inbounds y.diag[i] = max(y.diag[i], x[i, i])
             end
             return y
         else
-            idxs = diagind(x)
-            @.. broadcast=false y.diag=max(y.diag, @view(x[idxs]))
+            y .= max.(y.diag, @view(x[diagind(x)]))
             return y
         end
     else
-        idxs = diagind(x)
-        return Diagonal(@.. broadcast=false max(y.diag, @view(x[idxs])))
+        return Diagonal(max.(y.diag, @view(x[diagind(x)])))
     end
 end
 
diff --git a/src/descent/damped_newton.jl b/src/descent/damped_newton.jl
index 77b204b13..ba3e1d028 100644
--- a/src/descent/damped_newton.jl
+++ b/src/descent/damped_newton.jl
@@ -226,12 +226,12 @@ end
     if __can_setindex(J_cache)
         copyto!(J_cache, J)
         if fast_scalar_indexing(J_cache)
-            @inbounds for i in axes(J_cache, 1)
-                J_cache[i, i] += D[i, i]
+            @simd for i in axes(J_cache, 1)
+                @inbounds J_cache[i, i] += D[i, i]
             end
         else
             idxs = diagind(J_cache)
-            @.. broadcast=false @view(J_cache[idxs])=@view(J[idxs]) + @view(D[idxs])
+            J_cache[idxs] .= @view(J[idxs]) .+ @view(D[idxs])
         end
         return J_cache
     else
@@ -242,12 +242,12 @@ end
     if __can_setindex(J_cache)
         copyto!(J_cache, J)
         if fast_scalar_indexing(J_cache)
-            @inbounds for i in axes(J_cache, 1)
-                J_cache[i, i] += D
+            @simd for i in axes(J_cache, 1)
+                @inbounds J_cache[i, i] += D
             end
         else
             idxs = diagind(J_cache)
-            @.. broadcast=false @view(J_cache[idxs])=@view(J[idxs]) + D
+            J_cache[idxs] .= @view(J[idxs]) .+ D
         end
         return J_cache
     else
diff --git a/src/internal/approximate_initialization.jl b/src/internal/approximate_initialization.jl
index 5ce348ae3..985234ac4 100644
--- a/src/internal/approximate_initialization.jl
+++ b/src/internal/approximate_initialization.jl
@@ -20,11 +20,11 @@ end
 function (::DiagonalStructure)(J::AbstractVector, J_new::AbstractMatrix)
     if __can_setindex(J)
         if fast_scalar_indexing(J)
-            @inbounds for i in eachindex(J)
-                J[i] = J_new[i, i]
+            @simd for i in eachindex(J)
+                @inbounds J[i] = J_new[i, i]
             end
         else
-            @.. broadcast=false J=@view(J_new[diagind(J_new)])
+            J .= @view(J_new[diagind(J_new)])
         end
         return J
     end
diff --git a/test/runtests.jl b/test/runtests.jl
index bf2a8ecdb..33ca5da7a 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -7,7 +7,7 @@ const GROUP = lowercase(get(ENV, "GROUP", "All"))
 const EXTRA_PKGS = Pkg.PackageSpec[]
 (GROUP == "all" || GROUP == "downstream") &&
     push!(EXTRA_PKGS, Pkg.PackageSpec("ModelingToolkit"))
-Pkg.add(EXTRA_PKGS)
+length(EXTRA_PKGS) ≥ 1 && Pkg.add(EXTRA_PKGS)
 
 const RETESTITEMS_NWORKERS = parse(
     Int, get(ENV, "RETESTITEMS_NWORKERS", string(min(Hwloc.num_physical_cores(), 4))))
@@ -19,4 +19,4 @@ const RETESTITEMS_NWORKER_THREADS = parse(Int,
 
 ReTestItems.runtests(NonlinearSolve; tags = (GROUP == "all" ? nothing : [Symbol(GROUP)]),
     nworkers = RETESTITEMS_NWORKERS, nworker_threads = RETESTITEMS_NWORKER_THREADS,
-    testitem_timeout = 3600, retries = 4)
+    testitem_timeout = 3600)