From 5de311981aa24b1b43f58af261692b71c37f154e Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 17 May 2025 22:44:30 +0800
Subject: [PATCH 1/5] Add maxsol and minsol

---
 docs/pages.jl                                 |  3 +-
 docs/src/tutorials/inequality.md              | 28 +++++++
 .../src/solution_utils.jl                     |  3 +
 .../src/BoundaryValueDiffEqMIRK.jl            |  1 +
 .../src/interpolation.jl                      | 83 +++++++++++++++++++
 .../test/mirk_basic_tests.jl                  | 17 ++++
 6 files changed, 134 insertions(+), 1 deletion(-)
 create mode 100644 docs/src/tutorials/inequality.md

diff --git a/docs/pages.jl b/docs/pages.jl
index 0acc961a..b9e078ae 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -2,7 +2,8 @@
 
 pages = ["index.md",
     "Getting Started with BVP solving in Julia" => "tutorials/getting_started.md",
-    "Tutorials" => Any["tutorials/continuation.md", "tutorials/solve_nlls_bvp.md"],
+    "Tutorials" => Any["tutorials/continuation.md",
+        "tutorials/solve_nlls_bvp.md", "tutorials/inequality.md"],
     "Basics" => Any["basics/bvp_problem.md", "basics/bvp_functions.md",
         "basics/solve.md", "basics/autodiff.md", "basics/error_control.md"],
     "Solver Summaries and Recommendations" => Any[
diff --git a/docs/src/tutorials/inequality.md b/docs/src/tutorials/inequality.md
new file mode 100644
index 00000000..45314c8e
--- /dev/null
+++ b/docs/src/tutorials/inequality.md
@@ -0,0 +1,28 @@
+# Solve BVP with Inequality Boundary Conditions
+
+When dealing with scenarios that boundary conditions are stronger than just explicit values at specific points—such as inequality boundary conditions, where the solution must satisfy constraints like staying within upper and lower bounds—we need a more flexible approach. In such cases, we can use the `maxsol` and `minsol` functions provided by BoundaryValueDiffEq.jl to specify such inequality conditions.
+
+Let's walk through this functionlity with an intuitive example. We still revisit the simple pendulum example here, but this time, we’ll impose upper and lower bound constraints on the solution, specified as:
+
+```math
+lb \leq u \leq ub
+```
+
+where `lb=-4.8161991710010925` and `ub=5.0496477654230745`. So the states must be bigger than `lb` but smaller than `ub`. To solve such problems, we can simply use the `minsol` and `maxsol` functions when defining the boundary value problem in BoundaryValueDiffEq.jl.
+
+```@example inequality
+tspan = (0.0, pi / 2)
+function simplependulum!(du, u, p, t)
+    θ = u[1]
+    dθ = u[2]
+    du[1] = dθ
+    du[2] = -9.81 * sin(θ)
+end
+function bc!(residual, u, p, t)
+    residual[1] = maxsol(u, (0.0, pi / 2)) - 5.0496477654230745
+    residual[2] = minsol(u, (0.0, pi / 2)) + 4.8161991710010925
+end
+prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
+sol = solve(prob, MIRK4(), dt = 0.05)
+plot(sol)
+```
diff --git a/lib/BoundaryValueDiffEqCore/src/solution_utils.jl b/lib/BoundaryValueDiffEqCore/src/solution_utils.jl
index ed0002a9..ae973a04 100644
--- a/lib/BoundaryValueDiffEqCore/src/solution_utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/solution_utils.jl
@@ -29,3 +29,6 @@ Base.eltype(e::EvalSol) = eltype(e.u)
 function Base.show(io::IO, m::MIME"text/plain", e::EvalSol)
     (print(io, "t: "); show(io, m, e.t); println(io); print(io, "u: "); show(io, m, e.u))
 end
+
+Base.maximum(sol::EvalSol) = maximum(Iterators.flatten(sol.u))
+Base.minimum(sol::EvalSol) = minimum(Iterators.flatten(sol.u))
diff --git a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
index 1291c8d0..629ab929 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
@@ -160,5 +160,6 @@ end
 
 export MIRK2, MIRK3, MIRK4, MIRK5, MIRK6, MIRK6I
 export BVPJacobianAlgorithm
+export maxsol, minsol
 
 end
diff --git a/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl b/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
index 4f3b1c37..8f1341a1 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
@@ -158,6 +158,89 @@ function (s::EvalSol{C})(tval::Number) where {C <: MIRKCache}
     return z
 end
 
+# Interpolate intermidiate solution at multiple points
+function (s::EvalSol{C})(tvals::AbstractArray{<:Number}) where {C <: MIRKCache}
+    (; t, u, cache) = s
+    (; alg, stage, k_discrete, mesh_dt) = cache
+    # Quick handle for the case where tval is at the boundary
+    zvals = [zero(last(u)) for _ in tvals]
+    for (i, tval) in enumerate(tvals)
+        (tval == t[1]) && return first(u)
+        (tval == t[end]) && return last(u)
+        ii = interval(t, tval)
+        dt = mesh_dt[ii]
+        τ = (tval - t[ii]) / dt
+        w, _ = evalsol_interp_weights(τ, alg)
+        K = __needs_diffcache(alg.jac_alg) ? @view(k_discrete[ii].du[:, 1:stage]) :
+            @view(k_discrete[ii][:, 1:stage])
+        __maybe_matmul!(zvals[i], K, @view(w[1:stage]))
+        zvals[i] .= zvals[i] .* dt .+ u[ii]
+    end
+    return zvals
+end
+
+# Intermidiate derivative solution for evaluating boundry conditions
+function (s::EvalSol{C})(tval::Number, ::Type{Val{1}}) where {C <: MIRKCache}
+    (; t, u, cache) = s
+    (; alg, stage, k_discrete, mesh_dt) = cache
+    z′ = zeros(typeof(tval), 2)
+    ii = interval(t, tval)
+    dt = mesh_dt[ii]
+    τ = (tval - t[ii]) / dt
+    _, w′ = interp_weights(τ, alg)
+    __maybe_matmul!(z′, @view(k_discrete[ii].du[:, 1:stage]), @view(w′[1:stage]))
+    return z′
+end
+
+"""
+n root-finding problems to find the critical points with continuous derivative polynomials
+"""
+function __construct_then_solve_root_problem(
+        sol::EvalSol{C}, tspan::Tuple) where {C <: MIRKCache}
+    (; alg) = sol.cache
+    n = first(size(sol))
+    nlprobs = Vector{NonlinearProblem}(undef, n)
+    for i in 1:n
+        f = @closure (t, p) -> sol(t, Val{1})[i]
+        nlprob = NonlinearProblem(f, sol.cache.prob.u0[i], tspan)
+        nlprobs[i] = nlprob
+    end
+    nlsols = Vector{SciMLBase.NonlinearSolution}(undef, length(nlprobs))
+    nlsolve_alg = __concrete_nonlinearsolve_algorithm(first(nlprobs), alg.nlsolve)
+    for (i, nlprob) in enumerate(nlprobs)
+        nlsols[i] = solve(nlprob, nlsolve_alg)
+    end
+    return nlsols
+end
+
+# It turns out the critical points can't cover all possible maximum/minimum values
+# especially when the solution are monotonic, we still need to compare the extremes with
+# value at critical points to find the maximum/minimum
+
+"""
+    maxsol(sol::EvalSol, tspan::Tuple)
+
+Find the maximum of the solution over the time span `tspan`.
+"""
+function maxsol(sol::EvalSol{C}, tspan::Tuple) where {C <: MIRKCache}
+    nlsols = __construct_then_solve_root_problem(sol, tspan)
+    tvals = map(nlsol -> (SciMLBase.successful_retcode(nlsol); return nlsol.u), nlsols)
+    u = sol(tvals)
+    return max(maximum(sol), maximum(Iterators.flatten(u)))
+end
+
+"""
+    minsol(sol::EvalSol, tspan::Tuple)
+
+Find the minimum of the solution over the time span `tspan`.
+"""
+function minsol(sol::EvalSol{C}, tspan::Tuple) where {C <: MIRKCache}
+    nlsols = __construct_then_solve_root_problem(sol, tspan)
+    tvals = map(nlsol -> (SciMLBase.successful_retcode(nlsol); return nlsol.u), nlsols)
+    u = sol(tvals)
+    return min(minimum(sol), minimum(Iterators.flatten(u)))
+end
+
 @inline function evalsol_interp_weights(τ::T, ::MIRK2) where {T}
     w = [0, τ * (1 - τ / 2), τ^2 / 2]
 
diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
index ab72a287..a3268ba1 100644
--- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
@@ -401,3 +401,20 @@ end
         @test SciMLBase.successful_retcode(sol)
     end
 end
+
+@testset "Test maxsol and minsol" begin
+    tspan = (0.0, pi / 2)
+    function simplependulum!(du, u, p, t)
+        θ = u[1]
+        dθ = u[2]
+        du[1] = dθ
+        du[2] = -9.81 * sin(θ)
+    end
+    function bc!(residual, u, p, t)
+        residual[1] = maxsol(u, (0.0, pi / 2)) - 5.0496477654230745
+        residual[2] = minsol(u, (0.0, pi / 2)) + 4.8161991710010925
+    end
+    prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
+    sol = solve(prob, MIRK4(), dt = 0.05)
+    @test SciMLBase.successful_retcode(sol)
+end

From 697b367c516eeca15769e8c57a0bc42b427d4349 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sun, 18 May 2025 00:36:46 +0800
Subject: [PATCH 2/5] Fix CI failings

---
 docs/src/tutorials/inequality.md                     | 1 +
 lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl | 8 +++++---
 2 files changed, 6 insertions(+), 3 deletions(-)

diff --git a/docs/src/tutorials/inequality.md b/docs/src/tutorials/inequality.md
index 45314c8e..5a07f647 100644
--- a/docs/src/tutorials/inequality.md
+++ b/docs/src/tutorials/inequality.md
@@ -11,6 +11,7 @@ lb \leq u \leq ub
 where `lb=-4.8161991710010925` and `ub=5.0496477654230745`. So the states must be bigger than `lb` but smaller than `ub`. To solve such problems, we can simply use the `minsol` and `maxsol` functions when defining the boundary value problem in BoundaryValueDiffEq.jl.
 
 ```@example inequality
+using BoundaryValueDiffEq
 tspan = (0.0, pi / 2)
 function simplependulum!(du, u, p, t)
     θ = u[1]
diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
index a3268ba1..349f61da 100644
--- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
@@ -402,7 +402,7 @@ end
     end
 end
 
-@testset "Test maxsol and minsol" begin
+@testitem "Test maxsol and minsol" setup=[MIRKConvergenceTests] begin
     tspan = (0.0, pi / 2)
     function simplependulum!(du, u, p, t)
         θ = u[1]
@@ -415,6 +415,8 @@ end
         residual[2] = minsol(u, (0.0, pi / 2)) + 4.8161991710010925
     end
     prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
-    sol = solve(prob, MIRK4(), dt = 0.05)
-    @test SciMLBase.successful_retcode(sol)
+    for order in (4, 6)
+        sol = solve(prob, mirk_solver(Val(order)), dt = 0.05)
+        @test SciMLBase.successful_retcode(sol)
+    end
 end

From 1a0fe09a0bc1dde4ece35a01cfe91715ef30124c Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sun, 18 May 2025 01:16:16 +0800
Subject: [PATCH 3/5] Parent package export maxsol and minsol

---
 src/BoundaryValueDiffEq.jl | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/src/BoundaryValueDiffEq.jl b/src/BoundaryValueDiffEq.jl
index 9775268b..28f47b64 100644
--- a/src/BoundaryValueDiffEq.jl
+++ b/src/BoundaryValueDiffEq.jl
@@ -32,4 +32,6 @@ export BVPM2, BVPSOL, COLNEW # From ODEInterface.jl
 
 export BVPJacobianAlgorithm
 
+export maxsol, minsol
+
 end

From f885af191bd1edb10f940af70695493b0346741f Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sun, 18 May 2025 02:02:33 +0800
Subject: [PATCH 4/5] Dont forget Plots

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

diff --git a/docs/src/tutorials/inequality.md b/docs/src/tutorials/inequality.md
index 5a07f647..ea8dbcb4 100644
--- a/docs/src/tutorials/inequality.md
+++ b/docs/src/tutorials/inequality.md
@@ -11,7 +11,7 @@ lb \leq u \leq ub
 where `lb=-4.8161991710010925` and `ub=5.0496477654230745`. So the states must be bigger than `lb` but smaller than `ub`. To solve such problems, we can simply use the `minsol` and `maxsol` functions when defining the boundary value problem in BoundaryValueDiffEq.jl.
 
 ```@example inequality
-using BoundaryValueDiffEq
+using BoundaryValueDiffEq, Plots
 tspan = (0.0, pi / 2)
 function simplependulum!(du, u, p, t)
     θ = u[1]

From 055e4f60336775fb273dd3839db1aa3fb7894133 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Wed, 21 May 2025 17:26:43 +0800
Subject: [PATCH 5/5] Actually extremum boundary conditions

---
 docs/pages.jl                                 |  4 +--
 docs/src/tutorials/extremum.md                | 30 +++++++++++++++++++
 docs/src/tutorials/inequality.md              | 29 ------------------
 .../src/default_nlsolve.jl                    | 12 ++++++++
 .../src/BoundaryValueDiffEqMIRK.jl            |  3 +-
 .../src/interpolation.jl                      | 11 +++----
 6 files changed, 50 insertions(+), 39 deletions(-)
 create mode 100644 docs/src/tutorials/extremum.md
 delete mode 100644 docs/src/tutorials/inequality.md

diff --git a/docs/pages.jl b/docs/pages.jl
index b9e078ae..dd11e26e 100644
--- a/docs/pages.jl
+++ b/docs/pages.jl
@@ -2,8 +2,8 @@
 
 pages = ["index.md",
     "Getting Started with BVP solving in Julia" => "tutorials/getting_started.md",
-    "Tutorials" => Any["tutorials/continuation.md",
-        "tutorials/solve_nlls_bvp.md", "tutorials/inequality.md"],
+    "Tutorials" => Any[
+        "tutorials/continuation.md", "tutorials/solve_nlls_bvp.md", "tutorials/extremum.md"],
     "Basics" => Any["basics/bvp_problem.md", "basics/bvp_functions.md",
         "basics/solve.md", "basics/autodiff.md", "basics/error_control.md"],
     "Solver Summaries and Recommendations" => Any[
diff --git a/docs/src/tutorials/extremum.md b/docs/src/tutorials/extremum.md
new file mode 100644
index 00000000..8cd502c5
--- /dev/null
+++ b/docs/src/tutorials/extremum.md
@@ -0,0 +1,30 @@
+# Solve BVP with Extremum Boundary Conditions
+
+In many physical systems, boundary conditions are not always defined at fixed points such as intitial or terminal ends of the domain. Instead, we may encounter scenarios where constraints are imposed on the maximum or minimum values that the solution must attain somewhere within the domain. In such cases, we can use the `maxsol` and `minsol` functions provided by BoundaryValueDiffEq.jl to specify such extremum boundary conditions.
+
+Let's walk through this functionality with an intuitive example. We still revisit the simple pendulum example here, but this time, suppose we need to impose the maximum and minimum value to our boundary conditions, specified as:
+
+```math
+\max{u}=ub\\
+\min{u}=lb
+```
+
+where `lb=-4.8161991710010925` and `ub=5.0496477654230745`. So the states must conform that the maximum value of the state should be `lb` while the minimum value of the state should be `ub`. To solve such problems, we can simply use the `maxsol` and `minsol` functions when defining the boundary value problem in BoundaryValueDiffEq.jl.
+
+```@example inequality
+using BoundaryValueDiffEq, Plots
+tspan = (0.0, pi / 2)
+function simplependulum!(du, u, p, t)
+    θ = u[1]
+    dθ = u[2]
+    du[1] = dθ
+    du[2] = -9.81 * sin(θ)
+end
+function bc!(residual, u, p, t)
+    residual[1] = maxsol(u, (0.0, pi / 2)) - 5.0496477654230745
+    residual[2] = minsol(u, (0.0, pi / 2)) + 4.8161991710010925
+end
+prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
+sol = solve(prob, MIRK4(), dt = 0.05)
+plot(sol)
+```
diff --git a/docs/src/tutorials/inequality.md b/docs/src/tutorials/inequality.md
deleted file mode 100644
index ea8dbcb4..00000000
--- a/docs/src/tutorials/inequality.md
+++ /dev/null
@@ -1,29 +0,0 @@
-# Solve BVP with Inequality Boundary Conditions
-
-When dealing with scenarios that boundary conditions are stronger than just explicit values at specific points—such as inequality boundary conditions, where the solution must satisfy constraints like staying within upper and lower bounds—we need a more flexible approach. In such cases, we can use the `maxsol` and `minsol` functions provided by BoundaryValueDiffEq.jl to specify such inequality conditions.
-
-Let's walk through this functionlity with an intuitive example. We still revisit the simple pendulum example here, but this time, we’ll impose upper and lower bound constraints on the solution, specified as:
-
-```math
-lb \leq u \leq ub
-```
-
-where `lb=-4.8161991710010925` and `ub=5.0496477654230745`. So the states must be bigger than `lb` but smaller than `ub`. To solve such problems, we can simply use the `minsol` and `maxsol` functions when defining the boundary value problem in BoundaryValueDiffEq.jl.
-
-```@example inequality
-using BoundaryValueDiffEq, Plots
-tspan = (0.0, pi / 2)
-function simplependulum!(du, u, p, t)
-    θ = u[1]
-    dθ = u[2]
-    du[1] = dθ
-    du[2] = -9.81 * sin(θ)
-end
-function bc!(residual, u, p, t)
-    residual[1] = maxsol(u, (0.0, pi / 2)) - 5.0496477654230745
-    residual[2] = minsol(u, (0.0, pi / 2)) + 4.8161991710010925
-end
-prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
-sol = solve(prob, MIRK4(), dt = 0.05)
-plot(sol)
-```
diff --git a/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl b/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl
index f62d61d9..cb1ebd42 100644
--- a/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl
+++ b/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl
@@ -36,6 +36,18 @@ function __FastShortcutBVPCompatibleNonlinearPolyalg(
     return NonlinearSolvePolyAlgorithm(algs)
 end
 
+function __FastShortcutNonlinearPolyalg(::Type{T} = Float64; concrete_jac = nothing,
+        linsolve = nothing, autodiff = nothing) where {T}
+    if T <: Complex
+        algs = (NewtonRaphson(; concrete_jac, linsolve, autodiff),)
+    else
+        algs = (NewtonRaphson(; concrete_jac, linsolve, autodiff),
+            NewtonRaphson(; concrete_jac, linsolve, linesearch = BackTracking(), autodiff),
+            TrustRegion(; concrete_jac, linsolve, autodiff))
+    end
+    return NonlinearSolvePolyAlgorithm(algs)
+end
+
 @inline __concrete_nonlinearsolve_algorithm(prob, alg) = alg
 @inline function __concrete_nonlinearsolve_algorithm(prob, ::Nothing)
     if prob isa NonlinearLeastSquaresProblem
diff --git a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
index 629ab929..4297e0e6 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
@@ -21,7 +21,8 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                DefectControl, GlobalErrorControl, SequentialErrorControl,
                                HybridErrorControl, HOErrorControl, __use_both_error_control,
                                __default_coloring_algorithm, DiffCacheNeeded,
-                               NoDiffCacheNeeded, __split_kwargs
+                               NoDiffCacheNeeded, __split_kwargs,
+                               __FastShortcutNonlinearPolyalg
 
 using ConcreteStructs: @concrete
 using DiffEqBase: DiffEqBase
diff --git a/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl b/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
index 8f1341a1..b17db952 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/interpolation.jl
@@ -193,21 +193,18 @@ function (s::EvalSol{C})(tval::Number, ::Type{Val{1}}) where {C <: MIRKCache}
 end
 
 """
-n root-finding problems to find the critical points with continuous derivative polynomials
+Construct n root-finding problems and solve them to find the critical points with continuous derivative polynomials
 """
 function __construct_then_solve_root_problem(
         sol::EvalSol{C}, tspan::Tuple) where {C <: MIRKCache}
     (; alg) = sol.cache
     n = first(size(sol))
-    nlprobs = Vector{NonlinearProblem}(undef, n)
+    nlprobs = Vector{SciMLBase.NonlinearProblem}(undef, n)
+    nlsols = Vector{SciMLBase.NonlinearSolution}(undef, length(nlprobs))
+    nlsolve_alg = __FastShortcutNonlinearPolyalg(eltype(sol.cache))
     for i in 1:n
         f = @closure (t, p) -> sol(t, Val{1})[i]
         nlprob = NonlinearProblem(f, sol.cache.prob.u0[i], tspan)
-        nlprobs[i] = nlprob
-    end
-    nlsols = Vector{SciMLBase.NonlinearSolution}(undef, length(nlprobs))
-    nlsolve_alg = __concrete_nonlinearsolve_algorithm(first(nlprobs), alg.nlsolve)
-    for (i, nlprob) in enumerate(nlprobs)
         nlsols[i] = solve(nlprob, nlsolve_alg)
     end
     return nlsols