From 030f7754306615bb6041361e18a58d503fa69476 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sat, 19 Jul 2025 15:59:08 -0400
Subject: [PATCH 01/15] Explicit imports

---
 docs/src/api/SciMLJacobianOperators.md        |  2 +-
 docs/src/api/fastlevenbergmarquardt.md        |  5 +-
 docs/src/api/fixedpointacceleration.md        |  5 +-
 docs/src/api/homotopycontinuation.md          |  5 +-
 docs/src/api/leastsquaresoptim.md             |  5 +-
 docs/src/api/minpack.md                       |  5 +-
 docs/src/api/nlsolve.md                       |  5 +-
 docs/src/api/nlsolvers.md                     |  5 +-
 docs/src/api/petsc.md                         |  5 +-
 docs/src/api/siamfanlequations.md             |  5 +-
 docs/src/api/speedmapping.md                  |  5 +-
 docs/src/api/sundials.md                      |  4 +-
 docs/src/basics/diagnostics_api.md            | 18 ++--
 docs/src/basics/faq.md                        | 47 ++++-----
 docs/src/native/steadystatediffeq.md          |  4 +-
 docs/src/solvers/steady_state_solvers.md      | 10 +-
 docs/src/tutorials/code_optimization.md       | 32 +++----
 docs/src/tutorials/getting_started.md         | 40 ++++----
 docs/src/tutorials/iterator_interface.md      |  8 +-
 docs/src/tutorials/large_systems.md           | 95 ++++++++++---------
 docs/src/tutorials/modelingtoolkit.md         | 35 +++----
 docs/src/tutorials/nonlinear_solve_gpus.md    | 49 +++++-----
 .../tutorials/optimizing_parameterized_ode.md | 33 +++----
 docs/src/tutorials/small_compile.md           | 14 +--
 docs/src/tutorials/snes_ex2.md                | 30 +++---
 25 files changed, 251 insertions(+), 220 deletions(-)

diff --git a/docs/src/api/SciMLJacobianOperators.md b/docs/src/api/SciMLJacobianOperators.md
index 91b7f61df..79b576ab8 100644
--- a/docs/src/api/SciMLJacobianOperators.md
+++ b/docs/src/api/SciMLJacobianOperators.md
@@ -10,7 +10,7 @@ operator built on top on DifferentiationInterface.jl.
 ```julia
 import Pkg
 Pkg.add("SciMLJacobianOperators")
-using SciMLJacobianOperators
+import SciMLJacobianOperators
 ```
 
 ## Jacobian API
diff --git a/docs/src/api/fastlevenbergmarquardt.md b/docs/src/api/fastlevenbergmarquardt.md
index 423c70bbc..c8f37f291 100644
--- a/docs/src/api/fastlevenbergmarquardt.md
+++ b/docs/src/api/fastlevenbergmarquardt.md
@@ -5,9 +5,10 @@ interface. Note that these solvers do not come by default, and thus one needs to
 the package before using these solvers:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("FastLevenbergMarquardt")
-using FastLevenbergMarquardt, NonlinearSolve
+import FastLevenbergMarquardt
+import NonlinearSolve as NLS
 ```
 
 ## Solver API
diff --git a/docs/src/api/fixedpointacceleration.md b/docs/src/api/fixedpointacceleration.md
index 38e6ae0b5..178f23c4f 100644
--- a/docs/src/api/fixedpointacceleration.md
+++ b/docs/src/api/fixedpointacceleration.md
@@ -5,9 +5,10 @@ interface. Note that these solvers do not come by default, and thus one needs to
 the package before using these solvers:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("FixedPointAcceleration")
-using FixedPointAcceleration, NonlinearSolve
+import FixedPointAcceleration
+import NonlinearSolve as NLS
 ```
 
 ## Solver API
diff --git a/docs/src/api/homotopycontinuation.md b/docs/src/api/homotopycontinuation.md
index 783eb6883..eabacaabd 100644
--- a/docs/src/api/homotopycontinuation.md
+++ b/docs/src/api/homotopycontinuation.md
@@ -5,9 +5,10 @@ HomotopyContinuation.jl. This solver is not included by default and needs
 to be installed separately:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("NonlinearSolveHomotopyContinuation")
-using NonlinearSolveHomotopyContinuation, NonlinearSolve
+import NonlinearSolveHomotopyContinuation
+import NonlinearSolve as NLS
 ```
 
 # Solver API
diff --git a/docs/src/api/leastsquaresoptim.md b/docs/src/api/leastsquaresoptim.md
index 0581ee7c2..d4f127628 100644
--- a/docs/src/api/leastsquaresoptim.md
+++ b/docs/src/api/leastsquaresoptim.md
@@ -5,9 +5,10 @@ interface. Note that these solvers do not come by default, and thus one needs to
 the package before using these solvers:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("LeastSquaresOptim")
-using LeastSquaresOptim, NonlinearSolve
+import LeastSquaresOptim
+import NonlinearSolve as NLS
 ```
 
 ## Solver API
diff --git a/docs/src/api/minpack.md b/docs/src/api/minpack.md
index d55491c0e..c4f200ace 100644
--- a/docs/src/api/minpack.md
+++ b/docs/src/api/minpack.md
@@ -5,9 +5,10 @@ these solvers do not come by default, and thus one needs to install the package
 these solvers:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("MINPACK")
-using MINPACK, NonlinearSolve
+import MINPACK
+import NonlinearSolve as NLS
 ```
 
 ## Solver API
diff --git a/docs/src/api/nlsolve.md b/docs/src/api/nlsolve.md
index bea379953..a4cec0ef5 100644
--- a/docs/src/api/nlsolve.md
+++ b/docs/src/api/nlsolve.md
@@ -5,9 +5,10 @@ that these solvers do not come by default, and thus one needs to install the pac
 using these solvers:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("NLsolve")
-using NLsolve, NonlinearSolve
+import NLsolve
+import NonlinearSolve as NLS
 ```
 
 ## Solver API
diff --git a/docs/src/api/nlsolvers.md b/docs/src/api/nlsolvers.md
index a1cc72656..deb5afeb2 100644
--- a/docs/src/api/nlsolvers.md
+++ b/docs/src/api/nlsolvers.md
@@ -5,9 +5,10 @@ that these solvers do not come by default, and thus one needs to install the pac
 using these solvers:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("NLSolvers")
-using NLSolvers, NonlinearSolve
+import NLSolvers
+import NonlinearSolve as NLS
 ```
 
 ## Solver API
diff --git a/docs/src/api/petsc.md b/docs/src/api/petsc.md
index 7593ac5dd..05ac23827 100644
--- a/docs/src/api/petsc.md
+++ b/docs/src/api/petsc.md
@@ -5,9 +5,10 @@ that these solvers do not come by default, and thus one needs to install the pac
 using these solvers:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("PETSc")
-using PETSc, NonlinearSolve
+import PETSc
+import NonlinearSolve as NLS
 ```
 
 ## Solver API
diff --git a/docs/src/api/siamfanlequations.md b/docs/src/api/siamfanlequations.md
index 5ee36ba12..d1b4db535 100644
--- a/docs/src/api/siamfanlequations.md
+++ b/docs/src/api/siamfanlequations.md
@@ -6,9 +6,10 @@ interface. Note that these solvers do not come by default, and thus one needs to
 the package before using these solvers:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("SIAMFANLEquations")
-using SIAMFANLEquations, NonlinearSolve
+import SIAMFANLEquations
+import NonlinearSolve as NLS
 ```
 
 ## Solver API
diff --git a/docs/src/api/speedmapping.md b/docs/src/api/speedmapping.md
index 76b13ba18..3a99801fc 100644
--- a/docs/src/api/speedmapping.md
+++ b/docs/src/api/speedmapping.md
@@ -5,9 +5,10 @@ interface. Note that these solvers do not come by default, and thus one needs to
 the package before using these solvers:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("SpeedMapping")
-using SpeedMapping, NonlinearSolve
+import SpeedMapping
+import NonlinearSolve as NLS
 ```
 
 ## Solver API
diff --git a/docs/src/api/sundials.md b/docs/src/api/sundials.md
index 80a6b367d..82f6f0302 100644
--- a/docs/src/api/sundials.md
+++ b/docs/src/api/sundials.md
@@ -5,9 +5,9 @@ Note that these solvers do not come by default, and thus one needs to install th
 before using these solvers:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("Sundials")
-using Sundials
+import Sundials
 ```
 
 These methods can be used independently of the rest of NonlinearSolve.jl.
diff --git a/docs/src/basics/diagnostics_api.md b/docs/src/basics/diagnostics_api.md
index b3f8dea4c..716c7faeb 100644
--- a/docs/src/basics/diagnostics_api.md
+++ b/docs/src/basics/diagnostics_api.md
@@ -27,7 +27,7 @@ Note that you will have to restart Julia to disable the timer outputs once enabl
 ## Example Usage
 
 ```@example diagnostics_example
-using NonlinearSolve
+import NonlinearSolve as NLS
 
 function nlfunc(resid, u0, p)
     resid[1] = u0[1] * u0[1] - p
@@ -36,23 +36,23 @@ function nlfunc(resid, u0, p)
     nothing
 end
 
-prob = NonlinearProblem(nlfunc, [1.0, 3.0, 5.0], 2.0)
+prob = NLS.NonlinearProblem(nlfunc, [1.0, 3.0, 5.0], 2.0)
 
-solve(prob)
+NLS.solve(prob)
 ```
 
 This produced the output, but it is hard to diagnose what is going on. We can turn on
 the trace to see what is happening:
 
 ```@example diagnostics_example
-solve(prob; show_trace = Val(true), trace_level = TraceAll(10))
+NLS.solve(prob; show_trace = Val(true), trace_level = NLS.TraceAll(10))
 nothing; # hide
 ```
 
 You can also store the trace in the solution object:
 
 ```@example diagnostics_example
-sol = solve(prob; trace_level = TraceAll(), store_trace = Val(true));
+sol = NLS.solve(prob; trace_level = NLS.TraceAll(), store_trace = Val(true));
 
 sol.trace
 ```
@@ -62,16 +62,16 @@ to use the `init` and `solve!` API for this. The `TimerOutput` will be present i
 `cache.timer`. However, note that for poly-algorithms this is currently not implemented.
 
 ```@example diagnostics_example
-cache = init(prob, NewtonRaphson(); show_trace = Val(true));
-solve!(cache)
+cache = NLS.init(prob, NLS.NewtonRaphson(); show_trace = Val(true));
+NLS.solve!(cache)
 cache.timer
 ```
 
 Let's try for some other solver:
 
 ```@example diagnostics_example
-cache = init(prob, DFSane(); show_trace = Val(true), trace_level = TraceMinimal(50));
-solve!(cache)
+cache = NLS.init(prob, NLS.DFSane(); show_trace = Val(true), trace_level = NLS.TraceMinimal(50));
+NLS.solve!(cache)
 cache.timer
 ```
 
diff --git a/docs/src/basics/faq.md b/docs/src/basics/faq.md
index 8411477fe..4e37612a4 100644
--- a/docs/src/basics/faq.md
+++ b/docs/src/basics/faq.md
@@ -6,7 +6,8 @@ This is addressed in a [Twitter thread with the author of the improved fzero](ht
 On the test example:
 
 ```@example
-using NonlinearSolve, BenchmarkTools
+import NonlinearSolve as NLS
+import BenchmarkTools
 
 const N = 100_000;
 levels = 1.5 .* rand(N);
@@ -15,23 +16,23 @@ myfun(x, lv) = x * sin(x) - lv
 
 function f(out, levels, u0)
     for i in 1:N
-        out[i] = solve(
-            IntervalNonlinearProblem{false}(
-                IntervalNonlinearFunction{false}(myfun), u0, levels[i]),
-            Falsi()).u
+        out[i] = NLS.solve(
+            NLS.IntervalNonlinearProblem{false}(
+                NLS.IntervalNonlinearFunction{false}(myfun), u0, levels[i]),
+            NLS.Falsi()).u
     end
 end
 
 function f2(out, levels, u0)
     for i in 1:N
-        out[i] = solve(
-            NonlinearProblem{false}(NonlinearFunction{false}(myfun), u0, levels[i]),
-            SimpleNewtonRaphson()).u
+        out[i] = NLS.solve(
+            NLS.NonlinearProblem{false}(NLS.NonlinearFunction{false}(myfun), u0, levels[i]),
+            NLS.SimpleNewtonRaphson()).u
     end
 end
 
-@btime f(out, levels, (0.0, 2.0))
-@btime f2(out, levels, 1.0)
+BenchmarkTools.@btime f(out, levels, (0.0, 2.0))
+BenchmarkTools.@btime f2(out, levels, 1.0)
 ```
 
 MATLAB 2022a achieves 1.66s. Try this code yourself: we receive 0.009 seconds, or a 184x
@@ -46,7 +47,8 @@ input types. For example, consider this example taken from
 [this issue](https://github.com/SciML/NonlinearSolve.jl/issues/298)
 
 ```@example dual_error_faq
-using NonlinearSolve, Random
+import NonlinearSolve as NLS
+import Random
 
 function fff_incorrect(var, p)
     v_true = [1.0, 0.1, 2.0, 0.5]
@@ -58,9 +60,9 @@ end
 v_true = [1.0, 0.1, 2.0, 0.5]
 v_init = v_true .+ randn!(similar(v_true)) * 0.1
 
-prob_oop = NonlinearLeastSquaresProblem{false}(fff_incorrect, v_init)
+prob_oop = NLS.NonlinearLeastSquaresProblem{false}(fff_incorrect, v_init)
 try
-    sol = solve(prob_oop, LevenbergMarquardt(); maxiters = 10000, abstol = 1e-8)
+    sol = NLS.solve(prob_oop, NLS.LevenbergMarquardt(); maxiters = 10000, abstol = 1e-8)
 catch e
     @error e
 end
@@ -91,8 +93,8 @@ be a Dual number. This causes the error. To fix it:
         return xx - v_true
     end
     
-    prob_oop = NonlinearLeastSquaresProblem{false}(fff_correct, v_init)
-    sol = solve(prob_oop, LevenbergMarquardt(); maxiters = 10000, abstol = 1e-8)
+    prob_oop = NLS.NonlinearLeastSquaresProblem{false}(fff_correct, v_init)
+    sol = NLS.solve(prob_oop, NLS.LevenbergMarquardt(); maxiters = 10000, abstol = 1e-8)
     ```
 
 ## I thought NonlinearSolve.jl was type-stable and fast. But it isn't, why?
@@ -106,13 +108,14 @@ static arrays, ForwardDiff will create type unstable code and lead to dynamic di
 internally. See this simple example:
 
 ```@example type_unstable
-using NonlinearSolve, InteractiveUtils
+import NonlinearSolve as NLS
+import InteractiveUtils
 
 f(u, p) = @. u^2 - p
 
-prob = NonlinearProblem{false}(f, 1.0, 2.0)
+prob = NLS.NonlinearProblem{false}(f, 1.0, 2.0)
 
-@code_warntype solve(prob, NewtonRaphson())
+InteractiveUtils.@code_warntype NLS.solve(prob, NLS.NewtonRaphson())
 nothing # hide
 ```
 
@@ -120,11 +123,11 @@ Notice that this was type-stable, since it is a scalar problem. Now what happens
 arrays
 
 ```@example type_unstable
-using StaticArrays
+import StaticArrays
 
-prob = NonlinearProblem{false}(f, @SVector([1.0, 2.0]), 2.0)
+prob = NLS.NonlinearProblem{false}(f, StaticArrays.@SVector([1.0, 2.0]), 2.0)
 
-@code_warntype solve(prob, NewtonRaphson())
+InteractiveUtils.@code_warntype NLS.solve(prob, NLS.NewtonRaphson())
 nothing # hide
 ```
 
@@ -133,7 +136,7 @@ Again Type-Stable! Now let's try using a regular array:
 ```@example type_unstable
 prob = NonlinearProblem(f, [1.0, 2.0], 2.0)
 
-@code_warntype solve(prob, NewtonRaphson())
+InteractiveUtils.@code_warntype NLS.solve(prob, NLS.NewtonRaphson())
 nothing # hide
 ```
 
diff --git a/docs/src/native/steadystatediffeq.md b/docs/src/native/steadystatediffeq.md
index 471fc3f01..cb0cb4741 100644
--- a/docs/src/native/steadystatediffeq.md
+++ b/docs/src/native/steadystatediffeq.md
@@ -6,9 +6,9 @@ interface. Note that these solvers do not come by default, and thus one needs to
 the package before using these solvers:
 
 ```julia
-using Pkg
+import Pkg
 Pkg.add("SteadyStateDiffEq")
-using SteadyStateDiffEq
+import SteadyStateDiffEq as SSDE
 ```
 
 These methods can be used independently of the rest of NonlinearSolve.jl
diff --git a/docs/src/solvers/steady_state_solvers.md b/docs/src/solvers/steady_state_solvers.md
index 91530c448..3dddddd33 100644
--- a/docs/src/solvers/steady_state_solvers.md
+++ b/docs/src/solvers/steady_state_solvers.md
@@ -56,11 +56,13 @@ though often computationally more expensive than direct methods.
 Example usage:
 
 ```julia
-using NonlinearSolve, SteadyStateDiffEq, OrdinaryDiffEq
-sol = solve(prob, DynamicSS(Tsit5()))
+import NonlinearSolve as NLS
+import SteadyStateDiffEq as SSDE
+import OrdinaryDiffEq as ODE
+sol = NLS.solve(prob, SSDE.DynamicSS(ODE.Tsit5()))
 
-using Sundials
-sol = solve(prob, DynamicSS(CVODE_BDF()), dt = 1.0)
+import Sundials
+sol = NLS.solve(prob, SSDE.DynamicSS(Sundials.CVODE_BDF()), dt = 1.0)
 ```
 
 !!! note
diff --git a/docs/src/tutorials/code_optimization.md b/docs/src/tutorials/code_optimization.md
index b7eb9c174..1602e3683 100644
--- a/docs/src/tutorials/code_optimization.md
+++ b/docs/src/tutorials/code_optimization.md
@@ -31,21 +31,21 @@ systems.
 Take for example a prototypical small nonlinear solver code in its out-of-place form:
 
 ```@example small_opt
-using NonlinearSolve
+import NonlinearSolve as NLS
 
 f(u, p) = u .* u .- p
 u0 = [1.0, 1.0]
 p = 2.0
-prob = NonlinearProblem(f, u0, p)
-sol = solve(prob, NewtonRaphson())
+prob = NLS.NonlinearProblem(f, u0, p)
+sol = NLS.solve(prob, NLS.NewtonRaphson())
 ```
 
 We can use BenchmarkTools.jl to get more precise timings:
 
 ```@example small_opt
-using BenchmarkTools
+import BenchmarkTools
 
-@benchmark solve(prob, NewtonRaphson())
+BenchmarkTools.@benchmark NLS.solve(prob, NLS.NewtonRaphson())
 ```
 
 Note that this way of writing the function is a shorthand for:
@@ -70,8 +70,8 @@ function f(du, u, p)
     return nothing
 end
 
-prob = NonlinearProblem(f, u0, p)
-@benchmark sol = solve(prob, NewtonRaphson())
+prob = NLS.NonlinearProblem(f, u0, p)
+@benchmark sol = NLS.solve(prob, NLS.NewtonRaphson())
 ```
 
 Notice how much faster this already runs! We can make this code even simpler by using
@@ -83,7 +83,7 @@ function f(du, u, p)
     return nothing
 end
 
-@benchmark sol = solve(prob, NewtonRaphson())
+@benchmark sol = NLS.solve(prob, NLS.NewtonRaphson())
 ```
 
 ## Further Optimizations for Small Nonlinear Solves with Static Arrays and SimpleNonlinearSolve
@@ -110,9 +110,9 @@ arrays have their length determined at compile-time. They are created using macr
 to normal array expressions, for example:
 
 ```@example small_opt
-using StaticArrays
+import StaticArrays
 
-A = SA[2.0, 3.0, 5.0]
+A = StaticArrays.SA[2.0, 3.0, 5.0]
 typeof(A)
 ```
 
@@ -135,20 +135,20 @@ array. Doing it with broadcasting looks like:
 ```@example small_opt
 f_SA(u, p) = u .* u .- p
 
-u0 = SA[1.0, 1.0]
+u0 = StaticArrays.SA[1.0, 1.0]
 p = 2.0
-prob = NonlinearProblem(f_SA, u0, p)
+prob = NLS.NonlinearProblem(f_SA, u0, p)
 
-@benchmark solve(prob, NewtonRaphson())
+BenchmarkTools.@benchmark NLS.solve(prob, NLS.NewtonRaphson())
 ```
 
 Note that only change here is that `u0` is made into a StaticArray! If we needed to write
 `f` out for a more complex nonlinear case, then we'd simply do the following:
 
 ```@example small_opt
-f_SA(u, p) = SA[u[1] * u[1] - p, u[2] * u[2] - p]
+f_SA(u, p) = StaticArrays.SA[u[1] * u[1] - p, u[2] * u[2] - p]
 
-@benchmark solve(prob, NewtonRaphson())
+BenchmarkTools.@benchmark NLS.solve(prob, NLS.NewtonRaphson())
 ```
 
 However, notice that this did not give us a speedup but rather a slowdown. This is because
@@ -159,7 +159,7 @@ of the methods from SimpleNonlinearSolve.jl which are designed for these small-s
 examples. Let's now use `SimpleNewtonRaphson`:
 
 ```@example small_opt
-@benchmark solve(prob, SimpleNewtonRaphson())
+BenchmarkTools.@benchmark NLS.solve(prob, NLS.SimpleNewtonRaphson())
 ```
 
 And there we go, around `40ns` from our starting point of almost `4μs`!
diff --git a/docs/src/tutorials/getting_started.md b/docs/src/tutorials/getting_started.md
index d0729eeaf..67d406bdf 100644
--- a/docs/src/tutorials/getting_started.md
+++ b/docs/src/tutorials/getting_started.md
@@ -31,13 +31,13 @@ the parameters of the system. For example, the following solves the vector
 equation $$f(u) = u^2 - p$$ for a vector of equations:
 
 ```@example 1
-using NonlinearSolve
+import NonlinearSolve as NLS
 
 f(u, p) = u .* u .- p
 u0 = [1.0, 1.0]
 p = 2.0
-prob = NonlinearProblem(f, u0, p)
-sol = solve(prob)
+prob = NLS.NonlinearProblem(f, u0, p)
+sol = NLS.solve(prob)
 ```
 
 where `u0` is the initial condition for the rootfinder. Native NonlinearSolve.jl solvers use
@@ -77,6 +77,7 @@ the general command `SciMLBase.successful_retcode` to check whether the solution
 exited as intended:
 
 ```@example 1
+import SciMLBase
 SciMLBase.successful_retcode(sol)
 ```
 
@@ -92,13 +93,13 @@ For more information on `NonlinearSolution`s, see the
 
 ### Interacting with the Solver Options
 
-While `sol = solve(prob)` worked for our case here, in many situations you may need to
+While `sol = NonlinearSolve.solve(prob)` worked for our case here, in many situations you may need to
 interact more deeply with how the solving process is done. First things first, you can
 specify the solver using the positional arguments. For example, let's set the solver to
 `TrustRegion()`:
 
 ```@example 1
-solve(prob, TrustRegion())
+NLS.solve(prob, NLS.TrustRegion())
 ```
 
 For a complete list of solver choices, see
@@ -108,7 +109,7 @@ Next we can modify the tolerances. Here let's set some really low tolerances to
 tight solution:
 
 ```@example 1
-solve(prob, TrustRegion(), reltol = 1e-12, abstol = 1e-12)
+NLS.solve(prob, NLS.TrustRegion(), reltol = 1e-12, abstol = 1e-12)
 ```
 
 There are many more options for doing this configuring. Specifically for handling
@@ -124,11 +125,11 @@ condition, you pass a `uspan (a,b)` bracket in which the zero is expected to liv
 example:
 
 ```@example 1
-using NonlinearSolve
+import NonlinearSolve as NLS
 f(u, p) = u * u - 2.0
 uspan = (1.0, 2.0) # brackets
-prob_int = IntervalNonlinearProblem(f, uspan)
-sol = solve(prob_int)
+prob_int = NLS.IntervalNonlinearProblem(f, uspan)
+sol = NLS.solve(prob_int)
 ```
 
 All of the same option handling from before works just as before, now just with different
@@ -136,7 +137,7 @@ solver choices (see the [bracketing solvers](@ref bracketing) page for more deta
 example, let's set the solver to `ITP()` and set a high absolute tolerance:
 
 ```@example 1
-sol = solve(prob_int, ITP(), abstol = 0.01)
+sol = NLS.solve(prob_int, NLS.ITP(), abstol = 0.01)
 ```
 
 ## Problem Type 3: Solving Steady State Problems
@@ -146,26 +147,27 @@ For Steady State Problems, we have a wrapper package
 automates handling SteadyStateProblems with NonlinearSolve and OrdinaryDiffEq.
 
 ```@example 1
-using NonlinearSolve, SteadyStateDiffEq
+import NonlinearSolve as NLS
+import SteadyStateDiffEq as SSDE
 
 f(u, p, t) = [2 - 2u[1]; u[1] - 4u[2]]
 u0 = [0.0, 0.0]
-prob = SteadyStateProblem(f, u0)
+prob = SSDE.SteadyStateProblem(f, u0)
 
-solve(prob, SSRootfind())
+NLS.solve(prob, SSDE.SSRootfind())
 ```
 
 If you don't provide a nonlinear solver to `SSRootfind` it uses a polyalgorithm from
 NonlinearSolve. We can also provide the actual nonlinear solver to use:
 
 ```@example 1
-solve(prob, SSRootfind(Broyden()))
+NLS.solve(prob, SSDE.SSRootfind(NLS.Broyden()))
 ```
 
 ## Problem Type 4: Solving Nonlinear Least Squares Problems
 
 ```@example 1
-using NonlinearSolve
+import NonlinearSolve as NLS
 
 function nlls!(du, u, p)
     du[1] = 2u[1] - 2
@@ -180,16 +182,16 @@ be skipped for out of place problems):
 
 ```@example 1
 u0 = [0.0, 0.0]
-prob = NonlinearLeastSquaresProblem(
-    NonlinearFunction(nlls!, resid_prototype = zeros(3)), u0)
+prob = NLS.NonlinearLeastSquaresProblem(
+    NLS.NonlinearFunction(nlls!, resid_prototype = zeros(3)), u0)
 
-solve(prob)
+NLS.solve(prob)
 ```
 
 Same as before, we can change the solver and tolerances:
 
 ```@example 1
-solve(prob, GaussNewton(), reltol = 1e-12, abstol = 1e-12)
+NLS.solve(prob, NLS.GaussNewton(), reltol = 1e-12, abstol = 1e-12)
 ```
 
 ## Going Beyond the Basics: How to use the Documentation
diff --git a/docs/src/tutorials/iterator_interface.md b/docs/src/tutorials/iterator_interface.md
index 26873fe35..44fb9a234 100644
--- a/docs/src/tutorials/iterator_interface.md
+++ b/docs/src/tutorials/iterator_interface.md
@@ -4,13 +4,13 @@ There is an iterator form of the nonlinear solver which somewhat mirrors the Dif
 integrator interface:
 
 ```@example iterator_interface
-using NonlinearSolve
+import NonlinearSolve as NLS
 
 f(u, p) = u .* u .- 2.0
 u0 = 1.5
-probB = NonlinearProblem(f, u0)
+probB = NLS.NonlinearProblem(f, u0)
 
-nlcache = init(probB, NewtonRaphson())
+nlcache = NLS.init(probB, NLS.NewtonRaphson())
 ```
 
 `init` takes the same keyword arguments as [`solve`](@ref solver_options), but it returns a
@@ -27,7 +27,7 @@ We can perform 10 steps of the Newton-Raphson solver with the following:
 
 ```@example iterator_interface
 for i in 1:10
-    step!(nlcache)
+    NLS.step!(nlcache)
 end
 ```
 
diff --git a/docs/src/tutorials/large_systems.md b/docs/src/tutorials/large_systems.md
index 434acb185..86cf6cbd0 100644
--- a/docs/src/tutorials/large_systems.md
+++ b/docs/src/tutorials/large_systems.md
@@ -60,7 +60,10 @@ broadcast). Use `dx=dy=1/32`.
 The resulting `NonlinearProblem` definition is:
 
 ```@example ill_conditioned_nlprob
-using NonlinearSolve, LinearAlgebra, SparseArrays, LinearSolve
+import NonlinearSolve as NLS
+import LinearAlgebra
+import SparseArrays
+import LinearSolve as LS
 
 const N = 32
 const xyd_brusselator = range(0, stop = 1, length = N)
@@ -99,7 +102,7 @@ function init_brusselator_2d(xyd)
 end
 
 u0 = init_brusselator_2d(xyd_brusselator)
-prob_brusselator_2d = NonlinearProblem(
+prob_brusselator_2d = NLS.NonlinearProblem(
     brusselator_2d_loop, u0, p; abstol = 1e-10, reltol = 1e-10
 )
 ```
@@ -130,28 +133,28 @@ In the next section, we will show how to specify `sparsity` to trigger automatic
 detection.
 
 ```@example ill_conditioned_nlprob
-using BenchmarkTools # for @btime
+import BenchmarkTools # for @btime
 
-@btime solve(prob_brusselator_2d, NewtonRaphson());
+BenchmarkTools.BenchmarkTools.@btime NLS.solve(prob_brusselator_2d, NLS.NewtonRaphson());
 nothing # hide
 ```
 
 ```@example ill_conditioned_nlprob
-using SparseConnectivityTracer
+import SparseConnectivityTracer
 
-prob_brusselator_2d_autosparse = NonlinearProblem(
-    NonlinearFunction(brusselator_2d_loop; sparsity = TracerSparsityDetector()),
+prob_brusselator_2d_autosparse = NLS.NonlinearProblem(
+    NLS.NonlinearFunction(brusselator_2d_loop; sparsity = SparseConnectivityTracer.TracerSparsityDetector()),
     u0, p; abstol = 1e-10, reltol = 1e-10
 )
 
-@btime solve(prob_brusselator_2d_autosparse,
-    NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 12)));
-@btime solve(prob_brusselator_2d_autosparse,
-    NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 12),
-        linsolve = KLUFactorization()));
-@btime solve(prob_brusselator_2d_autosparse,
-    NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 12),
-        linsolve = KrylovJL_GMRES()));
+BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_autosparse,
+    NLS.NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 12)));
+BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_autosparse,
+    NLS.NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 12),
+        linsolve = LS.KLUFactorization()));
+BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_autosparse,
+    NLS.NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 12),
+        linsolve = LS.KrylovJL_GMRES()));
 nothing # hide
 ```
 
@@ -176,38 +179,39 @@ arguments, and it will kick out a sparse matrix with our pattern, that we can tu
     actual implementation of sparsity detection.
 
 ```@example ill_conditioned_nlprob
-using SparseConnectivityTracer, ADTypes
+import SparseConnectivityTracer
+import ADTypes
 
 f! = (du, u) -> brusselator_2d_loop(du, u, p)
 du0 = similar(u0)
-jac_sparsity = ADTypes.jacobian_sparsity(f!, du0, u0, TracerSparsityDetector())
+jac_sparsity = ADTypes.jacobian_sparsity(f!, du0, u0, SparseConnectivityTracer.TracerSparsityDetector())
 ```
 
 Notice that Julia gives a nice print out of the sparsity pattern. That's neat, and would be
 tedious to build by hand! Now we just pass it to the `NonlinearFunction` like as before:
 
 ```@example ill_conditioned_nlprob
-ff = NonlinearFunction(brusselator_2d_loop; jac_prototype = jac_sparsity)
+ff = NLS.NonlinearFunction(brusselator_2d_loop; jac_prototype = jac_sparsity)
 ```
 
 Build the `NonlinearProblem`:
 
 ```@example ill_conditioned_nlprob
-prob_brusselator_2d_sparse = NonlinearProblem(ff, u0, p; abstol = 1e-10, reltol = 1e-10)
+prob_brusselator_2d_sparse = NLS.NonlinearProblem(ff, u0, p; abstol = 1e-10, reltol = 1e-10)
 ```
 
 Now let's see how the version with sparsity compares to the version without:
 
 ```@example ill_conditioned_nlprob
-@btime solve(prob_brusselator_2d, NewtonRaphson());
-@btime solve(prob_brusselator_2d_sparse, NewtonRaphson());
-@btime solve(prob_brusselator_2d_sparse, NewtonRaphson(linsolve = KLUFactorization()));
+BenchmarkTools.@btime NLS.solve(prob_brusselator_2d, NLS.NewtonRaphson());
+BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_sparse, NLS.NewtonRaphson());
+BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_sparse, NLS.NewtonRaphson(linsolve = LS.KLUFactorization()));
 nothing # hide
 ```
 
 Note that depending on the properties of the sparsity pattern, one may want to try
-alternative linear solvers such as `NewtonRaphson(linsolve = KLUFactorization())`
-or `NewtonRaphson(linsolve = UMFPACKFactorization())`
+alternative linear solvers such as `NLS.NewtonRaphson(linsolve = LS.KLUFactorization())`
+or `NLS.NewtonRaphson(linsolve = LS.UMFPACKFactorization())`
 
 ## Using Jacobian-Free Newton-Krylov
 
@@ -217,7 +221,7 @@ Krylov method. To swap the linear solver out, we use the `linsolve` command and
 GMRES linear solver.
 
 ```@example ill_conditioned_nlprob
-@btime solve(prob_brusselator_2d, NewtonRaphson(linsolve = KrylovJL_GMRES()));
+BenchmarkTools.@btime NLS.solve(prob_brusselator_2d, NLS.NewtonRaphson(linsolve = LS.KrylovJL_GMRES()));
 nothing # hide
 ```
 
@@ -244,12 +248,12 @@ used in the solution of the ODE. An example of this with using
 
 ```julia
 # FIXME: On 1.10+ this is broken. Skipping this for now.
-using IncompleteLU
+import IncompleteLU
 
-incompletelu(W, p = nothing) = ilu(W, τ = 50.0), LinearAlgebra.I
+incompletelu(W, p = nothing) = IncompleteLU.ilu(W, τ = 50.0), LinearAlgebra.I
 
-@btime solve(prob_brusselator_2d_sparse,
-    NewtonRaphson(linsolve = KrylovJL_GMRES(precs = incompletelu), concrete_jac = true)
+BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_sparse,
+    NLS.NewtonRaphson(linsolve = LS.KrylovJL_GMRES(precs = incompletelu), concrete_jac = true)
 );
 nothing # hide
 ```
@@ -271,15 +275,15 @@ parameter. Another option is to use
 which is more automatic. The setup is very similar to before:
 
 ```@example ill_conditioned_nlprob
-using AlgebraicMultigrid
+import AlgebraicMultigrid
 
 function algebraicmultigrid(W, p = nothing)
-    return aspreconditioner(ruge_stuben(convert(AbstractMatrix, W))), LinearAlgebra.I
+    return AlgebraicMultigrid.aspreconditioner(AlgebraicMultigrid.ruge_stuben(convert(AbstractMatrix, W))), LinearAlgebra.I
 end
 
-@btime solve(prob_brusselator_2d_sparse,
-    NewtonRaphson(
-        linsolve = KrylovJL_GMRES(; precs = algebraicmultigrid), concrete_jac = true
+BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_sparse,
+    NLS.NewtonRaphson(
+        linsolve = LS.KrylovJL_GMRES(; precs = algebraicmultigrid), concrete_jac = true
     )
 );
 nothing # hide
@@ -297,10 +301,10 @@ function algebraicmultigrid2(W, p = nothing)
     return Pl, LinearAlgebra.I
 end
 
-@btime solve(
+BenchmarkTools.@btime NLS.solve(
     prob_brusselator_2d_sparse,
-    NewtonRaphson(
-        linsolve = KrylovJL_GMRES(precs = algebraicmultigrid2), concrete_jac = true
+    NLS.NewtonRaphson(
+        linsolve = LS.KrylovJL_GMRES(precs = algebraicmultigrid2), concrete_jac = true
     )
 );
 nothing # hide
@@ -313,18 +317,19 @@ for the exact sparsity detection case, we left out the time it takes to perform
 sparsity detection. Let's compare the two by setting the sparsity detection algorithms.
 
 ```@example ill_conditioned_nlprob
-using DifferentiationInterface, SparseConnectivityTracer
+import DifferentiationInterface
+import SparseConnectivityTracer
 
-prob_brusselator_2d_exact_tracer = NonlinearProblem(
-    NonlinearFunction(brusselator_2d_loop; sparsity = TracerSparsityDetector()),
+prob_brusselator_2d_exact_tracer = NLS.NonlinearProblem(
+    NLS.NonlinearFunction(brusselator_2d_loop; sparsity = SparseConnectivityTracer.TracerSparsityDetector()),
     u0, p; abstol = 1e-10, reltol = 1e-10)
-prob_brusselator_2d_approx_di = NonlinearProblem(
-    NonlinearFunction(brusselator_2d_loop;
-        sparsity = DenseSparsityDetector(AutoForwardDiff(); atol = 1e-4)),
+prob_brusselator_2d_approx_di = NLS.NonlinearProblem(
+    NLS.NonlinearFunction(brusselator_2d_loop;
+        sparsity = DifferentiationInterface.DenseSparsityDetector(DifferentiationInterface.AutoForwardDiff(); atol = 1e-4)),
     u0, p; abstol = 1e-10, reltol = 1e-10)
 
-@btime solve(prob_brusselator_2d_exact_tracer, NewtonRaphson());
-@btime solve(prob_brusselator_2d_approx_di, NewtonRaphson());
+BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_exact_tracer, NLS.NewtonRaphson());
+BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_approx_di, NLS.NewtonRaphson());
 nothing # hide
 ```
 
diff --git a/docs/src/tutorials/modelingtoolkit.md b/docs/src/tutorials/modelingtoolkit.md
index e2016dde0..fde4f1033 100644
--- a/docs/src/tutorials/modelingtoolkit.md
+++ b/docs/src/tutorials/modelingtoolkit.md
@@ -6,14 +6,15 @@ for the numerical solvers which can make it easy to symbolically modify and gene
 equations to be solved. The basic form of using ModelingToolkit looks as follows:
 
 ```julia
-using ModelingToolkit, NonlinearSolve
+import ModelingToolkit as MTK
+import NonlinearSolve as NLS
 
-@variables x y z
-@parameters σ ρ β
+MTK.@variables x y z
+MTK.@parameters σ ρ β
 
 # Define a nonlinear system
 eqs = [0 ~ σ * (y - x), 0 ~ x * (ρ - z) - y, 0 ~ x * y - β * z]
-@mtkbuild ns = NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
+MTK.@mtkbuild ns = MTK.NonlinearSystem(eqs, [x, y, z], [σ, ρ, β])
 
 u0 = [x => 1.0, y => 0.0, z => 0.0]
 
@@ -21,8 +22,8 @@ ps = [σ => 10.0
       ρ => 26.0
       β => 8 / 3]
 
-prob = NonlinearProblem(ns, u0, ps)
-sol = solve(prob, NewtonRaphson())
+prob = NLS.NonlinearProblem(ns, u0, ps)
+sol = NLS.solve(prob, NLS.NewtonRaphson())
 ```
 
 ## Symbolic Derivations of Extra Functions
@@ -31,21 +32,21 @@ As a symbolic system, ModelingToolkit can be used to represent the equations and
 forms. For example, let's look at the equations:
 
 ```julia
-equations(ns)
+MTK.equations(ns)
 ```
 
 We can ask it what the Jacobian of our system is via `calculate_jacobian`:
 
 ```julia
-calculate_jacobian(ns)
+MTK.calculate_jacobian(ns)
 ```
 
 We can tell MTK to generate a computable form of this analytical Jacobian via `jac = true`
 to help the solver use efficient forms:
 
 ```julia
-prob = NonlinearProblem(ns, u0, ps, jac = true)
-sol = solve(prob, NewtonRaphson())
+prob = NLS.NonlinearProblem(ns, u0, ps, jac = true)
+sol = NLS.solve(prob, NLS.NewtonRaphson())
 ```
 
 ## Symbolic Simplification of Nonlinear Systems via Tearing
@@ -55,34 +56,34 @@ the systems. It's very easy to write down a mathematical model that, in theory,
 solved more simply. Let's take a look at a quick system:
 
 ```julia
-@variables u1 u2 u3 u4 u5
+MTK.@variables u1 u2 u3 u4 u5
 eqs = [0 ~ u1 - sin(u5), 0 ~ u2 - cos(u1), 0 ~ u3 - hypot(u1, u2),
     0 ~ u4 - hypot(u2, u3), 0 ~ u5 - hypot(u4, u1)]
-@named sys = NonlinearSystem(eqs, [u1, u2, u3, u4, u5], [])
+MTK.@named sys = MTK.NonlinearSystem(eqs, [u1, u2, u3, u4, u5], [])
 ```
 
 If we run structural simplification, we receive the following form:
 
 ```julia
-sys = structural_simplify(sys)
+sys = MTK.structural_simplify(sys)
 ```
 
 ```julia
-equations(sys)
+MTK.equations(sys)
 ```
 
 How did it do this? Let's look at the `observed` to see the relationships that it found:
 
 ```julia
-observed(sys)
+MTK.observed(sys)
 ```
 
 Using ModelingToolkit, we can build and solve the simplified system:
 
 ```julia
 u0 = [u5 .=> 1.0]
-prob = NonlinearProblem(sys, u0)
-sol = solve(prob, NewtonRaphson())
+prob = NLS.NonlinearProblem(sys, u0)
+sol = NLS.solve(prob, NLS.NewtonRaphson())
 ```
 
 We can then use symbolic indexing to retrieve any variable:
diff --git a/docs/src/tutorials/nonlinear_solve_gpus.md b/docs/src/tutorials/nonlinear_solve_gpus.md
index a2f0871ae..f48566d60 100644
--- a/docs/src/tutorials/nonlinear_solve_gpus.md
+++ b/docs/src/tutorials/nonlinear_solve_gpus.md
@@ -52,13 +52,14 @@ for example
 In practice when this comes together, it looks like:
 
 ```julia
-using NonlinearSolve, CUDA
+import NonlinearSolve as NLS
+import CUDA
 
 f(u, p) = u .* u .- p
-u0 = cu(ones(1000))
-p = cu(collect(1:1000))
-prob = NonlinearProblem(f, u0, p)
-sol = solve(prob, NewtonRaphson(), abstol=1f-4)
+u0 = CUDA.cu(ones(1000))
+p = CUDA.cu(collect(1:1000))
+prob = NLS.NonlinearProblem(f, u0, p)
+sol = NLS.solve(prob, NLS.NewtonRaphson(), abstol=1f-4)
 ```
 
 Notice a few things here. One, nothing is different except the input array types. But
@@ -85,17 +86,19 @@ it. This function must be defined using the `KernelAbstractions.@kernel` macro.
 like:
 
 ```julia
-using NonlinearSolve, StaticArrays
-using KernelAbstractions # For writing vendor-agnostic kernels
-using CUDA # For if you have an NVIDIA GPU
-using AMDGPU # For if you have an AMD GPU
-using Metal # For if you have a Mac M-series device and want to use the built-in GPU
-using OneAPI # For if you have an Intel GPU
-
-@kernel function parallel_nonlinearsolve_kernel!(result, @Const(prob), @Const(alg))
+import NonlinearSolve as NLS
+import StaticArrays
+import SciMLBase
+import KernelAbstractions # For writing vendor-agnostic kernels
+import CUDA # For if you have an NVIDIA GPU
+import AMDGPU # For if you have an AMD GPU
+import Metal # For if you have a Mac M-series device and want to use the built-in GPU
+import OneAPI # For if you have an Intel GPU
+
+@KernelAbstractions.kernel function parallel_nonlinearsolve_kernel!(result, @Const(prob), @Const(alg))
     i = @index(Global)
-    prob_i = remake(prob; p = prob.p[i])
-    sol = solve(prob_i, alg)
+    prob_i = SciMLBase.remake(prob; p = prob.p[i])
+    sol = NLS.solve(prob_i, alg)
     @inbounds result[i] = sol.u
 end
 ```
@@ -137,11 +140,11 @@ Now let's build a nonlinear system to test it on.
     out2 = sqrt(p[2]) * (x[3] - x[4])
     out3 = (x[2] - p[3] * x[3])^2
     out4 = sqrt(p[4]) * (x[1] - x[4]) * (x[1] - x[4])
-    SA[out1,out2,out3,out4]
+    StaticArrays.SA[out1,out2,out3,out4]
 end
 
-p = @SVector [@SVector(rand(Float32, 4)) for _ in 1:1024]
-u0 = SA[1f0, 2f0, 3f0, 4f0]
+p = StaticArrays.@SVector [StaticArrays.@SVector(rand(Float32, 4)) for _ in 1:1024]
+u0 = StaticArrays.SA[1f0, 2f0, 3f0, 4f0]
 prob = SciMLBase.ImmutableNonlinearProblem{false}(p2_f, u0, p)
 ```
 
@@ -161,15 +164,15 @@ and we then simply call our vectorized kernel to parallelize it:
 
 ```julia
 # Threaded CPU
-vectorized_solve(prob, SimpleNewtonRaphson(); backend = CPU())
+vectorized_solve(prob, NLS.SimpleNewtonRaphson(); backend = KernelAbstractions.CPU())
 # AMD ROCM GPU
-vectorized_solve(prob, SimpleNewtonRaphson(); backend = ROCBackend())
+vectorized_solve(prob, NLS.SimpleNewtonRaphson(); backend = AMDGPU.ROCBackend())
 # NVIDIA CUDA GPU
-vectorized_solve(prob, SimpleNewtonRaphson(); backend = CUDABackend())
+vectorized_solve(prob, NLS.SimpleNewtonRaphson(); backend = CUDA.CUDABackend())
 # Intel GPU
-vectorized_solve(prob, SimpleNewtonRaphson(); backend = oneAPI.oneAPIBackend())
+vectorized_solve(prob, NLS.SimpleNewtonRaphson(); backend = OneAPI.oneAPIBackend())
 # Mac M-Series, such as M3Max
-vectorized_solve(prob, SimpleNewtonRaphson(); backend = Metal.MetalBackend())
+vectorized_solve(prob, NLS.SimpleNewtonRaphson(); backend = Metal.MetalBackend())
 ```
 
 !!! warn
diff --git a/docs/src/tutorials/optimizing_parameterized_ode.md b/docs/src/tutorials/optimizing_parameterized_ode.md
index 69cca060d..49d6c5739 100644
--- a/docs/src/tutorials/optimizing_parameterized_ode.md
+++ b/docs/src/tutorials/optimizing_parameterized_ode.md
@@ -4,7 +4,9 @@ Let us fit a parameterized ODE to some data. We will use the Lotka-Volterra mode
 example. We will use Single Shooting to fit the parameters.
 
 ```@example parameterized_ode
-using OrdinaryDiffEqTsit5, NonlinearSolve, Plots
+import OrdinaryDiffEqTsit5 as ODE
+import NonlinearSolve as NLS
+import Plots
 ```
 
 Let us simulate some real data from the Lotka-Volterra model.
@@ -28,32 +30,31 @@ tsteps = 0.0:0.1:10.0
 p = [1.5, 1.0, 3.0, 1.0]
 
 # Setup the ODE problem, then solve
-prob = ODEProblem(lotka_volterra!, u0, tspan, p)
-sol = solve(prob, Tsit5(); saveat = tsteps)
+prob = ODE.ODEProblem(lotka_volterra!, u0, tspan, p)
+sol = ODE.solve(prob, ODE.Tsit5(); saveat = tsteps)
 
 # Plot the solution
-using Plots
-plot(sol; linewidth = 3)
+Plots.plot(sol; linewidth = 3)
 ```
 
 Let us now formulate the parameter estimation as a Nonlinear Least Squares Problem.
 
 ```@example parameterized_ode
 function loss_function(ode_param, data)
-    sol = solve(prob, Tsit5(); p = ode_param, saveat = tsteps)
+    sol = ODE.solve(prob, ODE.Tsit5(); p = ode_param, saveat = tsteps)
     return vec(reduce(hcat, sol.u)) .- data
 end
 
 p_init = zeros(4)
 
-nlls_prob = NonlinearLeastSquaresProblem(loss_function, p_init, vec(reduce(hcat, sol.u)))
+nlls_prob = NLS.NonlinearLeastSquaresProblem(loss_function, p_init, vec(reduce(hcat, sol.u)))
 ```
 
 Now, we can use any NLLS solver to solve this problem.
 
 ```@example parameterized_ode
-res = solve(nlls_prob, LevenbergMarquardt(); maxiters = 1000, show_trace = Val(true),
-    trace_level = TraceWithJacobianConditionNumber(25))
+res = NLS.solve(nlls_prob, NLS.LevenbergMarquardt(); maxiters = 1000, show_trace = Val(true),
+    trace_level = NLS.TraceWithJacobianConditionNumber(25))
 nothing # hide
 ```
 
@@ -64,8 +65,8 @@ res
 We can also use Trust Region methods.
 
 ```@example parameterized_ode
-res = solve(nlls_prob, TrustRegion(); maxiters = 1000, show_trace = Val(true),
-    trace_level = TraceWithJacobianConditionNumber(25))
+res = NLS.solve(nlls_prob, NLS.TrustRegion(); maxiters = 1000, show_trace = Val(true),
+    trace_level = NLS.TraceWithJacobianConditionNumber(25))
 nothing # hide
 ```
 
@@ -76,9 +77,9 @@ res
 Let's plot the solution.
 
 ```@example parameterized_ode
-prob2 = remake(prob; tspan = (0.0, 10.0))
-sol_fit = solve(prob2, Tsit5(); p = res.u)
-sol_true = solve(prob2, Tsit5(); p = p)
-plot(sol_true; linewidth = 3)
-plot!(sol_fit; linewidth = 3, linestyle = :dash)
+prob2 = ODE.remake(prob; tspan = (0.0, 10.0))
+sol_fit = ODE.solve(prob2, ODE.Tsit5(); p = res.u)
+sol_true = ODE.solve(prob2, ODE.Tsit5(); p = p)
+Plots.plot(sol_true; linewidth = 3)
+Plots.plot!(sol_fit; linewidth = 3, linestyle = :dash)
 ```
diff --git a/docs/src/tutorials/small_compile.md b/docs/src/tutorials/small_compile.md
index ec28361f4..70999caab 100644
--- a/docs/src/tutorials/small_compile.md
+++ b/docs/src/tutorials/small_compile.md
@@ -8,13 +8,13 @@ NonlinearSolve.jl can be done with SimpleNonlinearSolve.jl. Thus for example, we
 the core tutorial problem with just SimpleNonlinearSolve.jl as follows:
 
 ```@example simple
-using SimpleNonlinearSolve
+import SimpleNonlinearSolve as SNLS
 
 f(u, p) = u .* u .- p
 u0 = [1.0, 1.0]
 p = 2.0
-prob = NonlinearProblem(f, u0, p)
-sol = solve(prob, SimpleNewtonRaphson())
+prob = SNLS.NonlinearProblem(f, u0, p)
+sol = SNLS.solve(prob, SNLS.SimpleNewtonRaphson())
 ```
 
 However, there are a few downsides to SimpleNonlinearSolve's `SimpleX` style algorithms to
@@ -37,12 +37,12 @@ As such, you can use the code as shown above to have very low startup with good
 for more scaling and debuggability we recommend the full NonlinearSolve.jl. But that said,
 
 ```@example simple
-using StaticArrays
+import StaticArrays
 
-u0 = SA[1.0, 1.0]
+u0 = StaticArrays.SA[1.0, 1.0]
 p = 2.0
-prob = NonlinearProblem(f, u0, p)
-sol = solve(prob, SimpleNewtonRaphson())
+prob = SNLS.NonlinearProblem(f, u0, p)
+sol = SNLS.solve(prob, SNLS.SimpleNewtonRaphson())
 ```
 
 using StaticArrays.jl is also the fastest form for small equations, so if you know your
diff --git a/docs/src/tutorials/snes_ex2.md b/docs/src/tutorials/snes_ex2.md
index 050a39611..62ac81e9d 100644
--- a/docs/src/tutorials/snes_ex2.md
+++ b/docs/src/tutorials/snes_ex2.md
@@ -8,7 +8,11 @@ This solves the equations sequentially. Newton method to solve
 `u'' + u^{2} = f`, sequentially.
 
 ```@example snes_ex2
-using NonlinearSolve, PETSc, LinearAlgebra, SparseConnectivityTracer, BenchmarkTools
+import NonlinearSolve as NLS
+import PETSc
+import LinearAlgebra
+import SparseConnectivityTracer
+import BenchmarkTools
 
 u0 = fill(0.5, 128)
 
@@ -33,25 +37,25 @@ To use automatic sparsity detection, we need to specify `sparsity` keyword argum
 details.
 
 ```@example snes_ex2
-nlfunc_dense = NonlinearFunction(form_residual!)
-nlfunc_sparse = NonlinearFunction(form_residual!; sparsity = TracerSparsityDetector())
+nlfunc_dense = NLS.NonlinearFunction(form_residual!)
+nlfunc_sparse = NLS.NonlinearFunction(form_residual!; sparsity = SparseConnectivityTracer.TracerSparsityDetector())
 
-nlprob_dense = NonlinearProblem(nlfunc_dense, u0)
-nlprob_sparse = NonlinearProblem(nlfunc_sparse, u0)
+nlprob_dense = NLS.NonlinearProblem(nlfunc_dense, u0)
+nlprob_sparse = NLS.NonlinearProblem(nlfunc_sparse, u0)
 ```
 
 Now we can solve the problem using `PETScSNES` or with one of the native `NonlinearSolve.jl`
 solvers.
 
 ```@example snes_ex2
-sol_dense_nr = solve(nlprob_dense, NewtonRaphson(); abstol = 1e-8)
-sol_dense_snes = solve(nlprob_dense, PETScSNES(); abstol = 1e-8)
+sol_dense_nr = NLS.solve(nlprob_dense, NLS.NewtonRaphson(); abstol = 1e-8)
+sol_dense_snes = NLS.solve(nlprob_dense, NLS.PETScSNES(); abstol = 1e-8)
 sol_dense_nr .- sol_dense_snes
 ```
 
 ```@example snes_ex2
-sol_sparse_nr = solve(nlprob_sparse, NewtonRaphson(); abstol = 1e-8)
-sol_sparse_snes = solve(nlprob_sparse, PETScSNES(); abstol = 1e-8)
+sol_sparse_nr = NLS.solve(nlprob_sparse, NLS.NewtonRaphson(); abstol = 1e-8)
+sol_sparse_snes = NLS.solve(nlprob_sparse, NLS.PETScSNES(); abstol = 1e-8)
 sol_sparse_nr .- sol_sparse_snes
 ```
 
@@ -63,19 +67,19 @@ runtimes.
 ### Dense Jacobian
 
 ```@example snes_ex2
-@benchmark solve($(nlprob_dense), $(NewtonRaphson()); abstol = 1e-8)
+BenchmarkTools.@benchmark NLS.solve($(nlprob_dense), $(NLS.NewtonRaphson()); abstol = 1e-8)
 ```
 
 ```@example snes_ex2
-@benchmark solve($(nlprob_dense), $(PETScSNES()); abstol = 1e-8)
+BenchmarkTools.@benchmark NLS.solve($(nlprob_dense), $(NLS.PETScSNES()); abstol = 1e-8)
 ```
 
 ### Sparse Jacobian
 
 ```@example snes_ex2
-@benchmark solve($(nlprob_sparse), $(NewtonRaphson()); abstol = 1e-8)
+BenchmarkTools.@benchmark NLS.solve($(nlprob_sparse), $(NLS.NewtonRaphson()); abstol = 1e-8)
 ```
 
 ```@example snes_ex2
-@benchmark solve($(nlprob_sparse), $(PETScSNES()); abstol = 1e-8)
+BenchmarkTools.@benchmark NLS.solve($(nlprob_sparse), $(NLS.PETScSNES()); abstol = 1e-8)
 ```

From 23f5a671cda8c0dd4308ab83ee378a696cb3d930 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 01:07:00 -0400
Subject: [PATCH 02/15] Update faq.md

---
 docs/src/basics/faq.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/basics/faq.md b/docs/src/basics/faq.md
index 4e37612a4..5c95586b8 100644
--- a/docs/src/basics/faq.md
+++ b/docs/src/basics/faq.md
@@ -58,7 +58,7 @@ function fff_incorrect(var, p)
 end
 
 v_true = [1.0, 0.1, 2.0, 0.5]
-v_init = v_true .+ randn!(similar(v_true)) * 0.1
+v_init = v_true .+ Random.randn!(similar(v_true)) * 0.1
 
 prob_oop = NLS.NonlinearLeastSquaresProblem{false}(fff_incorrect, v_init)
 try
@@ -145,7 +145,7 @@ it will be dynamic and lead to dynamic dispatch. To fix this, we directly specif
 chunksize:
 
 ```@example type_unstable
-@code_warntype solve(
+InteractiveUtils.@code_warntype solve(
     prob,
     NewtonRaphson(;
         autodiff = AutoForwardDiff(; chunksize = NonlinearSolve.pickchunksize(prob.u0))

From 19713c44dd11a5dad2694a2627fd2826b646431e Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 01:07:35 -0400
Subject: [PATCH 03/15] Update code_optimization.md

---
 docs/src/tutorials/code_optimization.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/tutorials/code_optimization.md b/docs/src/tutorials/code_optimization.md
index 1602e3683..b26875ad4 100644
--- a/docs/src/tutorials/code_optimization.md
+++ b/docs/src/tutorials/code_optimization.md
@@ -71,7 +71,7 @@ function f(du, u, p)
 end
 
 prob = NLS.NonlinearProblem(f, u0, p)
-@benchmark sol = NLS.solve(prob, NLS.NewtonRaphson())
+BenchmarkTools.@benchmark sol = NLS.solve(prob, NLS.NewtonRaphson())
 ```
 
 Notice how much faster this already runs! We can make this code even simpler by using
@@ -83,7 +83,7 @@ function f(du, u, p)
     return nothing
 end
 
-@benchmark sol = NLS.solve(prob, NLS.NewtonRaphson())
+BenchmarkTools.@benchmark sol = NLS.solve(prob, NLS.NewtonRaphson())
 ```
 
 ## Further Optimizations for Small Nonlinear Solves with Static Arrays and SimpleNonlinearSolve

From 1be70fb0d200878c6c24cf6f4347aae46a1ea374 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 01:08:30 -0400
Subject: [PATCH 04/15] Update large_systems.md

---
 docs/src/tutorials/large_systems.md | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/docs/src/tutorials/large_systems.md b/docs/src/tutorials/large_systems.md
index 86cf6cbd0..86a5ee67a 100644
--- a/docs/src/tutorials/large_systems.md
+++ b/docs/src/tutorials/large_systems.md
@@ -148,12 +148,12 @@ prob_brusselator_2d_autosparse = NLS.NonlinearProblem(
 )
 
 BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_autosparse,
-    NLS.NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 12)));
+    NLS.NewtonRaphson(; autodiff = ADTypes.AutoForwardDiff(; chunksize = 12)));
 BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_autosparse,
-    NLS.NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 12),
+    NLS.NewtonRaphson(; autodiff = ADTypes.AutoForwardDiff(; chunksize = 12),
         linsolve = LS.KLUFactorization()));
 BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_autosparse,
-    NLS.NewtonRaphson(; autodiff = AutoForwardDiff(; chunksize = 12),
+    NLS.NewtonRaphson(; autodiff = ADTypes.AutoForwardDiff(; chunksize = 12),
         linsolve = LS.KrylovJL_GMRES()));
 nothing # hide
 ```
@@ -325,7 +325,7 @@ prob_brusselator_2d_exact_tracer = NLS.NonlinearProblem(
     u0, p; abstol = 1e-10, reltol = 1e-10)
 prob_brusselator_2d_approx_di = NLS.NonlinearProblem(
     NLS.NonlinearFunction(brusselator_2d_loop;
-        sparsity = DifferentiationInterface.DenseSparsityDetector(DifferentiationInterface.AutoForwardDiff(); atol = 1e-4)),
+        sparsity = DifferentiationInterface.DenseSparsityDetector(ADTypes.AutoForwardDiff(); atol = 1e-4)),
     u0, p; abstol = 1e-10, reltol = 1e-10)
 
 BenchmarkTools.@btime NLS.solve(prob_brusselator_2d_exact_tracer, NLS.NewtonRaphson());

From 9794b529553358ed9708b778130d28437f4cd509 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 01:09:54 -0400
Subject: [PATCH 05/15] Update Project.toml

---
 docs/Project.toml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/Project.toml b/docs/Project.toml
index a98646942..ce2fa3792 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -17,6 +17,7 @@ NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
 NonlinearSolveHomotopyContinuation = "2ac3b008-d579-4536-8c91-a1a5998c2f8b"
 NonlinearSolveQuasiNewton = "9a2c21bd-3a47-402d-9113-8faf9a0ee114"
+NonlinearSolveSciPy = "4827a3aa-8a82-4c61-8bd0-3c7d3e464ee5"
 NonlinearSolveSpectralMethods = "26075421-4e9a-44e1-8bd1-420ed7ad02b2"
 OrdinaryDiffEqTsit5 = "b1df2697-797e-41e3-8120-5422d3b24e4a"
 PETSc = "ace2c81b-2b5f-4b1e-a30d-d662738edfe0"
@@ -49,6 +50,7 @@ NonlinearSolveBase = "1"
 NonlinearSolveFirstOrder = "1"
 NonlinearSolveHomotopyContinuation = "0.1"
 NonlinearSolveQuasiNewton = "1"
+NonlinearSolveSciPy = "1"
 NonlinearSolveSpectralMethods = "1"
 OrdinaryDiffEqTsit5 = "1.1.0"
 PETSc = "0.3"

From 6e65c17c5653a22d878abb5aff1668056c13cfd6 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 01:10:19 -0400
Subject: [PATCH 06/15] Update make.jl

---
 docs/make.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/make.jl b/docs/make.jl
index d2915275a..746c26657 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -5,7 +5,7 @@ using Sundials
 using NonlinearSolveBase, SciMLBase, DiffEqBase
 using SimpleNonlinearSolve, BracketingNonlinearSolve
 using NonlinearSolveFirstOrder, NonlinearSolveQuasiNewton, NonlinearSolveSpectralMethods
-using NonlinearSolveHomotopyContinuation
+using NonlinearSolveHomotopyContinuation, NonlinearSolveSciPy
 using SciMLJacobianOperators
 using NonlinearSolve, SteadyStateDiffEq
 
@@ -38,7 +38,7 @@ makedocs(;
         NonlinearSolveFirstOrder, NonlinearSolveQuasiNewton, NonlinearSolveSpectralMethods,
         NonlinearSolveHomotopyContinuation,
         Sundials,
-        SciMLJacobianOperators,
+        SciMLJacobianOperators, NonlinearSolveSciPy,
         NonlinearSolve, SteadyStateDiffEq
     ],
     clean = true,

From db33ec9cc56433dcc41fc732d25fa8522eb19f0d Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 01:13:25 -0400
Subject: [PATCH 07/15] Update Project.toml

---
 docs/Project.toml | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/docs/Project.toml b/docs/Project.toml
index ce2fa3792..8c2f09069 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -12,6 +12,7 @@ DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
 IncompleteLU = "40713840-3770-5561-ab4c-a76e7d0d7895"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
+LineSearch = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
 NonlinearSolveBase = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
 NonlinearSolveFirstOrder = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
@@ -45,6 +46,7 @@ DocumenterInterLinks = "1.0.0"
 IncompleteLU = "0.2"
 InteractiveUtils = "<0.0.1, 1"
 LinearSolve = "2, 3"
+LineSearch = "0.1"
 NonlinearSolve = "4"
 NonlinearSolveBase = "1"
 NonlinearSolveFirstOrder = "1"

From 1e0c7d5b5d0f78d1256218d6f670596a3a9bc04e Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 01:13:39 -0400
Subject: [PATCH 08/15] Update make.jl

---
 docs/make.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/make.jl b/docs/make.jl
index 746c26657..ce2876953 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -6,7 +6,7 @@ using NonlinearSolveBase, SciMLBase, DiffEqBase
 using SimpleNonlinearSolve, BracketingNonlinearSolve
 using NonlinearSolveFirstOrder, NonlinearSolveQuasiNewton, NonlinearSolveSpectralMethods
 using NonlinearSolveHomotopyContinuation, NonlinearSolveSciPy
-using SciMLJacobianOperators
+using SciMLJacobianOperators, LineSearch
 using NonlinearSolve, SteadyStateDiffEq
 
 cp(
@@ -37,7 +37,7 @@ makedocs(;
         SimpleNonlinearSolve, BracketingNonlinearSolve,
         NonlinearSolveFirstOrder, NonlinearSolveQuasiNewton, NonlinearSolveSpectralMethods,
         NonlinearSolveHomotopyContinuation,
-        Sundials,
+        Sundials, LineSearch,
         SciMLJacobianOperators, NonlinearSolveSciPy,
         NonlinearSolve, SteadyStateDiffEq
     ],

From 246f0b1e574fb3dca83ec77719fd5922f85e7a86 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 08:50:23 -0400
Subject: [PATCH 09/15] Update faq.md

---
 docs/src/basics/faq.md | 12 +++++++-----
 1 file changed, 7 insertions(+), 5 deletions(-)

diff --git a/docs/src/basics/faq.md b/docs/src/basics/faq.md
index 5c95586b8..aa76c7a34 100644
--- a/docs/src/basics/faq.md
+++ b/docs/src/basics/faq.md
@@ -76,7 +76,8 @@ be a Dual number. This causes the error. To fix it:
  1. Specify the `autodiff` to be `AutoFiniteDiff`
     
     ```@example dual_error_faq
-    sol = solve(prob_oop, LevenbergMarquardt(; autodiff = AutoFiniteDiff());
+    import ADTypes
+    sol = solve(prob_oop, LevenbergMarquardt(; autodiff = ADTypes.AutoFiniteDiff());
         maxiters = 10000, abstol = 1e-8)
     ```
     
@@ -110,6 +111,7 @@ internally. See this simple example:
 ```@example type_unstable
 import NonlinearSolve as NLS
 import InteractiveUtils
+import ADTypes
 
 f(u, p) = @. u^2 - p
 
@@ -134,7 +136,7 @@ nothing # hide
 Again Type-Stable! Now let's try using a regular array:
 
 ```@example type_unstable
-prob = NonlinearProblem(f, [1.0, 2.0], 2.0)
+prob = NLS.NonlinearProblem(f, [1.0, 2.0], 2.0)
 
 InteractiveUtils.@code_warntype NLS.solve(prob, NLS.NewtonRaphson())
 nothing # hide
@@ -145,10 +147,10 @@ it will be dynamic and lead to dynamic dispatch. To fix this, we directly specif
 chunksize:
 
 ```@example type_unstable
-InteractiveUtils.@code_warntype solve(
+InteractiveUtils.@code_warntype NLS.solve(
     prob,
-    NewtonRaphson(;
-        autodiff = AutoForwardDiff(; chunksize = NonlinearSolve.pickchunksize(prob.u0))
+    NLS.NewtonRaphson(;
+        autodiff = ADTypes.AutoForwardDiff(; chunksize = NonlinearSolve.pickchunksize(prob.u0))
     )
 )
 nothing # hide

From b7e8bbefff07b4984a6255c785ac8219a9d2864f Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 08:50:56 -0400
Subject: [PATCH 10/15] Update large_systems.md

---
 docs/src/tutorials/large_systems.md | 1 +
 1 file changed, 1 insertion(+)

diff --git a/docs/src/tutorials/large_systems.md b/docs/src/tutorials/large_systems.md
index 86a5ee67a..294cc32b0 100644
--- a/docs/src/tutorials/large_systems.md
+++ b/docs/src/tutorials/large_systems.md
@@ -64,6 +64,7 @@ import NonlinearSolve as NLS
 import LinearAlgebra
 import SparseArrays
 import LinearSolve as LS
+import ADTypes
 
 const N = 32
 const xyd_brusselator = range(0, stop = 1, length = N)

From af7064ffd5c2c2cae2312a5b08599437e6a1cec7 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 10:44:43 -0400
Subject: [PATCH 11/15] Update faq.md

---
 docs/src/basics/faq.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/basics/faq.md b/docs/src/basics/faq.md
index aa76c7a34..979cb4f08 100644
--- a/docs/src/basics/faq.md
+++ b/docs/src/basics/faq.md
@@ -150,7 +150,7 @@ chunksize:
 InteractiveUtils.@code_warntype NLS.solve(
     prob,
     NLS.NewtonRaphson(;
-        autodiff = ADTypes.AutoForwardDiff(; chunksize = NonlinearSolve.pickchunksize(prob.u0))
+        autodiff = ADTypes.AutoForwardDiff(; chunksize = NLS.pickchunksize(prob.u0))
     )
 )
 nothing # hide

From c2b697a009d08d2a29c69642b049038d9ca1947d Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 10:45:16 -0400
Subject: [PATCH 12/15] Update faq.md

---
 docs/src/basics/faq.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/basics/faq.md b/docs/src/basics/faq.md
index 979cb4f08..f54283023 100644
--- a/docs/src/basics/faq.md
+++ b/docs/src/basics/faq.md
@@ -77,7 +77,7 @@ be a Dual number. This causes the error. To fix it:
     
     ```@example dual_error_faq
     import ADTypes
-    sol = solve(prob_oop, LevenbergMarquardt(; autodiff = ADTypes.AutoFiniteDiff());
+    sol = solve(prob_oop, NLS.LevenbergMarquardt(; autodiff = ADTypes.AutoFiniteDiff());
         maxiters = 10000, abstol = 1e-8)
     ```
     

From c0cb859000cef8eab08bc5c246a76cc2186e22e9 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 10:53:48 -0400
Subject: [PATCH 13/15] Add SciPy page

---
 docs/pages.jl         |  1 +
 docs/src/api/scipy.md | 25 +++++++++++++++++++++++++
 2 files changed, 26 insertions(+)
 create mode 100644 docs/src/api/scipy.md

diff --git a/docs/pages.jl b/docs/pages.jl
index 5a95aed37..23e449bba 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -48,6 +48,7 @@ pages = [
         "api/nlsolve.md",
         "api/nlsolvers.md",
         "api/petsc.md",
+        "api/scipy.md",
         "api/siamfanlequations.md",
         "api/speedmapping.md",
         "api/sundials.md",
diff --git a/docs/src/api/scipy.md b/docs/src/api/scipy.md
new file mode 100644
index 000000000..bae9ba417
--- /dev/null
+++ b/docs/src/api/scipy.md
@@ -0,0 +1,25 @@
+# SciPy
+
+This is an extension for importing solvers from
+[SciPy](https://scipy.org/) 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:
+
+```julia
+import Pkg
+Pkg.add("NonlinearSolveSciPy")
+import NonlinearSolveSciPy as NLSP
+```
+
+Note that using this package requires Python and SciPy to be available via PythonCall.jl.
+
+## Solver API
+
+```@docs
+NonlinearSolveSciPy.SciPyLeastSquares
+NonlinearSolveSciPy.SciPyLeastSquaresTRF
+NonlinearSolveSciPy.SciPyLeastSquaresDogbox
+NonlinearSolveSciPy.SciPyLeastSquaresLM
+NonlinearSolveSciPy.SciPyRoot
+NonlinearSolveSciPy.SciPyRootScalar
+```
\ No newline at end of file

From 6979276f612cde5f2ad2eb4dce7dcc482e432429 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 13:16:47 -0400
Subject: [PATCH 14/15] Update faq.md

---
 docs/src/basics/faq.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/basics/faq.md b/docs/src/basics/faq.md
index f54283023..da2cbf8ba 100644
--- a/docs/src/basics/faq.md
+++ b/docs/src/basics/faq.md
@@ -77,7 +77,7 @@ be a Dual number. This causes the error. To fix it:
     
     ```@example dual_error_faq
     import ADTypes
-    sol = solve(prob_oop, NLS.LevenbergMarquardt(; autodiff = ADTypes.AutoFiniteDiff());
+    sol = NLS.solve(prob_oop, NLS.LevenbergMarquardt(; autodiff = ADTypes.AutoFiniteDiff());
         maxiters = 10000, abstol = 1e-8)
     ```
     

From 529aea8699f31881fa29e33a027b1b65a793e100 Mon Sep 17 00:00:00 2001
From: Christopher Rackauckas <accounts@chrisrackauckas.com>
Date: Sun, 20 Jul 2025 13:18:08 -0400
Subject: [PATCH 15/15] Update docs/src/api/scipy.md

---
 docs/src/api/scipy.md | 3 ---
 1 file changed, 3 deletions(-)

diff --git a/docs/src/api/scipy.md b/docs/src/api/scipy.md
index bae9ba417..271b5591c 100644
--- a/docs/src/api/scipy.md
+++ b/docs/src/api/scipy.md
@@ -17,9 +17,6 @@ Note that using this package requires Python and SciPy to be available via Pytho
 
 ```@docs
 NonlinearSolveSciPy.SciPyLeastSquares
-NonlinearSolveSciPy.SciPyLeastSquaresTRF
-NonlinearSolveSciPy.SciPyLeastSquaresDogbox
-NonlinearSolveSciPy.SciPyLeastSquaresLM
 NonlinearSolveSciPy.SciPyRoot
 NonlinearSolveSciPy.SciPyRootScalar
 ```
\ No newline at end of file