From 53dc9f270baa275369f7ba3ef81a1b99a9ab313e Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Tue, 15 Jul 2025 00:43:05 +0800
Subject: [PATCH 01/50] Optimization based solvers

---
 .../src/BoundaryValueDiffEqCore.jl            |   2 +-
 lib/BoundaryValueDiffEqCore/src/algorithms.jl |   9 ++
 ...t_nlsolve.jl => default_internal_solve.jl} |  16 ++-
 lib/BoundaryValueDiffEqCore/src/utils.jl      |  61 ++++++++-
 lib/BoundaryValueDiffEqMIRK/Project.toml      |   4 +-
 .../src/BoundaryValueDiffEqMIRK.jl            | 125 ++----------------
 lib/BoundaryValueDiffEqMIRK/src/algorithms.jl |  12 +-
 lib/BoundaryValueDiffEqMIRK/src/mirk.jl       |  41 +++---
 .../test/mirk_basic_tests.jl                  |  20 +++
 9 files changed, 148 insertions(+), 142 deletions(-)
 rename lib/BoundaryValueDiffEqCore/src/{default_nlsolve.jl => default_internal_solve.jl} (76%)

diff --git a/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl b/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl
index b4335aa30..f3ff338a2 100644
--- a/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl
+++ b/lib/BoundaryValueDiffEqCore/src/BoundaryValueDiffEqCore.jl
@@ -29,7 +29,7 @@ include("utils.jl")
 include("algorithms.jl")
 include("abstract_types.jl")
 include("alg_utils.jl")
-include("default_nlsolve.jl")
+include("default_internal_solve.jl")
 include("calc_errors.jl")
 
 function SciMLBase.__solve(prob::AbstractBVProblem,
diff --git a/lib/BoundaryValueDiffEqCore/src/algorithms.jl b/lib/BoundaryValueDiffEqCore/src/algorithms.jl
index 0001712f6..7d6cdf384 100644
--- a/lib/BoundaryValueDiffEqCore/src/algorithms.jl
+++ b/lib/BoundaryValueDiffEqCore/src/algorithms.jl
@@ -26,3 +26,12 @@ function Base.show(io::IO, alg::AbstractBoundaryValueDiffEqAlgorithm)
     print(io, join(modifiers, ", "))
     print(io, ")")
 end
+
+# Check what's the internal solver, nonlinear or optimization?
+function __internal_solver(alg::AbstractBoundaryValueDiffEqAlgorithm)
+    # We don't allow both `nlsolve` and `optimize` to be specified at the same time
+    (isnothing(alg.nlsolve) && isnothing(alg.optimize)) &&
+        error("Either `nlsolve` or `optimize` must be specified in the algorithm, but not both.")
+    isnothing(alg.nlsolve) && return alg.optimize
+    isnothing(alg.optimize) && return alg.nlsolve
+end
diff --git a/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl b/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl
similarity index 76%
rename from lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl
rename to lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl
index cb1ebd42a..6ec958cd3 100644
--- a/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl
+++ b/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl
@@ -48,8 +48,20 @@ function __FastShortcutNonlinearPolyalg(::Type{T} = Float64; concrete_jac = noth
     return NonlinearSolvePolyAlgorithm(algs)
 end
 
-@inline __concrete_nonlinearsolve_algorithm(prob, alg) = alg
-@inline function __concrete_nonlinearsolve_algorithm(prob, ::Nothing)
+"""
+    __concrete_solve_algorithm(prob, nlsolve_alg, optimize_alg)
+
+Automatic solver choosing according to the input solver.
+If none of the solvers are specified, we use nonlinear solvers from NonlinearSolve.jl.
+If both of the nonlinear solver and optimization solver are specified, we throw an error.
+If only one of the nonlinear solver and optimization solver is specified, we use that solver.
+"""
+@inline __concrete_solve_algorithm(prob, alg, ::Nothing) = alg
+@inline __concrete_solve_algorithm(prob, ::Nothing, alg) = alg
+@inline __concrete_solve_algorithm(prob,
+    alg1,
+    alg2) = error("Both `nlsolve` and `optimize` are specified in the algorithm, but only one of them is allowed. Please specify only one of them.")
+@inline function __concrete_solve_algorithm(prob, ::Nothing, ::Nothing)
     if prob isa NonlinearLeastSquaresProblem
         return __FastShortcutBVPCompatibleNLLSPolyalg(eltype(prob.u0))
     else
diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 812f39929..2ecd4883f 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -464,6 +464,13 @@ end
     return NonlinearProblem(args...; kwargs...)
 end
 
+# Construct the internal OptimizationProblem
+@inline function __internal_optimization_problem(
+        ::BVProblem{uType, tType, iip}, args...; kwargs...) where {uType, tType, iip}
+    prob = OptimizationProblem(args...; kwargs...)
+    return prob
+end
+
 # Handling Initial Guesses
 """
     __extract_u0(u₀, t₀)
@@ -559,10 +566,17 @@ end
 end
 
 # Construct BVP Solution
-function __build_solution(prob::AbstractBVProblem, odesol, nlsol)
+function __build_solution(
+        prob::AbstractBVProblem, odesol, nlsol::SciMLBase.NonlinearSolution)
     retcode = ifelse(SciMLBase.successful_retcode(nlsol), odesol.retcode, nlsol.retcode)
     return SciMLBase.solution_new_original_retcode(odesol, nlsol, retcode, nlsol.resid)
 end
+function __build_solution(
+        prob::AbstractBVProblem, odesol, optsol::SciMLBase.OptimizationSolution)
+    retcode = ifelse(SciMLBase.successful_retcode(optsol), odesol.retcode, optsol.retcode)
+    return SciMLBase.solution_new_original_retcode(
+        odesol, optsol, retcode, zeros(length(first(odesol)))) # Need a patch in SciMLBase
+end
 
 # Fix3
 @concrete struct __Fix3
@@ -596,3 +610,48 @@ end
 function __split_kwargs(; abstol, adaptive, controller, kwargs...)
     return ((abstol, adaptive, controller), (; abstol, adaptive, kwargs...))
 end
+
+@inline __concrete_kwargs(nlsolve, ::Nothing, nlsolve_kwargs,
+    optimize_kwargs) = (nlsolve_kwargs..., alias_u0 = true)
+@inline __concrete_kwargs(::Nothing, optimize, nlsolve_kwargs, optimize_kwargs) = (;) # Doesn't support for now
+@inline __concrete_kwargs(::Nothing, ::Nothing, nlsolve_kwargs,
+    optimize_kwargs) = (nlsolve_kwargs..., alias_u0 = true)
+
+function __construct_internal_problem(
+        prob::BVProblem, alg, loss, jac, jac_prototype, resid_prototype, y, p, M, N)
+    T = eltype(y)
+    iip = SciMLBase.isinplace(prob)
+    if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize))
+        nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype,
+            jac_prototype = jac_prototype)
+        return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
+    else
+        optf = OptimizationFunction{iip}((x, p) -> 0.0, AutoFiniteDiff(), # Need to investigate the ForwardDiff dual problem
+            cons = loss,
+            cons_j = jac, cons_jac_prototype = jac_prototype)
+        lcons = zeros(T, N*M)
+        ucons = zeros(T, N*M)
+        println("p: ", p)
+        return __internal_optimization_problem(
+            prob, optf, y, p; lcons = lcons, ucons = ucons)
+    end
+end
+
+function __construct_internal_problem(
+        prob::TwoPointBVProblem, alg, loss, jac, jac_prototype, resid_prototype, y, p, M, N)
+    T = eltype(y)
+    iip = SciMLBase.isinplace(prob)
+    if !isnothing(alg.nlsolve)
+        nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype,
+            jac_prototype = jac_prototype)
+        return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
+    else
+        optf = OptimizationFunction{iip}((x, p) -> 0.0, get_dense_ad(diffmode), cons = loss,
+            cons_j = jac, cons_jac_prototype = Matrix(jac_prototype))
+        lcons = zeros(T, N*M)
+        ucons = zeros(T, N*M)
+
+        return __internal_optimization_problem(
+            pro, optf, y, p; lcons = lcons, ucons = ucons)
+    end
+end
diff --git a/lib/BoundaryValueDiffEqMIRK/Project.toml b/lib/BoundaryValueDiffEqMIRK/Project.toml
index 48ba01792..68c2f78d7 100644
--- a/lib/BoundaryValueDiffEqMIRK/Project.toml
+++ b/lib/BoundaryValueDiffEqMIRK/Project.toml
@@ -23,8 +23,8 @@ SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
 Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
-[sources.BoundaryValueDiffEqCore]
-path = "../BoundaryValueDiffEqCore"
+[sources]
+BoundaryValueDiffEqCore = {path = "../BoundaryValueDiffEqCore"}
 
 [compat]
 ADTypes = "1.14"
diff --git a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
index 9d6bb6d7a..bf6d499bf 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
@@ -6,7 +6,7 @@ using BandedMatrices: BandedMatrix, Ones
 using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                AbstractBoundaryValueDiffEqCache, BVPJacobianAlgorithm,
                                recursive_flatten, recursive_flatten!, recursive_unflatten!,
-                               __concrete_nonlinearsolve_algorithm, diff!, EvalSol,
+                               __concrete_solve_algorithm, diff!, EvalSol,
                                concrete_jacobian_algorithm, eval_bc_residual,
                                eval_bc_residual!, get_tmp, __maybe_matmul!, __resize!,
                                __extract_problem_details, __initial_guess, interval,
@@ -14,15 +14,16 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                __restructure_sol, __cache_trait, __get_bcresid_prototype,
                                safe_similar, __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!,
                                recursive_flatten_twopoint!, __internal_nlsolve_problem,
-                               __extract_mesh, __extract_u0, __has_initial_guess,
-                               __initial_guess_length, __initial_guess_on_mesh,
-                               __flatten_initial_guess, __build_solution, __Fix3,
-                               get_dense_ad, _sparse_like, AbstractErrorControl,
-                               DefectControl, GlobalErrorControl, SequentialErrorControl,
-                               HybridErrorControl, HOErrorControl, __use_both_error_control,
-                               __default_coloring_algorithm, DiffCacheNeeded,
-                               NoDiffCacheNeeded, __split_kwargs,
-                               __FastShortcutNonlinearPolyalg
+                               __internal_optimization_problem, __extract_mesh,
+                               __extract_u0, __has_initial_guess, __initial_guess_length,
+                               __initial_guess_on_mesh, __flatten_initial_guess,
+                               __build_solution, __Fix3, get_dense_ad, _sparse_like,
+                               AbstractErrorControl, DefectControl, GlobalErrorControl,
+                               SequentialErrorControl, HybridErrorControl, HOErrorControl,
+                               __use_both_error_control, __default_coloring_algorithm,
+                               DiffCacheNeeded, NoDiffCacheNeeded, __split_kwargs,
+                               __concrete_kwargs, __FastShortcutNonlinearPolyalg,
+                               __construct_internal_problem
 
 using ConcreteStructs: @concrete
 using DiffEqBase: DiffEqBase
@@ -57,110 +58,6 @@ include("interpolation.jl")
 include("mirk_tableaus.jl")
 include("sparse_jacobians.jl")
 
-@setup_workload begin
-    function f1!(du, u, p, t)
-        du[1] = u[2]
-        du[2] = 0
-    end
-    f1 = (u, p, t) -> [u[2], 0]
-
-    function bc1!(residual, u, p, t)
-        residual[1] = u(0.0)[1] - 5
-        residual[2] = u(5.0)[1]
-    end
-
-    bc1 = (u, p, t) -> [u(0.0)[1] - 5, u(5.0)[1]]
-
-    bc1_a! = (residual, ua, p) -> (residual[1] = ua[1] - 5)
-    bc1_b! = (residual, ub, p) -> (residual[1] = ub[1])
-
-    bc1_a = (ua, p) -> [ua[1] - 5]
-    bc1_b = (ub, p) -> [ub[1]]
-
-    tspan = (0.0, 5.0)
-    u0 = [5.0, -3.5]
-    bcresid_prototype = (Array{Float64}(undef, 1), Array{Float64}(undef, 1))
-
-    probs = [BVProblem(f1!, bc1!, u0, tspan; nlls = Val(false)),
-        BVProblem(f1, bc1, u0, tspan; nlls = Val(false)),
-        TwoPointBVProblem(
-            f1!, (bc1_a!, bc1_b!), u0, tspan; bcresid_prototype, nlls = Val(false)),
-        TwoPointBVProblem(
-            f1, (bc1_a, bc1_b), u0, tspan; bcresid_prototype, nlls = Val(false))]
-
-    algs = []
-
-    jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))
-
-    if Preferences.@load_preference("PrecompileMIRK", true)
-        append!(algs, [MIRK2(; jac_alg), MIRK4(; jac_alg), MIRK6(; jac_alg)])
-    end
-
-    @compile_workload begin
-        @sync for prob in probs, alg in algs
-
-            Threads.@spawn solve(prob, alg; dt = 0.2)
-        end
-    end
-
-    f1_nlls! = (du, u, p, t) -> begin
-        du[1] = u[2]
-        du[2] = -u[1]
-    end
-
-    f1_nlls = (u, p, t) -> [u[2], -u[1]]
-
-    bc1_nlls! = (resid, sol, p, t) -> begin
-        solₜ₁ = sol(0.0)
-        solₜ₂ = sol(100.0)
-        resid[1] = solₜ₁[1]
-        resid[2] = solₜ₂[1] - 1
-        resid[3] = solₜ₂[2] + 1.729109
-        return nothing
-    end
-    bc1_nlls = (sol, p, t) -> [sol(0.0)[1], sol(100.0)[1] - 1, sol(100.0)[2] + 1.729109]
-
-    bc1_nlls_a! = (resid, ua, p) -> (resid[1] = ua[1])
-    bc1_nlls_b! = (resid, ub, p) -> (resid[1] = ub[1] - 1; resid[2] = ub[2] + 1.729109)
-
-    bc1_nlls_a = (ua, p) -> [ua[1]]
-    bc1_nlls_b = (ub, p) -> [ub[1] - 1, ub[2] + 1.729109]
-
-    tspan = (0.0, 100.0)
-    u0 = [0.0, 1.0]
-    bcresid_prototype1 = Array{Float64}(undef, 3)
-    bcresid_prototype2 = (Array{Float64}(undef, 1), Array{Float64}(undef, 2))
-
-    probs = [
-        BVProblem(BVPFunction(f1_nlls!, bc1_nlls!; bcresid_prototype = bcresid_prototype1),
-            u0, tspan, nlls = Val(true)),
-        BVProblem(BVPFunction(f1_nlls, bc1_nlls; bcresid_prototype = bcresid_prototype1),
-            u0, tspan, nlls = Val(true)),
-        TwoPointBVProblem(f1_nlls!, (bc1_nlls_a!, bc1_nlls_b!), u0, tspan;
-            bcresid_prototype = bcresid_prototype2, nlls = Val(true)),
-        TwoPointBVProblem(f1_nlls, (bc1_nlls_a, bc1_nlls_b), u0, tspan;
-            bcresid_prototype = bcresid_prototype2, nlls = Val(true))]
-
-    jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))
-
-    nlsolvers = [LevenbergMarquardt(; disable_geodesic = Val(true)), GaussNewton()]
-
-    algs = []
-
-    if Preferences.@load_preference("PrecompileMIRKNLLS", false)
-        for nlsolve in nlsolvers
-            append!(algs, [MIRK2(; jac_alg, nlsolve), MIRK6(; jac_alg, nlsolve)])
-        end
-    end
-
-    @compile_workload begin
-        @sync for prob in probs, alg in algs
-
-            Threads.@spawn solve(prob, alg; dt = 0.2, abstol = 1e-2)
-        end
-    end
-end
-
 export MIRK2, MIRK3, MIRK4, MIRK5, MIRK6, MIRK6I
 export BVPJacobianAlgorithm
 export maxsol, minsol
diff --git a/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl b/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl
index 996d161ca..44a3a11a9 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl
@@ -17,6 +17,9 @@ for order in (2, 3, 4, 5, 6)
             `NonlinearProblem` interface can be used. Note that any autodiff argument for
             the solver will be ignored and a custom jacobian algorithm will be used.
 
+          - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML
+            `OptimizationProblem` interface can be used. Note that any autodiff argument for
+            the solver will be ignored and a custom jacobian algorithm will be used.
           - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to
             `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to
             use based on the input types and problem type.
@@ -48,8 +51,9 @@ for order in (2, 3, 4, 5, 6)
         }
         ```
         """
-        @kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractMIRK
+        @kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractMIRK
             nlsolve::N = nothing
+            optimize::O = nothing
             jac_alg::J = BVPJacobianAlgorithm()
             defect_threshold::T = 0.1
             max_num_subintervals::Int = 3000
@@ -73,6 +77,9 @@ for order in (6)
             `NonlinearProblem` interface can be used. Note that any autodiff argument for
             the solver will be ignored and a custom jacobian algorithm will be used.
 
+          - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML
+            `OptimizationProblem` interface can be used. Note that any autodiff argument for
+            the solver will be ignored and a custom jacobian algorithm will be used.
           - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to
             `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to
             use based on the input types and problem type.
@@ -104,8 +111,9 @@ for order in (6)
         }
         ```
         """
-        @kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractMIRK
+        @kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractMIRK
             nlsolve::N = nothing
+            optimize::O = nothing
             jac_alg::J = BVPJacobianAlgorithm()
             defect_threshold::T = 0.1
             max_num_subintervals::Int = 3000
diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
index 336f6788c..fd714b087 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
@@ -28,6 +28,7 @@
     new_stages
     resid_size
     nlsolve_kwargs
+    optimize_kwargs
     kwargs
 end
 
@@ -35,7 +36,8 @@ Base.eltype(::MIRKCache{iip, T, use_both}) where {iip, T, use_both} = T
 
 function SciMLBase.__init(
         prob::BVProblem, alg::AbstractMIRK; dt = 0.0, abstol = 1e-6, adaptive = true,
-        controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), kwargs...)
+        controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol),
+        optimize_kwargs = (; abstol = abstol), kwargs...)
     @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
     iip = isinplace(prob)
     diffcache = __cache_trait(alg.jac_alg)
@@ -127,9 +129,9 @@ function SciMLBase.__init(
 
     return MIRKCache{iip, T, use_both, typeof(diffcache), fit_parameters}(
         alg_order(alg), stage, N, size(X), f, bc, prob_, prob.problem_type, prob.p,
-        alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete, k_interp,
-        y, y₀, residual, fᵢ_cache, fᵢ₂_cache, errors, new_stages, resid₁_size,
-        nlsolve_kwargs, (; abstol, dt, adaptive, controller, fit_parameters, kwargs...))
+        alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete, k_interp, y, y₀,
+        residual, fᵢ_cache, fᵢ₂_cache, errors, new_stages, resid₁_size, nlsolve_kwargs,
+        optimize_kwargs, (; abstol, dt, adaptive, controller, fit_parameters, kwargs...))
 end
 
 """
@@ -187,9 +189,11 @@ end
 
 function __perform_mirk_iteration(
         cache::MIRKCache, abstol, adaptive::Bool, controller::AbstractErrorControl)
-    nlprob = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀))
-    nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve)
-    sol_nlprob = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true)
+    nlprob = __construct_problem(cache, vec(cache.y₀), copy(cache.y₀))
+    solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize)
+    kwargs = __concrete_kwargs(
+        cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs)
+    sol_nlprob = solve(nlprob, solve_alg, kwargs...)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
     error_norm = 2 * abstol
@@ -202,7 +206,7 @@ function __perform_mirk_iteration(
     if info == ReturnCode.Success # Nonlinear Solve was successful
         error_norm,
         info = error_estimate!(
-            cache, controller, cache.errors, sol_nlprob, nlsolve_alg, abstol)
+            cache, controller, cache.errors, sol_nlprob, solve_alg, abstol)
     end
 
     if info == ReturnCode.Success # Nonlinear Solve Successful and defect norm is acceptable
@@ -235,7 +239,7 @@ function __perform_mirk_iteration(
 end
 
 # Constructing the Nonlinear Problem
-function __construct_nlproblem(
+function __construct_problem(
         cache::MIRKCache{iip}, y::AbstractVector, y₀::AbstractVectorOfArray) where {iip}
     pt = cache.problem_type
     (; jac_alg) = cache.alg
@@ -275,7 +279,7 @@ function __construct_nlproblem(
             u, p, cache.y, pt, cache.bc, cache.mesh, cache, eval_sol, trait)
     end
 
-    return __construct_nlproblem(cache, y, loss_bc, loss_collocation, loss, pt)
+    return __construct_problem(cache, y, loss_bc, loss_collocation, loss, pt)
 end
 
 @views function __mirk_loss!(resid, u, p, y, pt::StandardBVProblem, bc!::BC, residual,
@@ -381,7 +385,7 @@ end
     return mapreduce(vec, vcat, resids)
 end
 
-function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collocation::C,
+function __construct_problem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collocation::C,
         loss::LF, ::StandardBVProblem) where {iip, BC, C, LF}
     (; jac_alg) = cache.alg
     (; bc_diffmode) = jac_alg
@@ -457,10 +461,8 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo
     end
 
     resid_prototype = vcat(resid_bc, resid_collocation)
-    nlf = NonlinearFunction{iip}(
-        loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
-
-    return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p)
+    return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype,
+        resid_prototype, y, cache.p, cache.M, N)
 end
 
 function __mirk_mpoint_jacobian!(
@@ -505,9 +507,9 @@ function __mirk_mpoint_jacobian(
     return J
 end
 
-function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collocation::C,
+function __construct_problem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collocation::C,
         loss::LF, ::TwoPointBVProblem) where {iip, BC, C, LF}
-    (; nlsolve, jac_alg) = cache.alg
+    (; jac_alg) = cache.alg
     N = length(cache.mesh)
 
     resid = vcat(@view(cache.bcresid_prototype[1:prod(cache.resid_size[1])]),
@@ -548,9 +550,8 @@ function __construct_nlproblem(cache::MIRKCache{iip}, y, loss_bc::BC, loss_collo
     end
 
     resid_prototype = copy(resid)
-    nlf = NonlinearFunction{iip}(
-        loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
-    return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p)
+    return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype,
+        resid_prototype, y, cache.p, cache.M, N)
 end
 
 function __mirk_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid, p) where {L}
diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
index 471df76af..99f053210 100644
--- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
@@ -462,3 +462,23 @@ end
 
     @test sol.prob.p ≈ [17.09658] atol=1e-5
 end
+
+@testset "Convergence with optimization based solver" setup=[MIRKConvergenceTests] begin
+    using LinearAlgebra, DiffEqDevTools, OptimizationMOI, Ipopt
+
+    @testset "Problem: $i" for i in (3, 4, 5, 6, 9, 10)
+        prob = probArr[i]
+        @testset "MIRK$order" for (_, order) in enumerate((2, 3, 4, 5, 6))
+            sim = test_convergence(
+                dts, prob, mirk_solver(Val(order), optimize = Ipopt.Optimizer());
+                abstol = 1e-8, reltol = 1e-8)
+            @test sim.𝒪est[:final]≈order atol=testTol
+        end
+
+        @testset "MIRK$(order)I" for (_, order) in enumerate((6,))
+            sim = test_convergence(dts, prob, MIRK6I(; optimize = Ipopt.Optimizer());
+                abstol = 1e-8, reltol = 1e-8)
+            @test sim.𝒪est[:final]≈order atol=testTol
+        end
+    end
+end

From 9b7ff46304c250974cbe2d5fb71ce834138f41c5 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Tue, 15 Jul 2025 00:47:56 +0800
Subject: [PATCH 02/50] Fix some typos

---
 .../src/BoundaryValueDiffEqMIRK.jl            | 104 ++++++++++++++++++
 1 file changed, 104 insertions(+)

diff --git a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
index bf6d499bf..b14341097 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
@@ -58,6 +58,110 @@ include("interpolation.jl")
 include("mirk_tableaus.jl")
 include("sparse_jacobians.jl")
 
+@setup_workload begin
+    function f1!(du, u, p, t)
+        du[1] = u[2]
+        du[2] = 0
+    end
+    f1 = (u, p, t) -> [u[2], 0]
+
+    function bc1!(residual, u, p, t)
+        residual[1] = u(0.0)[1] - 5
+        residual[2] = u(5.0)[1]
+    end
+
+    bc1 = (u, p, t) -> [u(0.0)[1] - 5, u(5.0)[1]]
+
+    bc1_a! = (residual, ua, p) -> (residual[1] = ua[1] - 5)
+    bc1_b! = (residual, ub, p) -> (residual[1] = ub[1])
+
+    bc1_a = (ua, p) -> [ua[1] - 5]
+    bc1_b = (ub, p) -> [ub[1]]
+
+    tspan = (0.0, 5.0)
+    u0 = [5.0, -3.5]
+    bcresid_prototype = (Array{Float64}(undef, 1), Array{Float64}(undef, 1))
+
+    probs = [BVProblem(f1!, bc1!, u0, tspan; nlls = Val(false)),
+        BVProblem(f1, bc1, u0, tspan; nlls = Val(false)),
+        TwoPointBVProblem(
+            f1!, (bc1_a!, bc1_b!), u0, tspan; bcresid_prototype, nlls = Val(false)),
+        TwoPointBVProblem(
+            f1, (bc1_a, bc1_b), u0, tspan; bcresid_prototype, nlls = Val(false))]
+
+    algs = []
+
+    jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))
+
+    if Preferences.@load_preference("PrecompileMIRK", true)
+        append!(algs, [MIRK2(; jac_alg), MIRK4(; jac_alg), MIRK6(; jac_alg)])
+    end
+
+    @compile_workload begin
+        @sync for prob in probs, alg in algs
+
+            Threads.@spawn solve(prob, alg; dt = 0.2)
+        end
+    end
+
+    f1_nlls! = (du, u, p, t) -> begin
+        du[1] = u[2]
+        du[2] = -u[1]
+    end
+
+    f1_nlls = (u, p, t) -> [u[2], -u[1]]
+
+    bc1_nlls! = (resid, sol, p, t) -> begin
+        solₜ₁ = sol(0.0)
+        solₜ₂ = sol(100.0)
+        resid[1] = solₜ₁[1]
+        resid[2] = solₜ₂[1] - 1
+        resid[3] = solₜ₂[2] + 1.729109
+        return nothing
+    end
+    bc1_nlls = (sol, p, t) -> [sol(0.0)[1], sol(100.0)[1] - 1, sol(100.0)[2] + 1.729109]
+
+    bc1_nlls_a! = (resid, ua, p) -> (resid[1] = ua[1])
+    bc1_nlls_b! = (resid, ub, p) -> (resid[1] = ub[1] - 1; resid[2] = ub[2] + 1.729109)
+
+    bc1_nlls_a = (ua, p) -> [ua[1]]
+    bc1_nlls_b = (ub, p) -> [ub[1] - 1, ub[2] + 1.729109]
+
+    tspan = (0.0, 100.0)
+    u0 = [0.0, 1.0]
+    bcresid_prototype1 = Array{Float64}(undef, 3)
+    bcresid_prototype2 = (Array{Float64}(undef, 1), Array{Float64}(undef, 2))
+
+    probs = [
+        BVProblem(BVPFunction(f1_nlls!, bc1_nlls!; bcresid_prototype = bcresid_prototype1),
+            u0, tspan, nlls = Val(true)),
+        BVProblem(BVPFunction(f1_nlls, bc1_nlls; bcresid_prototype = bcresid_prototype1),
+            u0, tspan, nlls = Val(true)),
+        TwoPointBVProblem(f1_nlls!, (bc1_nlls_a!, bc1_nlls_b!), u0, tspan;
+            bcresid_prototype = bcresid_prototype2, nlls = Val(true)),
+        TwoPointBVProblem(f1_nlls, (bc1_nlls_a, bc1_nlls_b), u0, tspan;
+            bcresid_prototype = bcresid_prototype2, nlls = Val(true))]
+
+    jac_alg = BVPJacobianAlgorithm(AutoForwardDiff(; chunksize = 2))
+
+    nlsolvers = [LevenbergMarquardt(; disable_geodesic = Val(true)), GaussNewton()]
+
+    algs = []
+
+    if Preferences.@load_preference("PrecompileMIRKNLLS", false)
+        for nlsolve in nlsolvers
+            append!(algs, [MIRK2(; jac_alg, nlsolve), MIRK6(; jac_alg, nlsolve)])
+        end
+    end
+
+    @compile_workload begin
+        @sync for prob in probs, alg in algs
+
+            Threads.@spawn solve(prob, alg; dt = 0.2, abstol = 1e-2)
+        end
+    end
+end
+
 export MIRK2, MIRK3, MIRK4, MIRK5, MIRK6, MIRK6I
 export BVPJacobianAlgorithm
 export maxsol, minsol

From 1cc250dc530867b14fa61c5a967e10d2839ff923 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Tue, 15 Jul 2025 01:56:53 +0800
Subject: [PATCH 03/50] Don't fail other solvers

---
 .../src/BoundaryValueDiffEqAscher.jl                      | 2 +-
 lib/BoundaryValueDiffEqAscher/src/ascher.jl               | 2 +-
 lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl | 8 ++++++++
 .../src/BoundaryValueDiffEqFIRK.jl                        | 2 +-
 lib/BoundaryValueDiffEqFIRK/src/firk.jl                   | 2 +-
 lib/BoundaryValueDiffEqMIRK/Project.toml                  | 4 +++-
 lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl      | 2 +-
 .../src/BoundaryValueDiffEqMIRKN.jl                       | 6 +++---
 lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl                 | 2 +-
 .../src/BoundaryValueDiffEqShooting.jl                    | 2 +-
 lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl  | 2 +-
 lib/BoundaryValueDiffEqShooting/src/single_shooting.jl    | 2 +-
 12 files changed, 23 insertions(+), 13 deletions(-)

diff --git a/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl b/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl
index 1c2725eae..e7da643e5 100644
--- a/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl
+++ b/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl
@@ -6,7 +6,7 @@ using AlmostBlockDiagonals: AlmostBlockDiagonals, IntermediateAlmostBlockDiagona
 using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                AbstractBoundaryValueDiffEqCache, BVPJacobianAlgorithm,
                                __extract_problem_details, concrete_jacobian_algorithm,
-                               __Fix3, __concrete_nonlinearsolve_algorithm,
+                               __Fix3, __concrete_solve_algorithm,
                                __internal_nlsolve_problem, __vec, __vec_f, __vec_f!,
                                __vec_bc, __vec_bc!, __extract_mesh, get_dense_ad,
                                __get_bcresid_prototype, __split_kwargs,
diff --git a/lib/BoundaryValueDiffEqAscher/src/ascher.jl b/lib/BoundaryValueDiffEqAscher/src/ascher.jl
index 4c51a956c..5c280d26c 100644
--- a/lib/BoundaryValueDiffEqAscher/src/ascher.jl
+++ b/lib/BoundaryValueDiffEqAscher/src/ascher.jl
@@ -176,7 +176,7 @@ function __perform_ascher_iteration(
         cache::AscherCache{iip, T}, abstol, adaptive::Bool) where {iip, T}
     info::ReturnCode.T = ReturnCode.Success
     nlprob = __construct_nlproblem(cache)
-    nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve)
+    nlsolve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve)
     nlsol = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs...)
     error_norm = 2 * abstol
     info = nlsol.retcode
diff --git a/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl b/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl
index 6ec958cd3..cdbfbddc9 100644
--- a/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl
+++ b/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl
@@ -56,11 +56,19 @@ If none of the solvers are specified, we use nonlinear solvers from NonlinearSol
 If both of the nonlinear solver and optimization solver are specified, we throw an error.
 If only one of the nonlinear solver and optimization solver is specified, we use that solver.
 """
+@inline __concrete_solve_algorithm(prob, alg) = alg
 @inline __concrete_solve_algorithm(prob, alg, ::Nothing) = alg
 @inline __concrete_solve_algorithm(prob, ::Nothing, alg) = alg
 @inline __concrete_solve_algorithm(prob,
     alg1,
     alg2) = error("Both `nlsolve` and `optimize` are specified in the algorithm, but only one of them is allowed. Please specify only one of them.")
+@inline function __concrete_solve_algorithm(prob, ::Nothing)
+    if prob isa NonlinearLeastSquaresProblem
+        return __FastShortcutBVPCompatibleNLLSPolyalg(eltype(prob.u0))
+    else
+        return __FastShortcutBVPCompatibleNonlinearPolyalg(eltype(prob.u0))
+    end
+end
 @inline function __concrete_solve_algorithm(prob, ::Nothing, ::Nothing)
     if prob isa NonlinearLeastSquaresProblem
         return __FastShortcutBVPCompatibleNLLSPolyalg(eltype(prob.u0))
diff --git a/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl b/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
index 824df0bb6..f906e6cf1 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
@@ -6,7 +6,7 @@ using BandedMatrices: BandedMatrix, Ones
 using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                AbstractBoundaryValueDiffEqCache, BVPJacobianAlgorithm,
                                recursive_flatten, recursive_flatten!, recursive_unflatten!,
-                               __concrete_nonlinearsolve_algorithm, diff!, EvalSol,
+                               __concrete_solve_algorithm, diff!, EvalSol,
                                concrete_jacobian_algorithm, eval_bc_residual, interval,
                                eval_bc_residual!, get_tmp, __maybe_matmul!, __resize!,
                                __extract_problem_details, __initial_guess, nodual_value,
diff --git a/lib/BoundaryValueDiffEqFIRK/src/firk.jl b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
index 9dd12493d..9df4ce4bd 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/firk.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
@@ -404,7 +404,7 @@ end
 function __perform_firk_iteration(
         cache::Union{FIRKCacheExpand, FIRKCacheNested}, abstol, adaptive::Bool)
     nlprob = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀))
-    nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve)
+    nlsolve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve)
     sol_nlprob = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
diff --git a/lib/BoundaryValueDiffEqMIRK/Project.toml b/lib/BoundaryValueDiffEqMIRK/Project.toml
index 68c2f78d7..818afe79b 100644
--- a/lib/BoundaryValueDiffEqMIRK/Project.toml
+++ b/lib/BoundaryValueDiffEqMIRK/Project.toml
@@ -67,9 +67,11 @@ DiffEqDevTools = "f3b72e0c-5b89-59e1-b016-84e28bfd966d"
 Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
 Hwloc = "0e44f5e4-bd66-52a0-8798-143a42290a1d"
 InteractiveUtils = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
+Ipopt = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
 JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b"
 LinearSolve = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
 Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
+OptimizationMOI = "fd9f6733-72f4-499f-8506-86b2bdd0dea1"
 OrdinaryDiffEqRosenbrock = "43230ef6-c299-4910-a778-202eb28ce4ce"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 ReTestItems = "817f1d60-ba6b-4fd5-9520-3cf149f6a823"
@@ -78,4 +80,4 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Aqua", "DiffEqDevTools", "Enzyme", "Hwloc", "InteractiveUtils", "JET", "LinearSolve", "Mooncake", "OrdinaryDiffEqRosenbrock", "Random", "ReTestItems", "RecursiveArrayTools", "StaticArrays", "Test"]
+test = ["Aqua", "DiffEqDevTools", "Enzyme", "Hwloc", "InteractiveUtils", "Ipopt", "JET", "LinearSolve", "Mooncake", "OptimizationMOI", "OrdinaryDiffEqRosenbrock", "Random", "ReTestItems", "RecursiveArrayTools", "StaticArrays", "Test"]
diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
index 99f053210..75fcb840e 100644
--- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
@@ -463,7 +463,7 @@ end
     @test sol.prob.p ≈ [17.09658] atol=1e-5
 end
 
-@testset "Convergence with optimization based solver" setup=[MIRKConvergenceTests] begin
+@testitem "Convergence with optimization based solver" setup=[MIRKConvergenceTests] begin
     using LinearAlgebra, DiffEqDevTools, OptimizationMOI, Ipopt
 
     @testset "Problem: $i" for i in (3, 4, 5, 6, 9, 10)
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl b/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
index 234b3f05d..702ecda96 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
@@ -6,9 +6,9 @@ using BandedMatrices: BandedMatrix, Ones
 using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                AbstractBoundaryValueDiffEqCache, BVPJacobianAlgorithm,
                                recursive_flatten, recursive_flatten!, recursive_unflatten!,
-                               __concrete_nonlinearsolve_algorithm, diff!, EvalSol,
-                               eval_bc_residual, eval_bc_residual!, get_tmp,
-                               __maybe_matmul!, __extract_problem_details, __initial_guess,
+                               __concrete_solve_algorithm, diff!, EvalSol, eval_bc_residual,
+                               eval_bc_residual!, get_tmp, __maybe_matmul!,
+                               __extract_problem_details, __initial_guess,
                                __maybe_allocate_diffcache, __restructure_sol,
                                __get_bcresid_prototype, safe_similar, __vec, __vec_f,
                                __vec_f!, __vec_bc, __vec_bc!, __vec_so_bc!, __vec_so_bc,
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
index ffd52d7e1..425ae8e89 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
@@ -113,7 +113,7 @@ end
 
 function __perform_mirkn_iteration(cache::MIRKNCache)
     nlprob::NonlinearProblem = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀))
-    nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve)
+    nlsolve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve)
     sol_nlprob = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
diff --git a/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl b/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
index 9292cc9f5..7b0169f20 100644
--- a/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
@@ -5,7 +5,7 @@ using ArrayInterface: fast_scalar_indexing
 using BandedMatrices: BandedMatrix, Ones
 using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm, BVPJacobianAlgorithm,
                                recursive_flatten, recursive_flatten!, recursive_unflatten!,
-                               __concrete_nonlinearsolve_algorithm, diff!, __any_sparse_ad,
+                               __concrete_solve_algorithm, diff!, __any_sparse_ad,
                                __cache_trait, concrete_jacobian_algorithm, eval_bc_residual,
                                eval_bc_residual!, get_tmp, __maybe_matmul!,
                                __extract_problem_details, __initial_guess,
diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
index 3fb8a958f..d68a9c2a9 100644
--- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
@@ -142,7 +142,7 @@ function __solve_nlproblem!(
 
     # NOTE: u_at_nodes is updated inplace
     nlprob = __internal_nlsolve_problem(prob, M, N, loss_function!, u_at_nodes, prob.p)
-    nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, alg.nlsolve)
+    nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve)
     __solve(nlprob, nlsolve_alg; kwargs..., alias_u0 = true)
 
     return nothing
diff --git a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
index 311f16a28..1a03d028f 100644
--- a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
@@ -76,7 +76,7 @@ function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (;
     nlf = NonlinearFunction{iip}(loss_fn; jac_prototype = jac_prototype,
         resid_prototype = resid_prototype, jac = jac_fn)
     nlprob = __internal_nlsolve_problem(prob, resid_prototype, u0, nlf, vec(u0), prob.p)
-    nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, alg.nlsolve)
+    nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve)
     nlsol::SciMLBase.NonlinearSolution = __solve(
         nlprob, nlsolve_alg; nlsolve_kwargs..., verbose, kwargs...)
 

From 59314300cc1325ffb0ded20803c95e86987b5c70 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Tue, 15 Jul 2025 02:18:48 +0800
Subject: [PATCH 04/50] Fix unrelated failings

---
 lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl            | 8 ++++----
 lib/BoundaryValueDiffEqFIRK/src/collocation.jl           | 8 ++++----
 lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl | 2 +-
 3 files changed, 9 insertions(+), 9 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl b/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl
index 15cb25bcf..04e3c3c7b 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl
@@ -106,7 +106,7 @@ end
     end
     τ = (t - mesh[j])
 
-    nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve)
+    nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve)
     nestprob_p = zeros(T, cache.M + 2)
 
     yᵢ = copy(cache.y[j].du)
@@ -152,7 +152,7 @@ end
     end
     τ = (t - mesh[j])
 
-    nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve)
+    nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve)
     nestprob_p = zeros(T, cache.M + 2)
 
     yᵢ = copy(cache.y[j])
@@ -457,7 +457,7 @@ end
     (; f, mesh, mesh_dt, defect, ITU, nest_prob, alg) = cache
     (; q_coeff, τ_star) = ITU
 
-    nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, cache.alg.nlsolve)
+    nlsolve_alg = __concrete_solve_algorithm(nest_prob, cache.alg.nlsolve)
     nestprob_p = zeros(T, cache.M + 2)
 
     for i in 1:(length(mesh) - 1)
@@ -509,7 +509,7 @@ end
     (; f, mesh, mesh_dt, defect, ITU, nest_prob, alg) = cache
     (; q_coeff, τ_star) = ITU
 
-    nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, cache.alg.nlsolve)
+    nlsolve_alg = __concrete_solve_algorithm(nest_prob, cache.alg.nlsolve)
     nestprob_p = zeros(T, cache.M + 2)
 
     for i in 1:(length(mesh) - 1)
diff --git a/lib/BoundaryValueDiffEqFIRK/src/collocation.jl b/lib/BoundaryValueDiffEqFIRK/src/collocation.jl
index 0c60b7dfa..bbdcf1350 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/collocation.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/collocation.jl
@@ -121,7 +121,7 @@ end
 
     T = eltype(u)
     nestprob_p = vcat(T(mesh[1]), T(mesh_dt[1]), get_tmp(y[1], u))
-    nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve)
+    nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve)
 
     for i in eachindex(k_discrete)
         residᵢ = residual[i]
@@ -151,7 +151,7 @@ end
 
     T = eltype(u)
     nestprob_p = vcat(T(mesh[1]), T(mesh_dt[1]), y[1])
-    nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, cache.alg.nlsolve)
+    nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, cache.alg.nlsolve)
 
     for i in eachindex(k_discrete)
         residᵢ = residual[i]
@@ -265,7 +265,7 @@ end
 
     T = eltype(u)
     nestprob_p = vcat(T(mesh[1]), T(mesh_dt[1]), get_tmp(y[1], u))
-    nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve)
+    nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve)
 
     for i in eachindex(k_discrete)
         residᵢ = residuals[i]
@@ -296,7 +296,7 @@ end
 
     T = eltype(u)
     nestprob_p = vcat(T(mesh[1]), T(mesh_dt[1]), y[1])
-    nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve)
+    nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve)
 
     for i in eachindex(k_discrete)
         residᵢ = residuals[i]
diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
index d68a9c2a9..0f5be76d3 100644
--- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
@@ -218,7 +218,7 @@ function __solve_nlproblem!(::StandardBVProblem, alg::MultipleShooting, bcresid_
 
     # NOTE: u_at_nodes is updated inplace
     nlprob = __internal_nlsolve_problem(prob, M, N, loss_function!, u_at_nodes, prob.p)
-    nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, alg.nlsolve)
+    nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve)
     __solve(nlprob, nlsolve_alg; kwargs..., alias_u0 = true)
 
     return nothing

From fac21aaa57c1281ddc34d575090582bf8a9c2621 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Tue, 15 Jul 2025 09:52:25 +0800
Subject: [PATCH 05/50] Fix more concrete solve algorithm

---
 lib/BoundaryValueDiffEqFIRK/src/interpolation.jl | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl b/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl
index 6de15e1a2..7fc755328 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/interpolation.jl
@@ -95,7 +95,7 @@ end
     τ = (t - mesh[j])
     length_z = length(z)
 
-    nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve)
+    nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve)
     nestprob_p = zeros(T, cache.M + 2)
 
     yᵢ = copy(cache.y[j].du)
@@ -141,7 +141,7 @@ end
     τ = (t - mesh[j])
     length_dz = length(dz)
 
-    nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve)
+    nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve)
     nestprob_p = zeros(T, cache.M + 2)
 
     yᵢ = copy(cache.y[j].du)
@@ -413,7 +413,7 @@ function (s::EvalSol{C})(tval::Number) where {C <: FIRKCacheNested}
     h = mesh_dt[j]
     τ = tval - t[j]
 
-    nest_nlsolve_alg = __concrete_nonlinearsolve_algorithm(nest_prob, alg.nlsolve)
+    nest_nlsolve_alg = __concrete_solve_algorithm(nest_prob, alg.nlsolve)
     nestprob_p = zeros(cache.M + 2)
 
     yᵢ = u[j]

From f9cd33fb2259d8d8a1ff7dd16e90da16df9ac285 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Thu, 17 Jul 2025 20:10:18 +0800
Subject: [PATCH 06/50] Add default cost

---
 lib/BoundaryValueDiffEqCore/src/utils.jl | 10 +++++++---
 1 file changed, 7 insertions(+), 3 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 2ecd4883f..777a87876 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -617,6 +617,10 @@ end
 @inline __concrete_kwargs(::Nothing, ::Nothing, nlsolve_kwargs,
     optimize_kwargs) = (nlsolve_kwargs..., alias_u0 = true)
 
+@inline __default_cost(::Nothing) = (x, p) -> 0.0
+@inline __default_cost(f) = f
+@inline __default_cost(fun::BVPFunction) = __default_cost(fun.cost)
+
 function __construct_internal_problem(
         prob::BVProblem, alg, loss, jac, jac_prototype, resid_prototype, y, p, M, N)
     T = eltype(y)
@@ -626,12 +630,11 @@ function __construct_internal_problem(
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
-        optf = OptimizationFunction{iip}((x, p) -> 0.0, AutoFiniteDiff(), # Need to investigate the ForwardDiff dual problem
+        optf = OptimizationFunction{iip}(__default_cost(prob.f), AutoFiniteDiff(), # Need to investigate the ForwardDiff dual problem
             cons = loss,
             cons_j = jac, cons_jac_prototype = jac_prototype)
         lcons = zeros(T, N*M)
         ucons = zeros(T, N*M)
-        println("p: ", p)
         return __internal_optimization_problem(
             prob, optf, y, p; lcons = lcons, ucons = ucons)
     end
@@ -646,7 +649,8 @@ function __construct_internal_problem(
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
-        optf = OptimizationFunction{iip}((x, p) -> 0.0, get_dense_ad(diffmode), cons = loss,
+        optf = OptimizationFunction{iip}(
+            __default_cost(prob.f), get_dense_ad(diffmode), cons = loss,
             cons_j = jac, cons_jac_prototype = Matrix(jac_prototype))
         lcons = zeros(T, N*M)
         ucons = zeros(T, N*M)

From c694c26ff7c651de979a3749ab91288aa20ff353 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Fri, 18 Jul 2025 02:06:02 +0800
Subject: [PATCH 07/50] Nested FIRK should work with optimization solvers

---
 lib/BoundaryValueDiffEqCore/src/utils.jl      |  2 +
 .../src/BoundaryValueDiffEqFIRK.jl            |  9 ++-
 lib/BoundaryValueDiffEqFIRK/src/algorithms.jl | 45 +++++++++--
 lib/BoundaryValueDiffEqFIRK/src/firk.jl       | 74 ++++++++++---------
 lib/BoundaryValueDiffEqMIRK/src/algorithms.jl |  3 +-
 lib/BoundaryValueDiffEqMIRK/src/mirk.jl       |  2 +
 6 files changed, 86 insertions(+), 49 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 777a87876..76de3867f 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -617,6 +617,8 @@ end
 @inline __concrete_kwargs(::Nothing, ::Nothing, nlsolve_kwargs,
     optimize_kwargs) = (nlsolve_kwargs..., alias_u0 = true)
 
+## Optimization solver related utils ##
+
 @inline __default_cost(::Nothing) = (x, p) -> 0.0
 @inline __default_cost(f) = f
 @inline __default_cost(fun::BVPFunction) = __default_cost(fun.cost)
diff --git a/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl b/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
index f906e6cf1..7d6b6f916 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
@@ -12,7 +12,7 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                __extract_problem_details, __initial_guess, nodual_value,
                                __maybe_allocate_diffcache, __restructure_sol,
                                __get_bcresid_prototype, __vec, __vec_f, __vec_f!, __vec_bc,
-                               __vec_bc!, recursive_flatten_twopoint!,
+                               __vec_bc!, recursive_flatten_twopoint!, __concrete_kwargs,
                                __internal_nlsolve_problem, __extract_mesh, __extract_u0,
                                __default_coloring_algorithm, __maybe_allocate_diffcache,
                                __restructure_sol, __get_bcresid_prototype, safe_similar,
@@ -20,9 +20,10 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                recursive_flatten_twopoint!, __internal_nlsolve_problem,
                                __extract_mesh, __extract_u0, DiffCacheNeeded,
                                NoDiffCacheNeeded, __has_initial_guess,
-                               __initial_guess_length, __initial_guess_on_mesh,
-                               __flatten_initial_guess, __build_solution, __Fix3,
-                               __split_kwargs, _sparse_like, get_dense_ad
+                               __construct_internal_problem, __initial_guess_length,
+                               __initial_guess_on_mesh, __flatten_initial_guess,
+                               __build_solution, __Fix3, __split_kwargs, _sparse_like,
+                               get_dense_ad, __internal_optimization_problem
 
 using ConcreteStructs: @concrete
 using DiffEqBase: DiffEqBase
diff --git a/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl b/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl
index c8156a868..1e093a45f 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl
@@ -17,6 +17,9 @@ for stage in (1, 2, 3, 5, 7)
             `NonlinearProblem` interface can be used. Note that any autodiff argument for
             the solver will be ignored and a custom jacobian algorithm will be used.
 
+          - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML
+            `OptimizationProblem` interface can be used. Note that any autodiff argument for
+            the solver will be ignored and a custom jacobian algorithm will be used.
           - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to
             `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to
             use based on the input types and problem type.
@@ -87,8 +90,9 @@ for stage in (1, 2, 3, 5, 7)
         }
         ```
         """
-        Base.@kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK
+        Base.@kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK
             nlsolve::N = nothing
+            optimize::O = nothing
             jac_alg::J = BVPJacobianAlgorithm()
             nested_nlsolve::Bool = false
             nested_nlsolve_kwargs::NamedTuple = (;)
@@ -96,13 +100,16 @@ for stage in (1, 2, 3, 5, 7)
             max_num_subintervals::Int = 3000
         end
         $(alg)(nlsolve::N,
+            optimize::O,
             jac_alg::J;
             nested = false,
             nested_nlsolve_kwargs::NamedTuple = (;),
             defect_threshold::T = 0.1,
             max_num_subintervals::Int = 3000) where {N,
+            O,
             J,
-            T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
+            T} = $(alg){N, O, J, T}(
+            nlsolve, optimize, jac_alg, nested, nested_nlsolve_kwargs,
             defect_threshold, max_num_subintervals)
     end
 end
@@ -123,6 +130,10 @@ for stage in (2, 3, 4, 5)
             `NonlinearProblem` interface can be used. Note that any autodiff argument for
               the solver will be ignored and a custom jacobian algorithm will be used.
 
+            - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML
+            `OptimizationProblem` interface can be used. Note that any autodiff argument for
+            the solver will be ignored and a custom jacobian algorithm will be used.
+
           - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to
             `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to
             use based on the input types and problem type.
@@ -192,8 +203,9 @@ for stage in (2, 3, 4, 5)
               }
               ```
         """
-        Base.@kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK
+        Base.@kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK
             nlsolve::N = nothing
+            optimize::O = nothing
             jac_alg::J = BVPJacobianAlgorithm()
             nested_nlsolve::Bool = false
             nested_nlsolve_kwargs::NamedTuple = (;)
@@ -201,13 +213,16 @@ for stage in (2, 3, 4, 5)
             max_num_subintervals::Int = 3000
         end
         $(alg)(nlsolve::N,
+            optimize::O,
             jac_alg::J;
             nested = false,
             nested_nlsolve_kwargs::NamedTuple = (;),
             defect_threshold::T = 0.1,
             max_num_subintervals::Int = 3000) where {N,
+            O,
             J,
-            T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
+            T} = $(alg){N, O, J, T}(
+            nlsolve, optimize, jac_alg, nested, nested_nlsolve_kwargs,
             defect_threshold, max_num_subintervals)
     end
 end
@@ -228,6 +243,9 @@ for stage in (2, 3, 4, 5)
             `NonlinearProblem` interface can be used. Note that any autodiff argument for
             the solver will be ignored and a custom jacobian algorithm will be used.
 
+          - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML
+            `OptimizationProblem` interface can be used. Note that any autodiff argument for
+            the solver will be ignored and a custom jacobian algorithm will be used.
           - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to
             `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to
             use based on the input types and problem type.
@@ -299,8 +317,9 @@ for stage in (2, 3, 4, 5)
         }
         ```
         """
-        Base.@kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK
+        Base.@kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK
             nlsolve::N = nothing
+            optimize::O = nothing
             jac_alg::J = BVPJacobianAlgorithm()
             nested_nlsolve::Bool = false
             nested_nlsolve_kwargs::NamedTuple = (;)
@@ -308,13 +327,16 @@ for stage in (2, 3, 4, 5)
             max_num_subintervals::Int = 3000
         end
         $(alg)(nlsolve::N,
+            optimize::O,
             jac_alg::J;
             nested = false,
             nested_nlsolve_kwargs::NamedTuple = (;),
             defect_threshold::T = 0.1,
             max_num_subintervals::Int = 3000) where {N,
+            O,
             J,
-            T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
+            T} = $(alg){N, O, J, T}(
+            nlsolve, optimize, jac_alg, nested, nested_nlsolve_kwargs,
             defect_threshold, max_num_subintervals)
     end
 end
@@ -335,6 +357,9 @@ for stage in (2, 3, 4, 5)
             `NonlinearProblem` interface can be used. Note that any autodiff argument for
             the solver will be ignored and a custom jacobian algorithm will be used.
 
+          - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML
+            `OptimizationProblem` interface can be used. Note that any autodiff argument for
+            the solver will be ignored and a custom jacobian algorithm will be used.
           - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to
             `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to
             use based on the input types and problem type.
@@ -406,8 +431,9 @@ for stage in (2, 3, 4, 5)
         }
         ```
         """
-        Base.@kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK
+        Base.@kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractFIRK
             nlsolve::N = nothing
+            optimize::O = nothing
             jac_alg::J = BVPJacobianAlgorithm()
             nested_nlsolve::Bool = false
             nested_nlsolve_kwargs::NamedTuple = (;)
@@ -415,13 +441,16 @@ for stage in (2, 3, 4, 5)
             max_num_subintervals::Int = 3000
         end
         $(alg)(nlsolve::N,
+            optimize::O,
             jac_alg::J;
             nested = false,
             nested_nlsolve_kwargs::NamedTuple = (;),
             defect_threshold::T = 0.1,
             max_num_subintervals::Int = 3000) where {N,
+            O,
             J,
-            T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
+            T} = $(alg){N, O, J, T}(
+            nlsolve, optimize, jac_alg, nested, nested_nlsolve_kwargs,
             defect_threshold, max_num_subintervals)
     end
 end
diff --git a/lib/BoundaryValueDiffEqFIRK/src/firk.jl b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
index 9df4ce4bd..13862b693 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/firk.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
@@ -27,6 +27,7 @@
     nest_prob
     resid_size
     nlsolve_kwargs
+    optimize_kwargs
     kwargs
 end
 
@@ -60,6 +61,7 @@ Base.eltype(::FIRKCacheNested{iip, T}) where {iip, T} = T
     defect
     resid_size
     nlsolve_kwargs
+    optimize_kwargs
     kwargs
 end
 
@@ -85,19 +87,23 @@ end
 
 function SciMLBase.__init(
         prob::BVProblem, alg::AbstractFIRK; dt = 0.0, abstol = 1e-6, adaptive = true,
-        controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), kwargs...)
+        controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol),
+        optimize_kwargs = (; abstol = abstol), kwargs...)
     if alg.nested_nlsolve
         return init_nested(prob, alg; dt = dt, abstol = abstol, adaptive = adaptive,
-            controller = controller, nlsolve_kwargs = nlsolve_kwargs, kwargs...)
+            controller = controller, nlsolve_kwargs = nlsolve_kwargs,
+            optimize_kwargs = optimize_kwargs, kwargs...)
     else
         return init_expanded(prob, alg; dt = dt, abstol = abstol, adaptive = adaptive,
-            controller = controller, nlsolve_kwargs = nlsolve_kwargs, kwargs...)
+            controller = controller, nlsolve_kwargs = nlsolve_kwargs,
+            optimize_kwargs = optimize_kwargs, kwargs...)
     end
 end
 
 function init_nested(
         prob::BVProblem, alg::AbstractFIRK; dt = 0.0, abstol = 1e-6, adaptive = true,
-        controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), kwargs...)
+        controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol),
+        optimize_kwargs = (; abstol = abstol), kwargs...)
     @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
 
     iip = isinplace(prob)
@@ -200,15 +206,16 @@ function init_nested(
     end
 
     return FIRKCacheNested{iip, T, typeof(diffcache), fit_parameters}(
-        alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type,
-        prob.p, alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete,
-        y, y₀, residual, fᵢ_cache, fᵢ₂_cache, defect, nestprob, resid₁_size,
-        nlsolve_kwargs, (; abstol, dt, adaptive, controller, kwargs...))
+        alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type, prob.p,
+        alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete, y, y₀, residual,
+        fᵢ_cache, fᵢ₂_cache, defect, nestprob, resid₁_size, nlsolve_kwargs,
+        optimize_kwargs, (; abstol, dt, adaptive, controller, kwargs...))
 end
 
 function init_expanded(
         prob::BVProblem, alg::AbstractFIRK; dt = 0.0, abstol = 1e-6, adaptive = true,
-        controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol), kwargs...)
+        controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol),
+        optimize_kwargs = (; abstol = abstol), kwargs...)
     @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
 
     if adaptive && isa(alg, FIRKNoAdaptivity)
@@ -301,9 +308,9 @@ function init_expanded(
 
     return FIRKCacheExpand{iip, T, typeof(diffcache), fit_parameters}(
         alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type,
-        prob.p, alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete,
-        y, y₀, residual, fᵢ_cache, fᵢ₂_cache, defect, resid₁_size,
-        nlsolve_kwargs, (; abstol, dt, adaptive, controller, kwargs...))
+        prob.p, alg, TU, ITU, bcresid_prototype, mesh, mesh_dt, k_discrete, y,
+        y₀, residual, fᵢ_cache, fᵢ₂_cache, defect, resid₁_size, nlsolve_kwargs,
+        optimize_kwargs, (; abstol, dt, adaptive, controller, kwargs...))
 end
 
 """
@@ -403,9 +410,11 @@ end
 
 function __perform_firk_iteration(
         cache::Union{FIRKCacheExpand, FIRKCacheNested}, abstol, adaptive::Bool)
-    nlprob = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀))
-    nlsolve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve)
-    sol_nlprob = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true)
+    nlprob = __construct_problem(cache, vec(cache.y₀), copy(cache.y₀))
+    solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize)
+    kwargs = __concrete_kwargs(
+        cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs)
+    sol_nlprob = solve(nlprob, solve_alg, kwargs...)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
     defect_norm = 2 * abstol
@@ -450,7 +459,7 @@ function __perform_firk_iteration(
 end
 
 # Constructing the Nonlinear Problem
-function __construct_nlproblem(cache::Union{FIRKCacheNested{iip}, FIRKCacheExpand{iip}},
+function __construct_problem(cache::Union{FIRKCacheNested{iip}, FIRKCacheExpand{iip}},
         y::AbstractVector, y₀::AbstractVectorOfArray) where {iip}
     pt = cache.problem_type
     (; jac_alg) = cache.alg
@@ -490,10 +499,10 @@ function __construct_nlproblem(cache::Union{FIRKCacheNested{iip}, FIRKCacheExpan
             u, p, cache.y, pt, cache.bc, cache.mesh, cache, eval_sol, trait)
     end
 
-    return __construct_nlproblem(cache, y, loss_bc, loss_collocation, loss, pt)
+    return __construct_problem(cache, y, loss_bc, loss_collocation, loss, pt)
 end
 
-function __construct_nlproblem(
+function __construct_problem(
         cache::FIRKCacheExpand{iip}, y, loss_bc::BC, loss_collocation::C,
         loss::LF, ::StandardBVProblem) where {iip, BC, C, LF}
     (; alg, stage) = cache
@@ -574,13 +583,11 @@ function __construct_nlproblem(
     end
 
     resid_prototype = vcat(resid_bc, resid_collocation)
-    nlf = NonlinearFunction{iip}(
-        loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
-
-    return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p)
+    return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype,
+        resid_prototype, y, cache.p, cache.M, N)
 end
 
-function __construct_nlproblem(
+function __construct_problem(
         cache::FIRKCacheExpand{iip}, y, loss_bc::BC, loss_collocation::C,
         loss::LF, ::TwoPointBVProblem) where {iip, BC, C, LF}
     (; jac_alg) = cache.alg
@@ -634,12 +641,11 @@ function __construct_nlproblem(
     end
 
     resid_prototype = copy(resid)
-    nlf = NonlinearFunction{iip}(
-        loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
-    return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p)
+    return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype,
+        resid_prototype, y, cache.p, cache.M, N)
 end
 
-function __construct_nlproblem(
+function __construct_problem(
         cache::FIRKCacheNested{iip}, y, loss_bc::BC, loss_collocation::C,
         loss::LF, ::StandardBVProblem) where {iip, BC, C, LF}
     (; jac_alg) = cache.alg
@@ -714,13 +720,11 @@ function __construct_nlproblem(
     end
 
     resid_prototype = vcat(resid_bc, resid_collocation)
-    nlf = NonlinearFunction{iip}(
-        loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
-
-    return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p)
+    return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype,
+        resid_prototype, y, cache.p, cache.M, N)
 end
 
-function __construct_nlproblem(
+function __construct_problem(
         cache::FIRKCacheNested{iip}, y, loss_bc::BC, loss_collocation::C,
         loss::LF, ::TwoPointBVProblem) where {iip, BC, C, LF}
     (; jac_alg) = cache.alg
@@ -765,9 +769,8 @@ function __construct_nlproblem(
     end
 
     resid_prototype = copy(resid)
-    nlf = NonlinearFunction{iip}(
-        loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
-    return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p)
+    return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype,
+        resid_prototype, y, cache.p, cache.M, N)
 end
 
 @views function __firk_loss!(resid, u, p, y, pt::StandardBVProblem, bc!::BC, residual,
@@ -829,7 +832,6 @@ end
 @views function __firk_loss(u, p, y::AbstractVector, pt::TwoPointBVProblem,
         bc::Tuple{BC1, BC2}, mesh, cache, _, trait) where {BC1, BC2}
     y_ = recursive_unflatten!(y, u)
-    soly_ = VectorOfArray(y_)
     resid_bca, resid_bcb = eval_bc_residual(pt, bc, y_, p, mesh)
     resid_co = Φ(cache, y_, u, trait)
     return vcat(resid_bca, mapreduce(vec, vcat, resid_co), resid_bcb)
diff --git a/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl b/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl
index 44a3a11a9..eae0b5474 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/algorithms.jl
@@ -19,7 +19,8 @@ for order in (2, 3, 4, 5, 6)
 
           - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML
             `OptimizationProblem` interface can be used. Note that any autodiff argument for
-            the solver will be ignored and a custom jacobian algorithm will be used.
+            the solver will be ignored and a custom jacobian algorithm will be used. Optimization
+            solvers should first be loaded to allow this functionality.
           - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to
             `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to
             use based on the input types and problem type.
diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
index fd714b087..0a646edb9 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
@@ -41,6 +41,7 @@ function SciMLBase.__init(
     @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
     iip = isinplace(prob)
     diffcache = __cache_trait(alg.jac_alg)
+    @assert (!iip & !isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
     fit_parameters = haskey(prob.kwargs, :fit_parameters)
 
     t₀, t₁ = prob.tspan
@@ -490,6 +491,7 @@ end
 function __mirk_mpoint_jacobian(
         J, _, x, bc_diffmode, nonbc_diffmode, bc_diffcache, nonbc_diffcache,
         loss_bc::BC, loss_collocation::C, L::Int, p) where {BC, C}
+    println("p: ", p)
     DI.jacobian!(loss_bc, @view(J[1:L, :]), bc_diffcache, bc_diffmode, x, Constant(p))
     DI.jacobian!(loss_collocation, @view(J[(L + 1):end, :]),
         nonbc_diffcache, nonbc_diffmode, x, Constant(p))

From 93ea85a796f6dad9257eb887e53ecbc878890abd Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Fri, 18 Jul 2025 14:03:41 +0800
Subject: [PATCH 08/50] More docstrings

---
 docs/src/basics/solve.md                      |  3 +-
 lib/BoundaryValueDiffEqFIRK/src/algorithms.jl | 28 ++++++++-----------
 lib/BoundaryValueDiffEqFIRK/src/firk.jl       |  6 ++--
 3 files changed, 17 insertions(+), 20 deletions(-)

diff --git a/docs/src/basics/solve.md b/docs/src/basics/solve.md
index 20c1ada4f..281fc94b8 100644
--- a/docs/src/basics/solve.md
+++ b/docs/src/basics/solve.md
@@ -5,6 +5,7 @@
   - `controller`: Error controller for collocation methods, default as `DefectControl()`, more controller options in [Error Control Adaptivity](@ref error_control).
   - `defect_threshold`: Monitor of the size of defect norm. Defaults to `0.1`.
   - `odesolve_kwargs`: OrdinaryDiffEq.jl solvers kwargs for passing to ODE solving in shooting methods. For more information, see the documentation for OrdinaryDiffEq: [Common Solver Options](https://docs.sciml.ai/DiffEqDocs/latest/basics/common_solver_opts/).
-  - `nlsolve_kwargs`: NonlinearSolve.jl solvers kwargs for passing to nonlinear solving in collocation methods and shooting methods. For more information, see the documentation for NonlinearSolve: [Commom Solver Options](https://docs.sciml.ai/NonlinearSolve/stable/basics/solve/). The default absolute tolerance of nonlinear solving in collocaio
+  - `nlsolve_kwargs`: NonlinearSolve.jl solvers kwargs for passing to nonlinear solving in collocation methods and shooting methods. For more information, see the documentation for NonlinearSolve: [Commom Solver Options](https://docs.sciml.ai/NonlinearSolve/stable/basics/solve/). The default internal nonlinear solver is [customized polyalgorithm](https://github.com/SciML/BoundaryValueDiffEq.jl/blob/master/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl) and the default absolute tolerance of nonlinear solving in collocation methods is `1e-6`.
+  - `optimize_kwargs`: Optimization.jl solvers kwargs for passing to optimizing in collocation methods and shooting methods. For more information, see the documentation for Optimization: [Commom Solver Options](https://docs.sciml.ai/Optimization/stable/API/solve/). The internal optimization solver should be specified and the default absolute tolerance of optimization problem solving in collocation methods is `1e-6`.
   - `verbose`:  Toggles whether warnings are thrown when the solver exits early. Defaults to `true`.
   - `ensemblealg`: Whether `MultipleShooting` uses multithreading, default as `EnsembleThreads()`. For more information, see the documentation for OrdinaryDiffEq: [EnsembleAlgorithms](https://docs.sciml.ai/DiffEqDocs/latest/features/ensemble/#EnsembleAlgorithms).
diff --git a/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl b/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl
index 1e093a45f..ec957fa4d 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl
@@ -102,15 +102,14 @@ for stage in (1, 2, 3, 5, 7)
         $(alg)(nlsolve::N,
             optimize::O,
             jac_alg::J;
-            nested = false,
+            nested_nlsolve = false,
             nested_nlsolve_kwargs::NamedTuple = (;),
             defect_threshold::T = 0.1,
             max_num_subintervals::Int = 3000) where {N,
             O,
             J,
-            T} = $(alg){N, O, J, T}(
-            nlsolve, optimize, jac_alg, nested, nested_nlsolve_kwargs,
-            defect_threshold, max_num_subintervals)
+            T} = $(alg){N, O, J, T}(nlsolve, optimize, jac_alg, nested_nlsolve,
+            nested_nlsolve_kwargs, defect_threshold, max_num_subintervals)
     end
 end
 
@@ -215,15 +214,14 @@ for stage in (2, 3, 4, 5)
         $(alg)(nlsolve::N,
             optimize::O,
             jac_alg::J;
-            nested = false,
+            nested_nlsolve = false,
             nested_nlsolve_kwargs::NamedTuple = (;),
             defect_threshold::T = 0.1,
             max_num_subintervals::Int = 3000) where {N,
             O,
             J,
-            T} = $(alg){N, O, J, T}(
-            nlsolve, optimize, jac_alg, nested, nested_nlsolve_kwargs,
-            defect_threshold, max_num_subintervals)
+            T} = $(alg){N, O, J, T}(nlsolve, optimize, jac_alg, nested_nlsolve,
+            nested_nlsolve_kwargs, defect_threshold, max_num_subintervals)
     end
 end
 
@@ -329,15 +327,14 @@ for stage in (2, 3, 4, 5)
         $(alg)(nlsolve::N,
             optimize::O,
             jac_alg::J;
-            nested = false,
+            nested_nlsolve = false,
             nested_nlsolve_kwargs::NamedTuple = (;),
             defect_threshold::T = 0.1,
             max_num_subintervals::Int = 3000) where {N,
             O,
             J,
-            T} = $(alg){N, O, J, T}(
-            nlsolve, optimize, jac_alg, nested, nested_nlsolve_kwargs,
-            defect_threshold, max_num_subintervals)
+            T} = $(alg){N, O, J, T}(nlsolve, optimize, jac_alg, nested_nlsolve,
+            nested_nlsolve_kwargs, defect_threshold, max_num_subintervals)
     end
 end
 
@@ -443,15 +440,14 @@ for stage in (2, 3, 4, 5)
         $(alg)(nlsolve::N,
             optimize::O,
             jac_alg::J;
-            nested = false,
+            nested_nlsolve = false,
             nested_nlsolve_kwargs::NamedTuple = (;),
             defect_threshold::T = 0.1,
             max_num_subintervals::Int = 3000) where {N,
             O,
             J,
-            T} = $(alg){N, O, J, T}(
-            nlsolve, optimize, jac_alg, nested, nested_nlsolve_kwargs,
-            defect_threshold, max_num_subintervals)
+            T} = $(alg){N, O, J, T}(nlsolve, optimize, jac_alg, nested_nlsolve,
+            nested_nlsolve_kwargs, defect_threshold, max_num_subintervals)
     end
 end
 
diff --git a/lib/BoundaryValueDiffEqFIRK/src/firk.jl b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
index 13862b693..408769c0d 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/firk.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
@@ -107,6 +107,7 @@ function init_nested(
     @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
 
     iip = isinplace(prob)
+    @assert (!iip & !isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
     if adaptive && isa(alg, FIRKNoAdaptivity)
         error("Algorithm doesn't support adaptivity. Please choose a higher order algorithm.")
     end
@@ -217,15 +218,14 @@ function init_expanded(
         controller = DefectControl(), nlsolve_kwargs = (; abstol = abstol),
         optimize_kwargs = (; abstol = abstol), kwargs...)
     @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
-
+    iip = isinplace(prob)
+    @assert (!iip & !isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
     if adaptive && isa(alg, FIRKNoAdaptivity)
         error("Algorithm $(alg) doesn't support adaptivity. Please choose a higher order algorithm.")
     end
     diffcache = __cache_trait(alg.jac_alg)
     fit_parameters = haskey(prob.kwargs, :fit_parameters)
 
-    iip = isinplace(prob)
-
     t₀, t₁ = prob.tspan
     ig, T,
     M,

From 54943e8bc026ee8495a9f2ffb7f43af06ce53d13 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Fri, 18 Jul 2025 14:40:23 +0800
Subject: [PATCH 09/50] Proper assert

---
 lib/BoundaryValueDiffEqFIRK/src/firk.jl | 4 ++--
 lib/BoundaryValueDiffEqMIRK/src/mirk.jl | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/src/firk.jl b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
index 408769c0d..cf53c0019 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/firk.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
@@ -107,7 +107,7 @@ function init_nested(
     @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
 
     iip = isinplace(prob)
-    @assert (!iip & !isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
+    @assert (iip || isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
     if adaptive && isa(alg, FIRKNoAdaptivity)
         error("Algorithm doesn't support adaptivity. Please choose a higher order algorithm.")
     end
@@ -219,7 +219,7 @@ function init_expanded(
         optimize_kwargs = (; abstol = abstol), kwargs...)
     @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
     iip = isinplace(prob)
-    @assert (!iip & !isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
+    @assert (iip || isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
     if adaptive && isa(alg, FIRKNoAdaptivity)
         error("Algorithm $(alg) doesn't support adaptivity. Please choose a higher order algorithm.")
     end
diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
index 0a646edb9..a77644734 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
@@ -41,7 +41,7 @@ function SciMLBase.__init(
     @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
     iip = isinplace(prob)
     diffcache = __cache_trait(alg.jac_alg)
-    @assert (!iip & !isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
+    @assert (iip || isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
     fit_parameters = haskey(prob.kwargs, :fit_parameters)
 
     t₀, t₁ = prob.tspan

From 5c2146295952463eedaae8835067d44b360c22d9 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Fri, 18 Jul 2025 16:07:08 +0800
Subject: [PATCH 10/50] Expanded FIRK should work

---
 lib/BoundaryValueDiffEqCore/src/utils.jl             | 6 +++---
 lib/BoundaryValueDiffEqFIRK/src/firk.jl              | 4 ++--
 lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl | 3 ++-
 3 files changed, 7 insertions(+), 6 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 76de3867f..8dffc9d02 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -632,7 +632,7 @@ function __construct_internal_problem(
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
-        optf = OptimizationFunction{iip}(__default_cost(prob.f), AutoFiniteDiff(), # Need to investigate the ForwardDiff dual problem
+        optf = OptimizationFunction{true}(__default_cost(prob.f), AutoFiniteDiff(), # Need to investigate the ForwardDiff dual problem
             cons = loss,
             cons_j = jac, cons_jac_prototype = jac_prototype)
         lcons = zeros(T, N*M)
@@ -651,9 +651,9 @@ function __construct_internal_problem(
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
-        optf = OptimizationFunction{iip}(
+        optf = OptimizationFunction{true}(
             __default_cost(prob.f), get_dense_ad(diffmode), cons = loss,
-            cons_j = jac, cons_jac_prototype = Matrix(jac_prototype))
+            cons_j = jac, cons_jac_prototype = jac_prototype)
         lcons = zeros(T, N*M)
         ucons = zeros(T, N*M)
 
diff --git a/lib/BoundaryValueDiffEqFIRK/src/firk.jl b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
index cf53c0019..2f513e034 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/firk.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
@@ -584,7 +584,7 @@ function __construct_problem(
 
     resid_prototype = vcat(resid_bc, resid_collocation)
     return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype,
-        resid_prototype, y, cache.p, cache.M, N)
+        resid_prototype, y, cache.p, cache.M, (N - 1) * (stage + 1) + 1)
 end
 
 function __construct_problem(
@@ -642,7 +642,7 @@ function __construct_problem(
 
     resid_prototype = copy(resid)
     return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype,
-        resid_prototype, y, cache.p, cache.M, N)
+        resid_prototype, y, cache.p, cache.M, (N - 1) * (stage + 1) + 1)
 end
 
 function __construct_problem(
diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
index 75fcb840e..f141c8ab5 100644
--- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
@@ -466,7 +466,8 @@ end
 @testitem "Convergence with optimization based solver" setup=[MIRKConvergenceTests] begin
     using LinearAlgebra, DiffEqDevTools, OptimizationMOI, Ipopt
 
-    @testset "Problem: $i" for i in (3, 4, 5, 6, 9, 10)
+    # Only test on inplace problems
+    @testset "Problem: $i" for i in (3, 5, 9)
         prob = probArr[i]
         @testset "MIRK$order" for (_, order) in enumerate((2, 3, 4, 5, 6))
             sim = test_convergence(

From ce2b3e0d6546004c4634fab15be6766bb2af793b Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 19 Jul 2025 00:30:16 +0800
Subject: [PATCH 11/50] Fix Aqua error

---
 lib/BoundaryValueDiffEqCore/src/utils.jl      | 6 ++++++
 lib/BoundaryValueDiffEqMIRK/Project.toml      | 2 ++
 lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl | 5 ++---
 3 files changed, 10 insertions(+), 3 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 8dffc9d02..4437dfce6 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -623,6 +623,12 @@ end
 @inline __default_cost(f) = f
 @inline __default_cost(fun::BVPFunction) = __default_cost(fun.cost)
 
+"""
+    __construct_internal_problem
+
+Constructs the internal problem based on the type of the boundary value problem and the
+algorithm used. It returns either a `NonlinearProblem` or an `OptimizationProblem`.
+"""
 function __construct_internal_problem(
         prob::BVProblem, alg, loss, jac, jac_prototype, resid_prototype, y, p, M, N)
     T = eltype(y)
diff --git a/lib/BoundaryValueDiffEqMIRK/Project.toml b/lib/BoundaryValueDiffEqMIRK/Project.toml
index 818afe79b..28e31c927 100644
--- a/lib/BoundaryValueDiffEqMIRK/Project.toml
+++ b/lib/BoundaryValueDiffEqMIRK/Project.toml
@@ -42,10 +42,12 @@ FastClosures = "0.3.2"
 ForwardDiff = "0.10.38, 1"
 Hwloc = "3"
 InteractiveUtils = "<0.0.1, 1"
+Ipopt = "1.10.6"
 JET = "0.9.18"
 LinearAlgebra = "1.10"
 LinearSolve = "2.36.2, 3"
 Mooncake = "0.4.108"
+OptimizationMOI = ""
 OrdinaryDiffEqRosenbrock = "1"
 PreallocationTools = "0.4.24"
 PrecompileTools = "1.2"
diff --git a/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl b/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl
index eff1d50ba..f77673a8b 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/adaptivity.jl
@@ -498,8 +498,7 @@ end
     new_prob = remake(prob, u0 = high_sol)
     high_cache = SciMLBase.__init(new_prob, alg, adaptive = false)
 
-    high_nlprob = __construct_nlproblem(
-        high_cache, vec(high_sol), VectorOfArray(high_sol.u))
+    high_nlprob = __construct_problem(high_cache, vec(high_sol), VectorOfArray(high_sol.u))
     high_sol_original = __solve(
         high_nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true)
     recursive_unflatten!(high_sol, high_sol_original.u)
@@ -518,7 +517,7 @@ end
     new_prob = remake(prob, u0 = high_sol)
     high_cache = SciMLBase.__init(new_prob, __high_order_method(alg), adaptive = false)
 
-    high_nlprob = __construct_nlproblem(high_cache, sol.u, high_sol)
+    high_nlprob = __construct_problem(high_cache, sol.u, high_sol)
     high_sol_nlprob = __solve(
         high_nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true)
     recursive_unflatten!(high_sol, high_sol_nlprob)

From 3b905c9ab6a4aa95143ceaf50ba063bb39a5758d Mon Sep 17 00:00:00 2001
From: Qingyu Qu <52615090+ErikQQY@users.noreply.github.com>
Date: Sat, 19 Jul 2025 04:40:28 +0800
Subject: [PATCH 12/50] Fix typo in compat

---
 lib/BoundaryValueDiffEqMIRK/Project.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/BoundaryValueDiffEqMIRK/Project.toml b/lib/BoundaryValueDiffEqMIRK/Project.toml
index 28e31c927..c5a8938a5 100644
--- a/lib/BoundaryValueDiffEqMIRK/Project.toml
+++ b/lib/BoundaryValueDiffEqMIRK/Project.toml
@@ -47,7 +47,7 @@ JET = "0.9.18"
 LinearAlgebra = "1.10"
 LinearSolve = "2.36.2, 3"
 Mooncake = "0.4.108"
-OptimizationMOI = ""
+OptimizationMOI = "0.5.5"
 OrdinaryDiffEqRosenbrock = "1"
 PreallocationTools = "0.4.24"
 PrecompileTools = "1.2"

From a041b901f267e157bd72cd4d57120b3e86c46ff3 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Mon, 21 Jul 2025 02:11:32 +0800
Subject: [PATCH 13/50] Shooting and MultipleShooting work fine

---
 lib/BoundaryValueDiffEqMIRK/src/mirk.jl       |  1 -
 .../src/BoundaryValueDiffEqShooting.jl        |  3 ++-
 .../src/algorithms.jl                         | 27 ++++++++++++-------
 .../src/multiple_shooting.jl                  | 27 ++++++++++---------
 .../src/single_shooting.jl                    | 14 +++++-----
 5 files changed, 42 insertions(+), 30 deletions(-)

diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
index a77644734..a1aed580c 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
@@ -491,7 +491,6 @@ end
 function __mirk_mpoint_jacobian(
         J, _, x, bc_diffmode, nonbc_diffmode, bc_diffcache, nonbc_diffcache,
         loss_bc::BC, loss_collocation::C, L::Int, p) where {BC, C}
-    println("p: ", p)
     DI.jacobian!(loss_bc, @view(J[1:L, :]), bc_diffcache, bc_diffmode, x, Constant(p))
     DI.jacobian!(loss_collocation, @view(J[(L + 1):end, :]),
         nonbc_diffcache, nonbc_diffmode, x, Constant(p))
diff --git a/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl b/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
index 7b0169f20..20ad0b2a5 100644
--- a/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
@@ -8,7 +8,8 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm, BVPJacobian
                                __concrete_solve_algorithm, diff!, __any_sparse_ad,
                                __cache_trait, concrete_jacobian_algorithm, eval_bc_residual,
                                eval_bc_residual!, get_tmp, __maybe_matmul!,
-                               __extract_problem_details, __initial_guess,
+                               __concrete_kwargs, __extract_problem_details,
+                               __initial_guess, __construct_internal_problem,
                                __default_coloring_algorithm, __default_sparsity_detector,
                                __maybe_allocate_diffcache, __get_bcresid_prototype,
                                safe_similar, __vec, __vec_f, __vec_f!, __vec_bc, __vec_bc!,
diff --git a/lib/BoundaryValueDiffEqShooting/src/algorithms.jl b/lib/BoundaryValueDiffEqShooting/src/algorithms.jl
index afcaa9d6d..1b9484b26 100644
--- a/lib/BoundaryValueDiffEqShooting/src/algorithms.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/algorithms.jl
@@ -4,7 +4,7 @@ abstract type AbstractShooting <: AbstractBoundaryValueDiffEqAlgorithm end
 """
     Shooting(ode_alg; kwargs...)
     Shooting(ode_alg, nlsolve; kwargs...)
-    Shooting(; ode_alg = nothing, nlsolve = nothing, jac_alg = nothing)
+    Shooting(; ode_alg = nothing, nlsolve = nothing, optimize = nothing, jac_alg = nothing)
 
 Single shooting method, reduces BVP to an initial value problem and solves the IVP.
 
@@ -16,6 +16,9 @@ Single shooting method, reduces BVP to an initial value problem and solves the I
   - `nlsolve`: Internal Nonlinear solver. Any solver which conforms to the SciML
     `NonlinearProblem` interface can be used. Note that any autodiff argument for the solver
     will be ignored and a custom jacobian algorithm will be used.
+  - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML
+    `OptimizationProblem` interface can be used. Note that any autodiff argument for the solver
+    will be ignored and a custom jacobian algorithm will be used.
   - `jac_alg`: Jacobian Algorithm used for the Nonlinear Solver. If this is not set, we
     check if `nlsolve.ad` exists and is not nothing. If it is, we use that to construct
     the jacobian. If not, we try to use the best algorithm based on the input types
@@ -25,11 +28,14 @@ Single shooting method, reduces BVP to an initial value problem and solves the I
 @concrete struct Shooting{J <: BVPJacobianAlgorithm} <: AbstractShooting
     ode_alg
     nlsolve
+    optimize
     jac_alg::J
 end
 
-function Shooting(; ode_alg = nothing, nlsolve = nothing, jac_alg = nothing)
-    return Shooting(ode_alg, nlsolve, __materialize_jacobian_algorithm(nlsolve, jac_alg))
+function Shooting(;
+        ode_alg = nothing, nlsolve = nothing, optimize = nothing, jac_alg = nothing)
+    return Shooting(
+        ode_alg, nlsolve, optimize, __materialize_jacobian_algorithm(nlsolve, jac_alg))
 end
 @inline Shooting(ode_alg; kwargs...) = Shooting(; ode_alg, kwargs...)
 @inline Shooting(ode_alg, nlsolve; kwargs...) = Shooting(; ode_alg, nlsolve, kwargs...)
@@ -42,7 +48,7 @@ end
 
 """
     MultipleShooting(; nshoots::Int, ode_alg = nothing, nlsolve = nothing,
-        grid_coarsening = true, jac_alg = nothing)
+        optimize = nothing, grid_coarsening = true, jac_alg = nothing)
     MultipleShooting(nshoots::Int; kwargs...)
     MultipleShooting(nshoots::Int, ode_alg; kwargs...)
     MultipleShooting(nshoots::Int, ode_alg, nlsolve; kwargs...)
@@ -59,6 +65,8 @@ Significantly more stable than Single Shooting.
     poly-algorithm if `DifferentialEquations.jl` is loaded else this must be supplied)
   - `nlsolve`: Internal Nonlinear solver. Any solver which conforms to the SciML
     `NonlinearProblem` interface can be used.
+  - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML
+    `OptimizationProblem` interface can be used.
   - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to
     `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to use based
     on the input types and problem type.
@@ -85,6 +93,7 @@ Significantly more stable than Single Shooting.
 @concrete struct MultipleShooting{J <: BVPJacobianAlgorithm} <: AbstractShooting
     ode_alg
     nlsolve
+    optimize
     jac_alg::J
     nshoots::Int
     grid_coarsening
@@ -93,17 +102,18 @@ end
 function concretize_jacobian_algorithm(alg::MultipleShooting, prob)
     jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
     return MultipleShooting(
-        alg.ode_alg, alg.nlsolve, jac_alg, alg.nshoots, alg.grid_coarsening)
+        alg.ode_alg, alg.nlsolve, alg.optimize, jac_alg, alg.nshoots, alg.grid_coarsening)
 end
 
 function update_nshoots(alg::MultipleShooting, nshoots::Int)
     return MultipleShooting(
-        alg.ode_alg, alg.nlsolve, alg.jac_alg, nshoots, alg.grid_coarsening)
+        alg.ode_alg, alg.nlsolve, alg.optimize, alg.jac_alg, nshoots, alg.grid_coarsening)
 end
 
 function MultipleShooting(; nshoots::Int,
         ode_alg = nothing,
         nlsolve = nothing,
+        optimize = nothing,
         grid_coarsening::Union{
             Bool, Function, <:AbstractVector{<:Integer}, Tuple{Vararg{Integer}}} = true,
         jac_alg = nothing)
@@ -112,9 +122,8 @@ function MultipleShooting(; nshoots::Int,
         sort!(grid_coarsening; rev = true)
         @assert all(grid_coarsening .> 0) && 1 ∉ grid_coarsening
     end
-    return MultipleShooting(
-        ode_alg, nlsolve, __materialize_jacobian_algorithm(nlsolve, jac_alg),
-        nshoots, grid_coarsening)
+    return MultipleShooting(ode_alg, nlsolve, optimize,
+        __materialize_jacobian_algorithm(nlsolve, jac_alg), nshoots, grid_coarsening)
 end
 @inline MultipleShooting(nshoots::Int; kwargs...) = MultipleShooting(; nshoots, kwargs...)
 @inline MultipleShooting(
diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
index 0f5be76d3..aa43fb877 100644
--- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
@@ -1,5 +1,6 @@
 function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwargs = (;),
-        nlsolve_kwargs = (;), ensemblealg = EnsembleThreads(), verbose = true, kwargs...)
+        nlsolve_kwargs = (;), optimize_kwargs = (;),
+        ensemblealg = EnsembleThreads(), verbose = true, kwargs...)
     (; f, tspan) = prob
 
     if !(ensemblealg isa EnsembleSerial) && !(ensemblealg isa EnsembleThreads)
@@ -81,11 +82,11 @@ function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwa
         diffmode_shooting = __get_non_sparse_ad(alg.jac_alg.bc_diffmode)
     end
     shooting_alg = Shooting(
-        alg.ode_alg, alg.nlsolve, BVPJacobianAlgorithm(diffmode_shooting))
+        alg.ode_alg, alg.nlsolve, alg.optimize, BVPJacobianAlgorithm(diffmode_shooting))
 
     single_shooting_prob = remake(prob; u0 = reshape(u_at_nodes[1:N], u0_size))
     return __solve(single_shooting_prob, shooting_alg; odesolve_kwargs,
-        nlsolve_kwargs, verbose, kwargs...)
+        nlsolve_kwargs, optimize_kwargs, verbose, kwargs...)
 end
 
 # TODO: We can save even more memory by hoisting the preallocated caches for the ODEs
@@ -141,9 +142,12 @@ function __solve_nlproblem!(
         jac_prototype = jac_prototype)
 
     # NOTE: u_at_nodes is updated inplace
-    nlprob = __internal_nlsolve_problem(prob, M, N, loss_function!, u_at_nodes, prob.p)
-    nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve)
-    __solve(nlprob, nlsolve_alg; kwargs..., alias_u0 = true)
+    nlprob = __construct_internal_problem(prob, alg, loss_fn, jac_fn, jac_prototype,
+        resid_prototype, u_at_nodes, prob.p, M, length(nodes))
+
+    nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize)
+    kwargs = __concrete_kwargs(alg.nlsolve, alg.optimize, kwargs...)
+    solve(nlprob, nlsolve_alg; kwargs...)
 
     return nothing
 end
@@ -213,13 +217,12 @@ function __solve_nlproblem!(::StandardBVProblem, alg::MultipleShooting, bcresid_
         J, u, p, similar(bcresid_prototype), resid_nodes, ode_jac_cache, bc_jac_cache,
         ode_fn, bc_fn, nonbc_diffmode, bc_diffmode, N, M, __cache_trait(alg.jac_alg))
 
-    loss_function! = NonlinearFunction{true}(loss_fn; resid_prototype = resid_prototype,
-        jac_prototype = jac_prototype, jac = jac_fn)
-
     # NOTE: u_at_nodes is updated inplace
-    nlprob = __internal_nlsolve_problem(prob, M, N, loss_function!, u_at_nodes, prob.p)
-    nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve)
-    __solve(nlprob, nlsolve_alg; kwargs..., alias_u0 = true)
+    nlprob = __construct_internal_problem(prob, alg, loss_fn, jac_fn, jac_prototype,
+        resid_prototype, u_at_nodes, prob.p, M, length(nodes))
+    nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize)
+    kwargs = __concrete_kwargs(alg.nlsolve, alg.optimize, kwargs...)
+    solve(nlprob, nlsolve_alg; kwargs...)
 
     return nothing
 end
diff --git a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
index 1a03d028f..13bb886c6 100644
--- a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
@@ -1,5 +1,5 @@
 function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (;),
-        nlsolve_kwargs = (;), verbose = true, kwargs...)
+        nlsolve_kwargs = (;), optimize_kwargs = (;), verbose = true, kwargs...)
     # Setup the problem
     if prob.u0 isa AbstractArray{<:Number}
         u0 = prob.u0
@@ -73,12 +73,12 @@ function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (;
             jac_prototype, u, jac_cache, diffmode, loss_fnₚ)
     end
 
-    nlf = NonlinearFunction{iip}(loss_fn; jac_prototype = jac_prototype,
-        resid_prototype = resid_prototype, jac = jac_fn)
-    nlprob = __internal_nlsolve_problem(prob, resid_prototype, u0, nlf, vec(u0), prob.p)
-    nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve)
-    nlsol::SciMLBase.NonlinearSolution = __solve(
-        nlprob, nlsolve_alg; nlsolve_kwargs..., verbose, kwargs...)
+    nlprob = __construct_internal_problem(prob, alg, loss_fn, jac_fn, jac_prototype,
+        resid_prototype, u0, prob.p, length(u0), 1)
+    solve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize)
+    kwargs = __concrete_kwargs(alg.nlsolve, alg.optimize, nlsolve_kwargs, optimize_kwargs)
+    #TODO: add verbose kwarg
+    nlsol = solve(nlprob, solve_alg; kwargs...)
 
     # There is no way to reinit with the same cache with different cache. But not saving
     # the internal values gives a significant speedup. So we just create a new cache

From b7f5762c6a19f3e356c697749392a407b6c34357 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Mon, 21 Jul 2025 16:17:28 +0800
Subject: [PATCH 14/50] MIRKN work fine

---
 lib/BoundaryValueDiffEqCore/src/utils.jl      | 12 ++++--
 .../test/expanded/firk_basic_tests.jl         | 40 +++++--------------
 .../src/BoundaryValueDiffEqMIRKN.jl           |  3 +-
 .../src/algorithms.jl                         |  9 ++++-
 lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl     | 34 +++++++++-------
 .../src/multiple_shooting.jl                  |  3 +-
 .../src/single_shooting.jl                    |  1 +
 7 files changed, 49 insertions(+), 53 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 4437dfce6..5bbf29477 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -471,6 +471,12 @@ end
     return prob
 end
 
+@inline function __internal_optimization_problem(::SecondOrderBVProblem{uType, tType, iip},
+        args...; kwargs...) where {uType, tType, iip}
+    prob = OptimizationProblem(args...; kwargs...)
+    return prob
+end
+
 # Handling Initial Guesses
 """
     __extract_u0(u₀, t₀)
@@ -630,7 +636,7 @@ Constructs the internal problem based on the type of the boundary value problem
 algorithm used. It returns either a `NonlinearProblem` or an `OptimizationProblem`.
 """
 function __construct_internal_problem(
-        prob::BVProblem, alg, loss, jac, jac_prototype, resid_prototype, y, p, M, N)
+        prob::AbstractBVProblem, alg, loss, jac, jac_prototype, resid_prototype, y, p, M, N)
     T = eltype(y)
     iip = SciMLBase.isinplace(prob)
     if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize))
@@ -638,7 +644,7 @@ function __construct_internal_problem(
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
-        optf = OptimizationFunction{true}(__default_cost(prob.f), AutoFiniteDiff(), # Need to investigate the ForwardDiff dual problem
+        optf = OptimizationFunction{true}((x, p) -> 0.0, AutoFiniteDiff(), # Need to investigate the ForwardDiff dual problem
             cons = loss,
             cons_j = jac, cons_jac_prototype = jac_prototype)
         lcons = zeros(T, N*M)
@@ -658,7 +664,7 @@ function __construct_internal_problem(
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
         optf = OptimizationFunction{true}(
-            __default_cost(prob.f), get_dense_ad(diffmode), cons = loss,
+            (x, p) -> 0.0, get_dense_ad(diffmode), cons = loss,
             cons_j = jac, cons_jac_prototype = jac_prototype)
         lcons = zeros(T, N*M)
         ucons = zeros(T, N*M)
diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
index 1ccf7c928..1fa9d356e 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
@@ -307,30 +307,6 @@ end
     testTol = 1e-6
     nested = false
 
-    for stage in (2, 3, 5, 7)
-        s = Symbol("RadauIIa$(stage)")
-        @eval radau_solver(::Val{$stage}) = $(s)(
-            NewtonRaphson(), BVPJacobianAlgorithm(); nested)
-    end
-
-    for stage in (3, 4, 5)
-        s = Symbol("LobattoIIIa$(stage)")
-        @eval lobattoIIIa_solver(::Val{$stage}) = $(s)(
-            NewtonRaphson(), BVPJacobianAlgorithm(); nested)
-    end
-
-    for stage in (3, 4, 5)
-        s = Symbol("LobattoIIIb$(stage)")
-        @eval lobattoIIIb_solver(::Val{$stage}) = $(s)(
-            NewtonRaphson(), BVPJacobianAlgorithm(); nested)
-    end
-
-    for stage in (3, 4, 5)
-        s = Symbol("LobattoIIIc$(stage)")
-        @eval lobattoIIIc_solver(::Val{$stage}) = $(s)(
-            NewtonRaphson(), BVPJacobianAlgorithm(); nested)
-    end
-
     @testset "Radau interpolations" begin
         @testset "Interpolation tests for RadauIIa$stage" for stage in (2, 3, 5, 7)
             @time sol = solve(prob_bvp_linear, radau_solver(Val(stage)); dt = 0.001)
@@ -356,10 +332,11 @@ end
             for (id, lobatto_solver) in zip(("a", "b", "c"),
                 (lobattoIIIa_solver, lobattoIIIb_solver, lobattoIIIc_solver))
                 begin
-                    @testset "Interpolation tests for LobattoIII$(id)$stage" for stage in
-                                                                                 (3, 4, 5)
-                        @time sol = solve(
-                            prob_bvp_linear, lobatto_solver(Val(stage)); dt = 0.001)
+                    @testset "Interpolation tests for LobattoIII$(id)$stage" for stage in (
+                        2, 3, 4, 5)
+                        adaptive = ifelse(stage == 2, false, true) # LobattoIIIa2 is not adaptive
+                        @time sol = solve(prob_bvp_linear, lobatto_solver(Val(stage));
+                            dt = 0.001, adaptive = adaptive)
                         @test sol(0.001)≈[0.998687464, -1.312035941] atol=testTol
                         @test sol(0.001; idxs = [1, 2])≈[0.998687464, -1.312035941] atol=testTol
                         @test sol(0.001; idxs = 1)≈0.998687464 atol=testTol
@@ -368,9 +345,10 @@ end
 
                     @testset "Derivative Interpolation tests for lobatto$(id)$stage" for stage in
                                                                                          (
-                        3, 4, 5)
-                        @time sol = solve(
-                            prob_bvp_linear, lobatto_solver(Val(stage)); dt = 0.001)
+                        2, 3, 4, 5)
+                        adaptive = ifelse(stage == 2, false, true) # LobattoIIIa2 is not adaptive
+                        @time sol = solve(prob_bvp_linear, lobatto_solver(Val(stage));
+                            dt = 0.001, adaptive = adaptive)
                         sol_analytic = prob_bvp_linear_analytic(nothing, λ, 0.04)
                         dsol_analytic = prob_bvp_linear_analytic_derivative(
                             nothing, λ, 0.04)
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl b/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
index 702ecda96..fd90fa1d5 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/BoundaryValueDiffEqMIRKN.jl
@@ -19,7 +19,8 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                __default_sparse_ad, __default_nonsparse_ad, get_dense_ad,
                                concrete_jacobian_algorithm, __default_coloring_algorithm,
                                __default_sparsity_detector, interval, __split_kwargs,
-                               NoErrorControl
+                               NoErrorControl, __construct_internal_problem,
+                               __concrete_kwargs
 
 using ConcreteStructs: @concrete
 using DiffEqBase: DiffEqBase
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/algorithms.jl b/lib/BoundaryValueDiffEqMIRKN/src/algorithms.jl
index ede6be63c..de36e6cf5 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/algorithms.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/algorithms.jl
@@ -6,7 +6,7 @@ for order in (4, 6)
 
     @eval begin
         """
-              $($alg)(; nlsolve = NewtonRaphson(), jac_alg = BVPJacobianAlgorithm(),
+              $($alg)(; nlsolve = NewtonRaphson(), optimize = nothing, jac_alg = BVPJacobianAlgorithm(),
                       defect_threshold = 0.1, max_num_subintervals = 3000)
 
         $($order)th order Monotonic Implicit Runge Kutta Nyström method.
@@ -16,6 +16,10 @@ for order in (4, 6)
             - `nlsolve`: Internal Nonlinear solver. Any solver which conforms to the SciML
               `NonlinearProblem` interface can be used. Note that any autodiff argument for
               the solver will be ignored and a custom jacobian algorithm will be used.
+            - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML
+              `OptimizationProblem` interface can be used. Note that any autodiff argument for
+              the solver will be ignored and a custom jacobian algorithm will be used. Optimization
+              solvers should first be loaded to allow this functionality.
             - `jac_alg`: Jacobian Algorithm used for the nonlinear solver. Defaults to
               `BVPJacobianAlgorithm()`, which automatically decides the best algorithm to
               use based on the input types and problem type.
@@ -46,8 +50,9 @@ for order in (4, 6)
         }
         ```
         """
-        @kwdef struct $(alg){N, J <: BVPJacobianAlgorithm, T} <: AbstractMIRKN
+        @kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm, T} <: AbstractMIRKN
             nlsolve::N = nothing
+            optimize::O = nothing
             jac_alg::J = BVPJacobianAlgorithm()
             defect_threshold::T = 0.1
             max_num_subintervals::Int = 3000
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
index 425ae8e89..25cd9508b 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
@@ -21,16 +21,19 @@
     fᵢ₂_cache
     resid_size
     nlsolve_kwargs
+    optimize_kwargs
     kwargs
 end
 
 Base.eltype(::MIRKNCache{iip, T}) where {iip, T} = T
 
-function SciMLBase.__init(prob::SecondOrderBVProblem, alg::AbstractMIRKN; dt = 0.0,
-        adaptive = false, abstol = 1e-6, controller = NoErrorControl(),
-        nlsolve_kwargs = (; abstol = abstol), kwargs...)
+function SciMLBase.__init(prob::SecondOrderBVProblem, alg::AbstractMIRKN;
+        dt = 0.0, adaptive = false, abstol = 1e-6,
+        controller = NoErrorControl(), nlsolve_kwargs = (; abstol = abstol),
+        optimize_kwargs = (; abstol = abstol), kwargs...)
     @set! alg.jac_alg = concrete_jacobian_algorithm(alg.jac_alg, prob, alg)
     iip = isinplace(prob)
+    @assert (iip || isnothing(alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
     t₀, t₁ = prob.tspan
     ig, T, M, Nig, X = __extract_problem_details(prob; dt, check_positive_dt = true)
     mesh = __extract_mesh(prob.u0, t₀, t₁, Nig)
@@ -93,9 +96,10 @@ function SciMLBase.__init(prob::SecondOrderBVProblem, alg::AbstractMIRKN; dt = 0
     prob_ = !(prob.u0 isa AbstractArray) ? remake(prob; u0 = X) : prob
 
     return MIRKNCache{iip, T}(
-        alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type, prob.p, alg, TU,
-        bcresid_prototype, mesh, mesh_dt, k_discrete, y, y₀, residual, fᵢ_cache, fᵢ₂_cache,
-        resid_size, nlsolve_kwargs, (; abstol, dt, adaptive, controller, kwargs...))
+        alg_order(alg), stage, M, size(X), f, bc, prob_, prob.problem_type,
+        prob.p, alg, TU, bcresid_prototype, mesh, mesh_dt, k_discrete,
+        y, y₀, residual, fᵢ_cache, fᵢ₂_cache, resid_size, nlsolve_kwargs,
+        optimize_kwargs, (; abstol, dt, adaptive, controller, kwargs...))
 end
 
 function SciMLBase.solve!(cache::MIRKNCache{iip, T}) where {iip, T}
@@ -112,9 +116,11 @@ function SciMLBase.solve!(cache::MIRKNCache{iip, T}) where {iip, T}
 end
 
 function __perform_mirkn_iteration(cache::MIRKNCache)
-    nlprob::NonlinearProblem = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀))
-    nlsolve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve)
-    sol_nlprob = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true)
+    nlprob = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀))
+    solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize)
+    kwargs = __concrete_kwargs(
+        cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs)
+    sol_nlprob = solve(nlprob, solve_alg, kwargs...)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
     return sol_nlprob, sol_nlprob.retcode
@@ -227,9 +233,8 @@ function __construct_nlproblem(cache::MIRKNCache{iip}, y, loss_bc::BC, loss_coll
             cache_collocation, loss_bc, loss_collocation, L, cache.p)
     end
     resid_prototype = vcat(resid_bc, resid_collocation)
-    nlf = NonlinearFunction{iip}(
-        loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
-    __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p)
+    return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype,
+        resid_prototype, y, cache.p, cache.M, 2 * N)
 end
 
 function __construct_nlproblem(cache::MIRKNCache{iip}, y, loss_bc::BC, loss_collocation::C,
@@ -270,9 +275,8 @@ function __construct_nlproblem(cache::MIRKNCache{iip}, y, loss_bc::BC, loss_coll
     end
 
     resid_prototype = copy(resid)
-    nlf = NonlinearFunction{iip}(
-        loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
-    return __internal_nlsolve_problem(cache.prob, resid_prototype, y, nlf, y, cache.p)
+    return __construct_internal_problem(cache.prob, cache.alg, loss, jac, jac_prototype,
+        resid_prototype, y, cache.p, cache.M, 2 * N)
 end
 
 function __mirkn_2point_jacobian!(J, x, diffmode, diffcache, loss_fn::L, resid, p) where {L}
diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
index aa43fb877..0e723b29e 100644
--- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
@@ -7,6 +7,7 @@ function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwa
         throw(ArgumentError("Currently MultipleShooting only supports `EnsembleSerial` and \
                              `EnsembleThreads`!"))
     end
+    @assert (iip || isnothing(_alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
 
     ig, T, N, Nig, u0 = __extract_problem_details(prob; dt = 0.1)
     has_initial_guess = _unwrap_val(ig)
@@ -221,7 +222,7 @@ function __solve_nlproblem!(::StandardBVProblem, alg::MultipleShooting, bcresid_
     nlprob = __construct_internal_problem(prob, alg, loss_fn, jac_fn, jac_prototype,
         resid_prototype, u_at_nodes, prob.p, M, length(nodes))
     nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize)
-    kwargs = __concrete_kwargs(alg.nlsolve, alg.optimize, kwargs...)
+    kwargs = __concrete_kwargs(alg.nlsolve, alg.optimize, nlsolve_kwargs, optimize_kwargs)
     solve(nlprob, nlsolve_alg; kwargs...)
 
     return nothing
diff --git a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
index 13bb886c6..f4d1319fd 100644
--- a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
@@ -7,6 +7,7 @@ function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (;
         verbose && @warn "Initial guess provided, but will be ignored for Shooting."
         u0 = __extract_u0(prob.u0, prob.p, first(prob.tspan))
     end
+    @assert (iip || isnothing(alg_.optimize)) "Out-of-place constraints don't allow optimization solvers "
 
     alg = concretize_jacobian_algorithm(alg_, prob)
     (; diffmode) = alg.jac_alg

From 130da34b87ed84b9017617e02e82385053fc0dea Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Mon, 21 Jul 2025 21:42:16 +0800
Subject: [PATCH 15/50] All solvers should be done

---
 .../src/BoundaryValueDiffEqAscher.jl          |  4 ++--
 .../src/algorithms.jl                         |  6 ++++-
 lib/BoundaryValueDiffEqAscher/src/ascher.jl   | 24 +++++++++++--------
 lib/BoundaryValueDiffEqAscher/src/utils.jl    |  2 +-
 lib/BoundaryValueDiffEqCore/src/utils.jl      |  4 ++--
 5 files changed, 24 insertions(+), 16 deletions(-)

diff --git a/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl b/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl
index e7da643e5..387527271 100644
--- a/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl
+++ b/lib/BoundaryValueDiffEqAscher/src/BoundaryValueDiffEqAscher.jl
@@ -9,8 +9,8 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                __Fix3, __concrete_solve_algorithm,
                                __internal_nlsolve_problem, __vec, __vec_f, __vec_f!,
                                __vec_bc, __vec_bc!, __extract_mesh, get_dense_ad,
-                               __get_bcresid_prototype, __split_kwargs,
-                               __default_nonsparse_ad
+                               __get_bcresid_prototype, __split_kwargs, __concrete_kwargs,
+                               __default_nonsparse_ad, __construct_internal_problem
 
 using ConcreteStructs: @concrete
 using DiffEqBase: DiffEqBase
diff --git a/lib/BoundaryValueDiffEqAscher/src/algorithms.jl b/lib/BoundaryValueDiffEqAscher/src/algorithms.jl
index 7e9d522de..784579ffb 100644
--- a/lib/BoundaryValueDiffEqAscher/src/algorithms.jl
+++ b/lib/BoundaryValueDiffEqAscher/src/algorithms.jl
@@ -14,6 +14,9 @@ for stage in (1, 2, 3, 4, 5, 6, 7)
           - `nlsolve`: Internal Nonlinear solver. Any solver which conforms to the SciML
             `NonlinearProblem` interface can be used. Note that any autodiff argument for
             the solver will be ignored and a custom jacobian algorithm will be used.
+          - `optimize`: Internal Optimization solver. Any solver which conforms to the SciML
+            `OptimizationProblem` interface can be used. Note that any autodiff argument for
+            the solver will be ignored and a custom jacobian algorithm will be used.
           - `max_num_subintervals`: Number of maximal subintervals, default as 3000.
           - `zeta`: side condition points, should always be provided.
 
@@ -46,8 +49,9 @@ for stage in (1, 2, 3, 4, 5, 6, 7)
         }
         ```
         """
-        @kwdef struct $(alg){N, J <: BVPJacobianAlgorithm} <: AbstractAscher
+        @kwdef struct $(alg){N, O, J <: BVPJacobianAlgorithm} <: AbstractAscher
             nlsolve::N = nothing
+            optimize::O = nothing
             zeta::Vector{Float64} = nothing
             jac_alg::J = BVPJacobianAlgorithm()
             max_num_subintervals::Int = 3000
diff --git a/lib/BoundaryValueDiffEqAscher/src/ascher.jl b/lib/BoundaryValueDiffEqAscher/src/ascher.jl
index 5c280d26c..e5f66fcc4 100644
--- a/lib/BoundaryValueDiffEqAscher/src/ascher.jl
+++ b/lib/BoundaryValueDiffEqAscher/src/ascher.jl
@@ -39,6 +39,7 @@
     TU
     valstr
     nlsolve_kwargs
+    optimize_kwargs
     kwargs
 end
 
@@ -60,7 +61,8 @@ end
 
 function SciMLBase.__init(
         prob::BVProblem, alg::AbstractAscher; dt = 0.0, controller = GlobalErrorControl(),
-        adaptive = true, abstol = 1e-4, nlsolve_kwargs = (; abstol = abstol), kwargs...)
+        adaptive = true, abstol = 1e-4, nlsolve_kwargs = (; abstol = abstol),
+        optimize_kwargs = (; abstol = abstol), kwargs...)
     (; tspan, p) = prob
     _, T, ncy, n, u0 = __extract_problem_details(prob; dt, check_positive_dt = true)
     t₀, t₁ = tspan
@@ -147,9 +149,9 @@ function SciMLBase.__init(
     g = build_almost_block_diagonals(zeta, ncomp, mesh, T)
     cache = AscherCache{iip, T}(
         prob, f, jac, bc, bcjac, k, copy(mesh), mesh, mesh_dt, ncomp, ny, p, zeta,
-        fixpnt, alg, prob.problem_type, bcresid_prototype, residual, zval, yval,
-        gval, err, g, w, v, lz, ly, dmz, delz, deldmz, dqdmz, dmv, pvtg, pvtw, TU,
-        valst, nlsolve_kwargs, (; abstol, dt, adaptive, controller, kwargs...))
+        fixpnt, alg, prob.problem_type, bcresid_prototype, residual, zval, yval, gval,
+        err, g, w, v, lz, ly, dmz, delz, deldmz, dqdmz, dmv, pvtg, pvtw, TU, valst,
+        nlsolve_kwargs, optimize_kwargs, (; abstol, dt, adaptive, controller, kwargs...))
     return cache
 end
 
@@ -176,8 +178,10 @@ function __perform_ascher_iteration(
         cache::AscherCache{iip, T}, abstol, adaptive::Bool) where {iip, T}
     info::ReturnCode.T = ReturnCode.Success
     nlprob = __construct_nlproblem(cache)
-    nlsolve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve)
-    nlsol = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs...)
+    solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize)
+    kwargs = __concrete_kwargs(
+        cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs)
+    nlsol = solve(nlprob, solve_alg; kwargs...)
     error_norm = 2 * abstol
     info = nlsol.retcode
 
@@ -203,7 +207,7 @@ function __perform_ascher_iteration(
         __expand_cache_for_error!(cache)
 
         _nlprob = __construct_nlproblem(cache)
-        nlsol = __solve(_nlprob, nlsolve_alg; cache.nlsolve_kwargs...)
+        nlsol = solve(_nlprob, solve_alg; kwargs...)
 
         error_norm = error_estimate!(cache)
         if norm(error_norm) > abstol
@@ -346,9 +350,9 @@ function __construct_nlproblem(cache::AscherCache{iip, T}) where {iip, T}
             jac_prototype, u, diffmode, jac_cache, loss, cache.p)
     end
 
-    nlf = NonlinearFunction{iip}(
-        loss; jac = jac, resid_prototype = resid_prototype, jac_prototype = jac_prototype)
-    return __internal_nlsolve_problem(cache.prob, similar(lz), lz, nlf, lz, cache.p)
+    return __construct_internal_problem(
+        cache.prob, alg, loss, jac, jac_prototype, resid_prototype,
+        lz, cache.p, cache.ncomp, length(cache.mesh))
 end
 
 function __ascher_mpoint_jacobian!(J, x, diffmode, diffcache, loss, resid, p)
diff --git a/lib/BoundaryValueDiffEqAscher/src/utils.jl b/lib/BoundaryValueDiffEqAscher/src/utils.jl
index fb75d23d0..1676b4fd3 100644
--- a/lib/BoundaryValueDiffEqAscher/src/utils.jl
+++ b/lib/BoundaryValueDiffEqAscher/src/utils.jl
@@ -112,7 +112,7 @@ end
     return nothing
 end
 
-@views function recursive_flatten!(y::Vector, x::Vector{Vector{T}}) where {T}
+@views function recursive_flatten!(y::AbstractArray, x::Vector{Vector{T}}) where {T}
     i = 0
     for xᵢ in x
         copyto!(y[(i + 1):(i + length(xᵢ))], xᵢ)
diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 5bbf29477..92c3c3e1f 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -644,7 +644,7 @@ function __construct_internal_problem(
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
-        optf = OptimizationFunction{true}((x, p) -> 0.0, AutoFiniteDiff(), # Need to investigate the ForwardDiff dual problem
+        optf = OptimizationFunction{true}(__default_cost(prob.f), AutoFiniteDiff(), # Need to investigate the ForwardDiff dual problem
             cons = loss,
             cons_j = jac, cons_jac_prototype = jac_prototype)
         lcons = zeros(T, N*M)
@@ -664,7 +664,7 @@ function __construct_internal_problem(
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
         optf = OptimizationFunction{true}(
-            (x, p) -> 0.0, get_dense_ad(diffmode), cons = loss,
+            __default_cost(prob.f), get_dense_ad(diffmode), cons = loss,
             cons_j = jac, cons_jac_prototype = jac_prototype)
         lcons = zeros(T, N*M)
         ucons = zeros(T, N*M)

From d7e4ef4b8bf1b6a1949e613a37ce6e2cd35a4271 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Mon, 21 Jul 2025 21:55:29 +0800
Subject: [PATCH 16/50] Fix kwargs issue with shooting

---
 .../src/multiple_shooting.jl                       | 14 ++++++--------
 1 file changed, 6 insertions(+), 8 deletions(-)

diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
index 0e723b29e..dbbb449af 100644
--- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
@@ -66,14 +66,14 @@ function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwa
 
         if prob.problem_type isa TwoPointBVProblem
             __solve_nlproblem!(prob.problem_type, alg, bcresid_prototype, u_at_nodes, nodes,
-                cur_nshoot, M, N, resida_len, residb_len, solve_internal_odes!,
-                bc[1], bc[2], prob, u0, ode_cache_loss_fn, ensemblealg,
-                internal_ode_kwargs; verbose, kwargs..., nlsolve_kwargs...)
+                cur_nshoot, M, N, resida_len, residb_len, solve_internal_odes!, bc[1],
+                bc[2], prob, u0, ode_cache_loss_fn, ensemblealg, internal_ode_kwargs;
+                verbose, kwargs..., nlsolve_kwargs, optimize_kwargs...)
         else
             __solve_nlproblem!(prob.problem_type, alg, bcresid_prototype, u_at_nodes, nodes,
-                cur_nshoot, M, N, prod(resid_size), solve_internal_odes!,
-                bc, prob, f, u0_size, u0, ode_cache_loss_fn, ensemblealg,
-                internal_ode_kwargs; verbose, kwargs..., nlsolve_kwargs...)
+                cur_nshoot, M, N, prod(resid_size), solve_internal_odes!, bc, prob,
+                f, u0_size, u0, ode_cache_loss_fn, ensemblealg, internal_ode_kwargs;
+                verbose, kwargs..., nlsolve_kwargs, optimize_kwargs...)
         end
     end
 
@@ -147,7 +147,6 @@ function __solve_nlproblem!(
         resid_prototype, u_at_nodes, prob.p, M, length(nodes))
 
     nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize)
-    kwargs = __concrete_kwargs(alg.nlsolve, alg.optimize, kwargs...)
     solve(nlprob, nlsolve_alg; kwargs...)
 
     return nothing
@@ -222,7 +221,6 @@ function __solve_nlproblem!(::StandardBVProblem, alg::MultipleShooting, bcresid_
     nlprob = __construct_internal_problem(prob, alg, loss_fn, jac_fn, jac_prototype,
         resid_prototype, u_at_nodes, prob.p, M, length(nodes))
     nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize)
-    kwargs = __concrete_kwargs(alg.nlsolve, alg.optimize, nlsolve_kwargs, optimize_kwargs)
     solve(nlprob, nlsolve_alg; kwargs...)
 
     return nothing

From 9901afbb4e49d560c7ab2656715b7cc83e31c3e1 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Mon, 21 Jul 2025 22:31:55 +0800
Subject: [PATCH 17/50] Fix inplace assert in multiple shooting

---
 lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
index dbbb449af..89bcf2f82 100644
--- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
@@ -7,7 +7,6 @@ function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwa
         throw(ArgumentError("Currently MultipleShooting only supports `EnsembleSerial` and \
                              `EnsembleThreads`!"))
     end
-    @assert (iip || isnothing(_alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
 
     ig, T, N, Nig, u0 = __extract_problem_details(prob; dt = 0.1)
     has_initial_guess = _unwrap_val(ig)
@@ -16,6 +15,7 @@ function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwa
 
     bcresid_prototype, resid_size = __get_bcresid_prototype(prob, u0)
     iip, bc, u0, u0_size = isinplace(prob), prob.f.bc, deepcopy(u0), size(u0)
+    @assert (iip || isnothing(_alg.optimize)) "Out-of-place constraints don't allow optimization solvers "
 
     __alg = concretize_jacobian_algorithm(_alg, prob)
     alg = if has_initial_guess && Nig != __alg.nshoots

From f710cd7494dc5f56c9cac6aff1fe85cce3920344 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Mon, 21 Jul 2025 22:40:26 +0800
Subject: [PATCH 18/50] Fix inplace assert in single shooting

---
 lib/BoundaryValueDiffEqShooting/src/single_shooting.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
index f4d1319fd..65180c1f9 100644
--- a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
@@ -7,13 +7,13 @@ function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (;
         verbose && @warn "Initial guess provided, but will be ignored for Shooting."
         u0 = __extract_u0(prob.u0, prob.p, first(prob.tspan))
     end
-    @assert (iip || isnothing(alg_.optimize)) "Out-of-place constraints don't allow optimization solvers "
 
     alg = concretize_jacobian_algorithm(alg_, prob)
     (; diffmode) = alg.jac_alg
 
     bcresid_prototype, resid_size = __get_bcresid_prototype(prob, u0)
     iip, bc, u0, u0_size = isinplace(prob), prob.f.bc, deepcopy(u0), size(u0)
+    @assert (iip || isnothing(alg_.optimize)) "Out-of-place constraints don't allow optimization solvers "
     resid_prototype = __vec(bcresid_prototype)
 
     # Construct the residual function

From 81da921f48d4e7d29eaa05444d1dcdbf083a9222 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Tue, 22 Jul 2025 22:47:30 +0800
Subject: [PATCH 19/50] Use correct setup in FIRK testing

---
 lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
index 3bd0607f1..d218f5c46 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
@@ -273,7 +273,7 @@ end
     @test_nowarn solve(bvp1, RadauIIa7(nl_solve, jac_alg; nested); dt = 0.05)
 end
 
-@testitem "Interpolation" begin
+@testitem "Interpolation" setup=[FIRKExpandedConvergenceTests] begin
     using LinearAlgebra
 
     λ = 1

From 4d677c7436b9873cc5bcd6c2085d55d041611ea7 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Thu, 24 Jul 2025 00:24:29 +0800
Subject: [PATCH 20/50] Unify testing environment

---
 lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
index 25848c684..b6b2a739f 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
@@ -294,7 +294,7 @@ end
     @test_nowarn solve(bvp1, RadauIIa7(nl_solve, jac_alg; nested); dt = 0.05)
 end
 
-@testitem "Interpolation" begin
+@testitem "Interpolation" setup=[FIRKNestedConvergenceTests] begin
     using LinearAlgebra
 
     λ = 1

From 66a6b899e93ac36e3ca0accf4a90b5a56913f33c Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Thu, 24 Jul 2025 00:49:36 +0800
Subject: [PATCH 21/50] Fix ensemble test for FIRK

---
 lib/BoundaryValueDiffEqFIRK/src/algorithms.jl | 44 -------------------
 .../test/expanded/ensemble_tests.jl           |  8 ++--
 .../test/nested/ensemble_tests.jl             |  8 ++--
 .../test/nested/firk_basic_tests.jl           | 24 ----------
 4 files changed, 8 insertions(+), 76 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl b/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl
index ec957fa4d..1220506c7 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl
@@ -99,17 +99,6 @@ for stage in (1, 2, 3, 5, 7)
             defect_threshold::T = 0.1
             max_num_subintervals::Int = 3000
         end
-        $(alg)(nlsolve::N,
-            optimize::O,
-            jac_alg::J;
-            nested_nlsolve = false,
-            nested_nlsolve_kwargs::NamedTuple = (;),
-            defect_threshold::T = 0.1,
-            max_num_subintervals::Int = 3000) where {N,
-            O,
-            J,
-            T} = $(alg){N, O, J, T}(nlsolve, optimize, jac_alg, nested_nlsolve,
-            nested_nlsolve_kwargs, defect_threshold, max_num_subintervals)
     end
 end
 
@@ -211,17 +200,6 @@ for stage in (2, 3, 4, 5)
             defect_threshold::T = 0.1
             max_num_subintervals::Int = 3000
         end
-        $(alg)(nlsolve::N,
-            optimize::O,
-            jac_alg::J;
-            nested_nlsolve = false,
-            nested_nlsolve_kwargs::NamedTuple = (;),
-            defect_threshold::T = 0.1,
-            max_num_subintervals::Int = 3000) where {N,
-            O,
-            J,
-            T} = $(alg){N, O, J, T}(nlsolve, optimize, jac_alg, nested_nlsolve,
-            nested_nlsolve_kwargs, defect_threshold, max_num_subintervals)
     end
 end
 
@@ -324,17 +302,6 @@ for stage in (2, 3, 4, 5)
             defect_threshold::T = 0.1
             max_num_subintervals::Int = 3000
         end
-        $(alg)(nlsolve::N,
-            optimize::O,
-            jac_alg::J;
-            nested_nlsolve = false,
-            nested_nlsolve_kwargs::NamedTuple = (;),
-            defect_threshold::T = 0.1,
-            max_num_subintervals::Int = 3000) where {N,
-            O,
-            J,
-            T} = $(alg){N, O, J, T}(nlsolve, optimize, jac_alg, nested_nlsolve,
-            nested_nlsolve_kwargs, defect_threshold, max_num_subintervals)
     end
 end
 
@@ -437,17 +404,6 @@ for stage in (2, 3, 4, 5)
             defect_threshold::T = 0.1
             max_num_subintervals::Int = 3000
         end
-        $(alg)(nlsolve::N,
-            optimize::O,
-            jac_alg::J;
-            nested_nlsolve = false,
-            nested_nlsolve_kwargs::NamedTuple = (;),
-            defect_threshold::T = 0.1,
-            max_num_subintervals::Int = 3000) where {N,
-            O,
-            J,
-            T} = $(alg){N, O, J, T}(nlsolve, optimize, jac_alg, nested_nlsolve,
-            nested_nlsolve_kwargs, defect_threshold, max_num_subintervals)
     end
 end
 
diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
index 47d4656db..bce665b95 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
@@ -27,7 +27,7 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested);
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested);
                 trajectories = 10, dt = 0.1)
             @test sol.converged
         end
@@ -40,7 +40,7 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested);
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg; nested_nlsolve = nested);
                 trajectories = 10, dt = 0.1)
             @test sol.converged
         end
@@ -52,7 +52,7 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested);
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested);
                 trajectories = 10, dt = 0.1)
             @test sol.converged
         end
@@ -64,7 +64,7 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested);
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested);
                 trajectories = 10, dt = 0.1)
             @test sol.converged
         end
diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl
index 65ff128d0..6f565b2da 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/nested/ensemble_tests.jl
@@ -27,7 +27,7 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested);
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested);
                 trajectories = 10, dt = 0.1)
             @test sol.converged
         end
@@ -40,7 +40,7 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested);
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested);
                 trajectories = 10, dt = 0.1)
             @test sol.converged
         end
@@ -52,7 +52,7 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested);
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested);
                 trajectories = 10, dt = 0.1)
             @test sol.converged
         end
@@ -64,7 +64,7 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested);
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested);
                 trajectories = 10, dt = 0.1)
             @test sol.converged
         end
diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
index b6b2a739f..8b9d6c273 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
@@ -328,30 +328,6 @@ end
     testTol = 1e-6
     nested = true
 
-    for stage in (2, 3, 5, 7)
-        s = Symbol("RadauIIa$(stage)")
-        @eval radau_solver(::Val{$stage}) = $(s)(
-            NewtonRaphson(), BVPJacobianAlgorithm(); nested)
-    end
-
-    for stage in (3, 4, 5)
-        s = Symbol("LobattoIIIa$(stage)")
-        @eval lobattoIIIa_solver(::Val{$stage}) = $(s)(
-            NewtonRaphson(), BVPJacobianAlgorithm(); nested)
-    end
-
-    for stage in (3, 4, 5)
-        s = Symbol("LobattoIIIb$(stage)")
-        @eval lobattoIIIb_solver(::Val{$stage}) = $(s)(
-            NewtonRaphson(), BVPJacobianAlgorithm(); nested)
-    end
-
-    for stage in (3, 4, 5)
-        s = Symbol("LobattoIIIc$(stage)")
-        @eval lobattoIIIc_solver(::Val{$stage}) = $(s)(
-            NewtonRaphson(), BVPJacobianAlgorithm(); nested)
-    end
-
     @testset "Radau interpolations" begin
         @testset "Interpolation tests for RadauIIa$stage" for stage in (2, 3, 5, 7)
             @time sol = solve(prob_bvp_linear, radau_solver(Val(stage)); dt = 0.001)

From 44d556f8568d7f3773243e666280fc7277dcc008 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Thu, 24 Jul 2025 16:04:25 +0800
Subject: [PATCH 22/50] Fix more FIRK tests

---
 .../test/expanded/ensemble_tests.jl           |  2 +-
 .../test/expanded/firk_basic_tests.jl         | 48 ++++++++++++-------
 2 files changed, 32 insertions(+), 18 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
index bce665b95..68b72b692 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
@@ -40,7 +40,7 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg; nested_nlsolve = nested);
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested);
                 trajectories = 10, dt = 0.1)
             @test sol.converged
         end
diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
index d218f5c46..fb53606e4 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
@@ -248,29 +248,43 @@ end
     nested = false
 
     # Using ForwardDiff might lead to Cache expansion warnings
-    @test_nowarn solve(bvp1, LobattoIIIa2(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa5(nl_solve, jac_alg; nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIa2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIa3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIa4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIa5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
+    @test_nowarn solve(bvp1, LobattoIIIb2(; nl_solve, jac_alg, nested_nlsolve = nested);
+        dt = 0.005, adaptive = false)
+    @test_nowarn solve(
+        bvp1, LobattoIIIb3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIb4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
     @test_nowarn solve(
-        bvp1, LobattoIIIb2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, LobattoIIIb3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIb4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIb5(nl_solve, jac_alg; nested); dt = 0.005)
+        bvp1, LobattoIIIb5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
+    @test_nowarn solve(bvp1, LobattoIIIc2(; nl_solve, jac_alg, nested_nlsolve = nested);
+        dt = 0.005, adaptive = false)
     @test_nowarn solve(
-        bvp1, LobattoIIIc2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, LobattoIIIc3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIc4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIc5(nl_solve, jac_alg; nested); dt = 0.005)
+        bvp1, LobattoIIIc3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIc4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIc5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
+    @test_nowarn solve(bvp1, RadauIIa1(; nl_solve, jac_alg, nested_nlsolve = nested);
+        dt = 0.005, adaptive = false)
+    @test_nowarn solve(
+        bvp1, RadauIIa2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, RadauIIa3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, RadauIIa5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.05)
     @test_nowarn solve(
-        bvp1, RadauIIa1(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, RadauIIa2(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, RadauIIa3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, RadauIIa5(nl_solve, jac_alg; nested); dt = 0.05)
-    @test_nowarn solve(bvp1, RadauIIa7(nl_solve, jac_alg; nested); dt = 0.05)
+        bvp1, RadauIIa7(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.05)
 end
 
 @testitem "Interpolation" setup=[FIRKExpandedConvergenceTests] begin

From 2237f7043fee2daf51a0cc295c3eb6a7f902095d Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Thu, 24 Jul 2025 21:58:41 +0800
Subject: [PATCH 23/50] Fix FIRK on simple pendulum

---
 .../test/expanded/firk_basic_tests.jl         | 36 +++++++++----------
 1 file changed, 18 insertions(+), 18 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
index fb53606e4..d0391b666 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
@@ -244,47 +244,47 @@ end
 
     jac_alg = BVPJacobianAlgorithm(;
         bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff()))
-    nl_solve = NewtonRaphson()
+    nlsolve = NewtonRaphson()
     nested = false
 
     # Using ForwardDiff might lead to Cache expansion warnings
     @test_nowarn solve(
-        bvp1, LobattoIIIa2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIa2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
     @test_nowarn solve(
-        bvp1, LobattoIIIa3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIa3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
     @test_nowarn solve(
-        bvp1, LobattoIIIa4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIa4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
     @test_nowarn solve(
-        bvp1, LobattoIIIa5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIa5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
-    @test_nowarn solve(bvp1, LobattoIIIb2(; nl_solve, jac_alg, nested_nlsolve = nested);
+    @test_nowarn solve(bvp1, LobattoIIIb2(; nlsolve, jac_alg, nested_nlsolve = nested);
         dt = 0.005, adaptive = false)
     @test_nowarn solve(
-        bvp1, LobattoIIIb3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIb3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
     @test_nowarn solve(
-        bvp1, LobattoIIIb4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIb4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
     @test_nowarn solve(
-        bvp1, LobattoIIIb5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIb5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
-    @test_nowarn solve(bvp1, LobattoIIIc2(; nl_solve, jac_alg, nested_nlsolve = nested);
+    @test_nowarn solve(bvp1, LobattoIIIc2(; nlsolve, jac_alg, nested_nlsolve = nested);
         dt = 0.005, adaptive = false)
     @test_nowarn solve(
-        bvp1, LobattoIIIc3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIc3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
     @test_nowarn solve(
-        bvp1, LobattoIIIc4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIc4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
     @test_nowarn solve(
-        bvp1, LobattoIIIc5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIc5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
-    @test_nowarn solve(bvp1, RadauIIa1(; nl_solve, jac_alg, nested_nlsolve = nested);
+    @test_nowarn solve(bvp1, RadauIIa1(; nlsolve, jac_alg, nested_nlsolve = nested);
         dt = 0.005, adaptive = false)
     @test_nowarn solve(
-        bvp1, RadauIIa2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, RadauIIa2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
     @test_nowarn solve(
-        bvp1, RadauIIa3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, RadauIIa3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
     @test_nowarn solve(
-        bvp1, RadauIIa5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.05)
+        bvp1, RadauIIa5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.05)
     @test_nowarn solve(
-        bvp1, RadauIIa7(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.05)
+        bvp1, RadauIIa7(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.05)
 end
 
 @testitem "Interpolation" setup=[FIRKExpandedConvergenceTests] begin

From 318f0159a80756766fd355db2bb98551b71f3b05 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Fri, 25 Jul 2025 02:18:45 +0800
Subject: [PATCH 24/50] Set default abstol for nonlinear solving in shooting

---
 lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl | 2 +-
 lib/BoundaryValueDiffEqShooting/src/single_shooting.jl   | 3 ++-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
index 758df5638..e3a6782f6 100644
--- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
@@ -1,5 +1,5 @@
 function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwargs = (;),
-        nlsolve_kwargs = (;), optimize_kwargs = (;),
+        nlsolve_kwargs = (; abstol = 1e-6), optimize_kwargs = (; abstol = 1e-6),
         ensemblealg = EnsembleThreads(), verbose = true, kwargs...)
     (; f, tspan) = prob
 
diff --git a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
index 65180c1f9..ecc0fce3a 100644
--- a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
@@ -1,5 +1,6 @@
 function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (;),
-        nlsolve_kwargs = (;), optimize_kwargs = (;), verbose = true, kwargs...)
+        nlsolve_kwargs = (; abstol = 1e-6),
+        optimize_kwargs = (; abstol = 1e-6), verbose = true, kwargs...)
     # Setup the problem
     if prob.u0 isa AbstractArray{<:Number}
         u0 = prob.u0

From ebc6dcaaca07f163d6f0eef994c0523491ed2a02 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Fri, 25 Jul 2025 02:46:01 +0800
Subject: [PATCH 25/50] Fix kwargs issue with Shooting

---
 lib/BoundaryValueDiffEqCore/src/utils.jl           | 14 ++++++++++----
 .../src/multiple_shooting.jl                       |  5 +++--
 .../src/single_shooting.jl                         |  6 +++---
 3 files changed, 16 insertions(+), 9 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index a3e2d8a7d..617de2ced 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -621,11 +621,17 @@ function __split_kwargs(; abstol, adaptive, controller, kwargs...)
     return ((abstol, adaptive, controller), (; abstol, adaptive, kwargs...))
 end
 
-@inline __concrete_kwargs(nlsolve, ::Nothing, nlsolve_kwargs,
-    optimize_kwargs) = (nlsolve_kwargs..., alias_u0 = true)
+@inline __concrete_kwargs(nlsolve,
+    ::Nothing,
+    nlsolve_kwargs,
+    optimize_kwargs) = (
+    nlsolve_kwargs..., alias = SciMLBase.NonlinearAliasSpecifier(alias_u0 = true))
 @inline __concrete_kwargs(::Nothing, optimize, nlsolve_kwargs, optimize_kwargs) = (;) # Doesn't support for now
-@inline __concrete_kwargs(::Nothing, ::Nothing, nlsolve_kwargs,
-    optimize_kwargs) = (nlsolve_kwargs..., alias_u0 = true)
+@inline __concrete_kwargs(::Nothing,
+    ::Nothing,
+    nlsolve_kwargs,
+    optimize_kwargs) = (
+    nlsolve_kwargs..., alias = SciMLBase.NonlinearAliasSpecifier(alias_u0 = true))
 
 ## Optimization solver related utils ##
 
diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
index e3a6782f6..83a74f62b 100644
--- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
@@ -1,5 +1,6 @@
-function SciMLBase.__solve(prob::BVProblem, _alg::MultipleShooting; odesolve_kwargs = (;),
-        nlsolve_kwargs = (; abstol = 1e-6), optimize_kwargs = (; abstol = 1e-6),
+function SciMLBase.__solve(
+        prob::BVProblem, _alg::MultipleShooting; abstol = 1e-6, odesolve_kwargs = (;),
+        nlsolve_kwargs = (; abstol = abstol), optimize_kwargs = (; abstol = abstol),
         ensemblealg = EnsembleThreads(), verbose = true, kwargs...)
     (; f, tspan) = prob
 
diff --git a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
index ecc0fce3a..3c8294111 100644
--- a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
@@ -1,6 +1,6 @@
-function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; odesolve_kwargs = (;),
-        nlsolve_kwargs = (; abstol = 1e-6),
-        optimize_kwargs = (; abstol = 1e-6), verbose = true, kwargs...)
+function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; abstol = 1e-6,
+        odesolve_kwargs = (;), nlsolve_kwargs = (; abstol = abstol),
+        optimize_kwargs = (; abstol = abstol), verbose = true, kwargs...)
     # Setup the problem
     if prob.u0 isa AbstractArray{<:Number}
         u0 = prob.u0

From 3137c0b3d467376766437a23b497072c02fa2e9e Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Fri, 25 Jul 2025 03:44:53 +0800
Subject: [PATCH 26/50] Fix multiple shooting only use inplace problem
 constructor

---
 lib/BoundaryValueDiffEqCore/src/utils.jl      | 30 +++++++++++++++----
 .../src/multiple_shooting.jl                  | 10 ++++---
 2 files changed, 31 insertions(+), 9 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 617de2ced..2659d117c 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -645,9 +645,10 @@ end
 Constructs the internal problem based on the type of the boundary value problem and the
 algorithm used. It returns either a `NonlinearProblem` or an `OptimizationProblem`.
 """
-function __construct_internal_problem(
-        prob::AbstractBVProblem, alg, loss, jac, jac_prototype, resid_prototype, y, p, M, N)
+function __construct_internal_problem(prob::AbstractBVProblem, alg, loss, jac,
+        jac_prototype, resid_prototype, y, p, M::Int, N::Int)
     T = eltype(y)
+    # multiple shooting always use iip
     iip = SciMLBase.isinplace(prob)
     if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize))
         nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype,
@@ -664,8 +665,8 @@ function __construct_internal_problem(
     end
 end
 
-function __construct_internal_problem(
-        prob::TwoPointBVProblem, alg, loss, jac, jac_prototype, resid_prototype, y, p, M, N)
+function __construct_internal_problem(prob::TwoPointBVProblem, alg, loss, jac,
+        jac_prototype, resid_prototype, y, p, M::Int, N::Int)
     T = eltype(y)
     iip = SciMLBase.isinplace(prob)
     if !isnothing(alg.nlsolve)
@@ -680,6 +681,25 @@ function __construct_internal_problem(
         ucons = zeros(T, N*M)
 
         return __internal_optimization_problem(
-            pro, optf, y, p; lcons = lcons, ucons = ucons)
+            prob, optf, y, p; lcons = lcons, ucons = ucons)
+    end
+end
+# Multiple shooting only use inplace version internal problem constructor
+function __construct_internal_problem(prob, alg, loss, jac, jac_prototype,
+        resid_prototype, y, p, M::Int, N::Int, ::Nothing)
+    T = eltype(y)
+    if !isnothing(alg.nlsolve)
+        nlf = NonlinearFunction{true}(loss; jac = jac, resid_prototype = resid_prototype,
+            jac_prototype = jac_prototype)
+        return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
+    else
+        optf = OptimizationFunction{true}(
+            __default_cost(prob.f), get_dense_ad(diffmode), cons = loss,
+            cons_j = jac, cons_jac_prototype = jac_prototype)
+        lcons = zeros(T, N*M)
+        ucons = zeros(T, N*M)
+
+        return __internal_optimization_problem(
+            prob, optf, y, p; lcons = lcons, ucons = ucons)
     end
 end
diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
index 83a74f62b..d6d6f2a01 100644
--- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
@@ -144,8 +144,9 @@ function __solve_nlproblem!(
         jac_prototype = jac_prototype)
 
     # NOTE: u_at_nodes is updated inplace
-    nlprob = __construct_internal_problem(prob, alg, loss_fn, jac_fn, jac_prototype,
-        resid_prototype, u_at_nodes, prob.p, M, length(nodes))
+    nlprob = __construct_internal_problem(
+        prob, alg, loss_fn, jac_fn, jac_prototype, resid_prototype,
+        u_at_nodes, prob.p, M, length(nodes), nothing)
 
     nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize)
     solve(nlprob, nlsolve_alg; kwargs...)
@@ -219,8 +220,9 @@ function __solve_nlproblem!(::StandardBVProblem, alg::MultipleShooting, bcresid_
         ode_fn, bc_fn, nonbc_diffmode, bc_diffmode, N, M, __cache_trait(alg.jac_alg))
 
     # NOTE: u_at_nodes is updated inplace
-    nlprob = __construct_internal_problem(prob, alg, loss_fn, jac_fn, jac_prototype,
-        resid_prototype, u_at_nodes, prob.p, M, length(nodes))
+    nlprob = __construct_internal_problem(
+        prob, alg, loss_fn, jac_fn, jac_prototype, resid_prototype,
+        u_at_nodes, prob.p, M, length(nodes), nothing)
     nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize)
     solve(nlprob, nlsolve_alg; kwargs...)
 

From 3d0dc521784b3d63d90bbb7d56c6bc53d13a7949 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Fri, 25 Jul 2025 03:58:16 +0800
Subject: [PATCH 27/50] Fix typo

---
 lib/BoundaryValueDiffEqCore/src/utils.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 2659d117c..f4d13b379 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -675,8 +675,8 @@ function __construct_internal_problem(prob::TwoPointBVProblem, alg, loss, jac,
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
         optf = OptimizationFunction{true}(
-            __default_cost(prob.f), get_dense_ad(diffmode), cons = loss,
-            cons_j = jac, cons_jac_prototype = jac_prototype)
+            __default_cost(prob.f), get_dense_ad(alg.diffmode),
+            cons = loss, cons_j = jac, cons_jac_prototype = jac_prototype)
         lcons = zeros(T, N*M)
         ucons = zeros(T, N*M)
 
@@ -694,8 +694,8 @@ function __construct_internal_problem(prob, alg, loss, jac, jac_prototype,
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
         optf = OptimizationFunction{true}(
-            __default_cost(prob.f), get_dense_ad(diffmode), cons = loss,
-            cons_j = jac, cons_jac_prototype = jac_prototype)
+            __default_cost(prob.f), get_dense_ad(alg.diffmode),
+            cons = loss, cons_j = jac, cons_jac_prototype = jac_prototype)
         lcons = zeros(T, N*M)
         ucons = zeros(T, N*M)
 

From f16dc9a4b7dbeb00ac38a6553e6db220e267eae3 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <52615090+ErikQQY@users.noreply.github.com>
Date: Fri, 25 Jul 2025 04:09:17 +0800
Subject: [PATCH 28/50] Use proper Jacobian algorithm

---
 lib/BoundaryValueDiffEqCore/src/utils.jl | 23 +++++++++++++++++++++--
 1 file changed, 21 insertions(+), 2 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index f4d13b379..a6b6a1e2c 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -675,7 +675,7 @@ function __construct_internal_problem(prob::TwoPointBVProblem, alg, loss, jac,
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
         optf = OptimizationFunction{true}(
-            __default_cost(prob.f), get_dense_ad(alg.diffmode),
+            __default_cost(prob.f), get_dense_ad(alg.jac_alg.diffmode),
             cons = loss, cons_j = jac, cons_jac_prototype = jac_prototype)
         lcons = zeros(T, N*M)
         ucons = zeros(T, N*M)
@@ -694,7 +694,26 @@ function __construct_internal_problem(prob, alg, loss, jac, jac_prototype,
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
     else
         optf = OptimizationFunction{true}(
-            __default_cost(prob.f), get_dense_ad(alg.diffmode),
+            __default_cost(prob.f), get_dense_ad(alg.jac_alg.diffmode),
+            cons = loss, cons_j = jac, cons_jac_prototype = jac_prototype)
+        lcons = zeros(T, N*M)
+        ucons = zeros(T, N*M)
+
+        return __internal_optimization_problem(
+            prob, optf, y, p; lcons = lcons, ucons = ucons)
+    end
+end
+function __construct_internal_problem(prob::TwoPointBVProblem, alg, loss, jac,
+        jac_prototype, resid_prototype, y, p, M::Int, N::Int,::Nothing)
+    T = eltype(y)
+    iip = SciMLBase.isinplace(prob)
+    if !isnothing(alg.nlsolve)
+        nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype,
+            jac_prototype = jac_prototype)
+        return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
+    else
+        optf = OptimizationFunction{true}(
+            __default_cost(prob.f), get_dense_ad(alg.jac_alg.nonbc_diffmode),
             cons = loss, cons_j = jac, cons_jac_prototype = jac_prototype)
         lcons = zeros(T, N*M)
         ucons = zeros(T, N*M)

From ed84cb0ac0f0cc02ed010e9cf10a38b44bd9ecd6 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Fri, 25 Jul 2025 11:30:08 +0800
Subject: [PATCH 29/50] Fix internal solver

---
 lib/BoundaryValueDiffEqCore/src/utils.jl | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index a6b6a1e2c..93212b764 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -669,7 +669,7 @@ function __construct_internal_problem(prob::TwoPointBVProblem, alg, loss, jac,
         jac_prototype, resid_prototype, y, p, M::Int, N::Int)
     T = eltype(y)
     iip = SciMLBase.isinplace(prob)
-    if !isnothing(alg.nlsolve)
+    if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize))
         nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype,
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
@@ -688,7 +688,7 @@ end
 function __construct_internal_problem(prob, alg, loss, jac, jac_prototype,
         resid_prototype, y, p, M::Int, N::Int, ::Nothing)
     T = eltype(y)
-    if !isnothing(alg.nlsolve)
+    if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize))
         nlf = NonlinearFunction{true}(loss; jac = jac, resid_prototype = resid_prototype,
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
@@ -703,11 +703,12 @@ function __construct_internal_problem(prob, alg, loss, jac, jac_prototype,
             prob, optf, y, p; lcons = lcons, ucons = ucons)
     end
 end
-function __construct_internal_problem(prob::TwoPointBVProblem, alg, loss, jac,
-        jac_prototype, resid_prototype, y, p, M::Int, N::Int,::Nothing)
+function __construct_internal_problem(
+        prob::TwoPointBVProblem, alg, loss, jac, jac_prototype,
+        resid_prototype, y, p, M::Int, N::Int, ::Nothing)
     T = eltype(y)
     iip = SciMLBase.isinplace(prob)
-    if !isnothing(alg.nlsolve)
+    if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize))
         nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype,
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)

From 15eb2def1ba9764300e0e7940c130c4daacdbd41 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 26 Jul 2025 01:20:36 +0800
Subject: [PATCH 30/50] Fix nested FIRK testing

---
 .../test/nested/firk_basic_tests.jl           | 48 ++++++++++++-------
 1 file changed, 31 insertions(+), 17 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
index 8b9d6c273..18cf6d3b5 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
@@ -269,29 +269,43 @@ end
     nested = true
 
     # Using ForwardDiff might lead to Cache expansion warnings
-    @test_nowarn solve(bvp1, LobattoIIIa2(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa5(nl_solve, jac_alg; nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIa2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIa3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIa4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIa5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
+    @test_nowarn solve(bvp1, LobattoIIIb2(; nl_solve, jac_alg, nested_nlsolve = nested);
+        dt = 0.005, adaptive = false)
+    @test_nowarn solve(
+        bvp1, LobattoIIIb3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIb4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
     @test_nowarn solve(
-        bvp1, LobattoIIIb2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, LobattoIIIb3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIb4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIb5(nl_solve, jac_alg; nested); dt = 0.005)
+        bvp1, LobattoIIIb5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
+    @test_nowarn solve(bvp1, LobattoIIIc2(; nl_solve, jac_alg, nested_nlsolve = nested);
+        dt = 0.005, adaptive = false)
     @test_nowarn solve(
-        bvp1, LobattoIIIc2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, LobattoIIIc3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIc4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIc5(nl_solve, jac_alg; nested); dt = 0.005)
+        bvp1, LobattoIIIc3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIc4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, LobattoIIIc5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
+    @test_nowarn solve(bvp1, RadauIIa1(; nl_solve, jac_alg, nested_nlsolve = nested);
+        dt = 0.005, adaptive = false)
+    @test_nowarn solve(
+        bvp1, RadauIIa2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, RadauIIa3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(
+        bvp1, RadauIIa5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.05)
     @test_nowarn solve(
-        bvp1, RadauIIa1(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, RadauIIa2(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, RadauIIa3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, RadauIIa5(nl_solve, jac_alg; nested); dt = 0.05)
-    @test_nowarn solve(bvp1, RadauIIa7(nl_solve, jac_alg; nested); dt = 0.05)
+        bvp1, RadauIIa7(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.05)
 end
 
 @testitem "Interpolation" setup=[FIRKNestedConvergenceTests] begin

From db1bddd2ab63bfb64b3e74f0e4f23ef17ab92a13 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 26 Jul 2025 02:12:28 +0800
Subject: [PATCH 31/50] Fix shooting in manifold tests

---
 test/misc/manifolds_tests.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/misc/manifolds_tests.jl b/test/misc/manifolds_tests.jl
index 9ae55eb3b..8bb061fab 100644
--- a/test/misc/manifolds_tests.jl
+++ b/test/misc/manifolds_tests.jl
@@ -65,7 +65,7 @@
 
         for alg in algs
             if alg isa Shooting || alg isa MultipleShooting
-                sol = solve(bvp, alg)
+                sol = solve(bvp, alg; abstol = 1e-8)
             else
                 sol = solve(bvp, alg; dt, abstol = 1e-8)
             end

From ac52720a289801bc38c4bba46192db9d4f73dc97 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sun, 27 Jul 2025 01:19:28 +0800
Subject: [PATCH 32/50] Remove some stale deps in shooting

---
 lib/BoundaryValueDiffEqShooting/Project.toml                | 6 ------
 .../src/BoundaryValueDiffEqShooting.jl                      | 4 ----
 2 files changed, 10 deletions(-)

diff --git a/lib/BoundaryValueDiffEqShooting/Project.toml b/lib/BoundaryValueDiffEqShooting/Project.toml
index ea5e83b11..6733e5000 100644
--- a/lib/BoundaryValueDiffEqShooting/Project.toml
+++ b/lib/BoundaryValueDiffEqShooting/Project.toml
@@ -10,13 +10,10 @@ BoundaryValueDiffEqCore = "56b672f2-a5fe-4263-ab2d-da677488eb3a"
 ConcreteStructs = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
-FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e"
 FastClosures = "9aa1b823-49e4-5ca5-8b0f-3971ec8bab6a"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
-PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-Preferences = "21216c6a-2e73-6563-6e65-726566657250"
 RecursiveArrayTools = "731186ca-8d62-57ce-b412-fbd966d074cd"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
@@ -36,7 +33,6 @@ ConcreteStructs = "0.2.3"
 DiffEqBase = "6.167"
 DiffEqDevTools = "2.44"
 DifferentiationInterface = "0.6.42, 0.7"
-FastAlmostBandedMatrices = "0.1.4"
 FastClosures = "0.3.2"
 ForwardDiff = "0.10.38, 1"
 Hwloc = "3"
@@ -49,8 +45,6 @@ OrdinaryDiffEqTsit5 = "1"
 OrdinaryDiffEqVerner = "1"
 Pkg = "1.10.0"
 PreallocationTools = "0.4.24"
-PrecompileTools = "1.2"
-Preferences = "1.4"
 Random = "1.10"
 ReTestItems = "1.23.1"
 RecursiveArrayTools = "3.27.0"
diff --git a/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl b/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
index 20ad0b2a5..d8d7faed9 100644
--- a/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/BoundaryValueDiffEqShooting.jl
@@ -24,13 +24,9 @@ using ConcreteStructs: @concrete
 using DiffEqBase: DiffEqBase, solve
 using DifferentiationInterface: DifferentiationInterface, Constant, prepare_jacobian,
                                 overloaded_input_type
-using FastAlmostBandedMatrices: AlmostBandedMatrix, fillpart, exclusive_bandpart,
-                                finish_part_setindex!
 using FastClosures: @closure
 using ForwardDiff: ForwardDiff, pickchunksize
 using LinearAlgebra
-using PrecompileTools: @compile_workload, @setup_workload
-using Preferences: Preferences
 using Reexport: @reexport
 using RecursiveArrayTools: ArrayPartition, DiffEqArray, VectorOfArray
 using SciMLBase: SciMLBase, AbstractDiffEqInterpolation, StandardBVProblem, __solve,

From c310ae6cfcbf47dc014fd4b416120a4542f95227 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sun, 27 Jul 2025 23:53:15 +0800
Subject: [PATCH 33/50] Improve docs for interface

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

diff --git a/docs/src/basics/solve.md b/docs/src/basics/solve.md
index 281fc94b8..a90e3b11c 100644
--- a/docs/src/basics/solve.md
+++ b/docs/src/basics/solve.md
@@ -5,7 +5,7 @@
   - `controller`: Error controller for collocation methods, default as `DefectControl()`, more controller options in [Error Control Adaptivity](@ref error_control).
   - `defect_threshold`: Monitor of the size of defect norm. Defaults to `0.1`.
   - `odesolve_kwargs`: OrdinaryDiffEq.jl solvers kwargs for passing to ODE solving in shooting methods. For more information, see the documentation for OrdinaryDiffEq: [Common Solver Options](https://docs.sciml.ai/DiffEqDocs/latest/basics/common_solver_opts/).
-  - `nlsolve_kwargs`: NonlinearSolve.jl solvers kwargs for passing to nonlinear solving in collocation methods and shooting methods. For more information, see the documentation for NonlinearSolve: [Commom Solver Options](https://docs.sciml.ai/NonlinearSolve/stable/basics/solve/). The default internal nonlinear solver is [customized polyalgorithm](https://github.com/SciML/BoundaryValueDiffEq.jl/blob/master/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl) and the default absolute tolerance of nonlinear solving in collocation methods is `1e-6`.
-  - `optimize_kwargs`: Optimization.jl solvers kwargs for passing to optimizing in collocation methods and shooting methods. For more information, see the documentation for Optimization: [Commom Solver Options](https://docs.sciml.ai/Optimization/stable/API/solve/). The internal optimization solver should be specified and the default absolute tolerance of optimization problem solving in collocation methods is `1e-6`.
+  - `nlsolve_kwargs`: NonlinearSolve.jl solvers kwargs for passing to nonlinear solving in collocation methods and shooting methods. For more information, see the documentation for NonlinearSolve: [Commom Solver Options](https://docs.sciml.ai/NonlinearSolve/stable/basics/solve/). The default internal nonlinear solver is [customized polyalgorithm](https://github.com/SciML/BoundaryValueDiffEq.jl/blob/master/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl) and the default absolute tolerance of nonlinear solving in collocation and shooting methods is `1e-6`.
+  - `optimize_kwargs`: Optimization.jl solvers kwargs for passing to optimizing in collocation methods and shooting methods. For more information, see the documentation for Optimization: [Commom Solver Options](https://docs.sciml.ai/Optimization/stable/API/solve/). The internal optimization solver should be specified and the default absolute tolerance of optimization problem solving in collocation and shooting methods is `1e-6`.
   - `verbose`:  Toggles whether warnings are thrown when the solver exits early. Defaults to `true`.
   - `ensemblealg`: Whether `MultipleShooting` uses multithreading, default as `EnsembleThreads()`. For more information, see the documentation for OrdinaryDiffEq: [EnsembleAlgorithms](https://docs.sciml.ai/DiffEqDocs/latest/features/ensemble/#EnsembleAlgorithms).

From 678e41319bb8e387a6ea852dff09acabb90f9fbb Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Wed, 30 Jul 2025 21:23:24 +0800
Subject: [PATCH 34/50] format

---
 lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl | 2 +-
 lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl   | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
index d01ebcfbc..d0b1040b5 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
@@ -331,7 +331,7 @@ end
         end
 
         @testset "Derivative Interpolation tests for RadauIIa$stage" for stage in
-                                                                        (2, 3, 5, 7)
+                                                                         (2, 3, 5, 7)
             @time sol = solve(prob_bvp_linear, radau_solver(Val(stage)); dt = 0.001)
             sol_analytic = prob_bvp_linear_analytic(nothing, λ, 0.04)
             dsol_analytic = prob_bvp_linear_analytic_derivative(nothing, λ, 0.04)
diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
index 59451495e..1fb1d69f4 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
@@ -352,7 +352,7 @@ end
         end
 
         @testset "Derivative Interpolation tests for RadauIIa$stage" for stage in
-                                                                        (2, 3, 5, 7)
+                                                                         (2, 3, 5, 7)
             @time sol = solve(prob_bvp_linear, radau_solver(Val(stage)); dt = 0.001)
             sol_analytic = prob_bvp_linear_analytic(nothing, λ, 0.04)
             dsol_analytic = prob_bvp_linear_analytic_derivative(nothing, λ, 0.04)

From e4486ab318354d1550099f46cd7b7d8828549084 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 2 Aug 2025 01:21:11 +0800
Subject: [PATCH 35/50] Fix merge conflicts

---
 lib/BoundaryValueDiffEqFIRK/src/firk.jl | 4 ++--
 lib/BoundaryValueDiffEqMIRK/src/mirk.jl | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/src/firk.jl b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
index 45fb9a452..50a5845a7 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/firk.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
@@ -405,8 +405,8 @@ function SciMLBase.solve!(cache::FIRKCacheNested{
 end
 
 function __perform_firk_iteration(cache::Union{FIRKCacheExpand, FIRKCacheNested}, abstol, adaptive::Bool)
-    nlprob = __construct_nlproblem(cache, vec(cache.y₀), copy(cache.y₀))
-    nlsolve_alg = __concrete_nonlinearsolve_algorithm(nlprob, cache.alg.nlsolve)
+    nlprob = __construct_problem(cache, vec(cache.y₀), copy(cache.y₀))
+    nlsolve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve)
     sol_nlprob = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
index 3e7d1b78d..def010411 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
@@ -238,7 +238,7 @@ function __perform_mirk_iteration(cache::MIRKCache, abstol, adaptive::Bool, cont
 end
 
 # Constructing the Nonlinear Problem
-function __construct_nlproblem(cache::MIRKCache{iip}, y::AbstractVector, y₀::AbstractVectorOfArray) where {iip}
+function __construct_problem(cache::MIRKCache{iip}, y::AbstractVector, y₀::AbstractVectorOfArray) where {iip}
     pt = cache.problem_type
     (; jac_alg) = cache.alg
 

From 5c00fe9dab8dfdcbb83c8bf62985613679b8acfe Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Mon, 4 Aug 2025 23:55:05 +0800
Subject: [PATCH 36/50] Update some docs

---
 docs/src/basics/solve.md                 | 2 +-
 docs/src/tutorials/unknown_parameters.md | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/basics/solve.md b/docs/src/basics/solve.md
index a90e3b11c..e65b6b7fb 100644
--- a/docs/src/basics/solve.md
+++ b/docs/src/basics/solve.md
@@ -6,6 +6,6 @@
   - `defect_threshold`: Monitor of the size of defect norm. Defaults to `0.1`.
   - `odesolve_kwargs`: OrdinaryDiffEq.jl solvers kwargs for passing to ODE solving in shooting methods. For more information, see the documentation for OrdinaryDiffEq: [Common Solver Options](https://docs.sciml.ai/DiffEqDocs/latest/basics/common_solver_opts/).
   - `nlsolve_kwargs`: NonlinearSolve.jl solvers kwargs for passing to nonlinear solving in collocation methods and shooting methods. For more information, see the documentation for NonlinearSolve: [Commom Solver Options](https://docs.sciml.ai/NonlinearSolve/stable/basics/solve/). The default internal nonlinear solver is [customized polyalgorithm](https://github.com/SciML/BoundaryValueDiffEq.jl/blob/master/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl) and the default absolute tolerance of nonlinear solving in collocation and shooting methods is `1e-6`.
-  - `optimize_kwargs`: Optimization.jl solvers kwargs for passing to optimizing in collocation methods and shooting methods. For more information, see the documentation for Optimization: [Commom Solver Options](https://docs.sciml.ai/Optimization/stable/API/solve/). The internal optimization solver should be specified and the default absolute tolerance of optimization problem solving in collocation and shooting methods is `1e-6`.
+  - `optimize_kwargs`: Optimization.jl solvers kwargs for passing to optimization problem solving in collocation methods and shooting methods. For more information, see the documentation for Optimization: [Commom Solver Options](https://docs.sciml.ai/Optimization/stable/API/solve/). The internal optimization solver should be specified and the default absolute tolerance of optimization problem solving in collocation and shooting methods is `1e-6`.
   - `verbose`:  Toggles whether warnings are thrown when the solver exits early. Defaults to `true`.
   - `ensemblealg`: Whether `MultipleShooting` uses multithreading, default as `EnsembleThreads()`. For more information, see the documentation for OrdinaryDiffEq: [EnsembleAlgorithms](https://docs.sciml.ai/DiffEqDocs/latest/features/ensemble/#EnsembleAlgorithms).
diff --git a/docs/src/tutorials/unknown_parameters.md b/docs/src/tutorials/unknown_parameters.md
index b064d5ba5..5bb315300 100644
--- a/docs/src/tutorials/unknown_parameters.md
+++ b/docs/src/tutorials/unknown_parameters.md
@@ -54,7 +54,7 @@ sol = solve(bvp, MIRK4(), dt = 0.05)
 plot(sol)
 ```
 
-after solving the boundary value problem, the estimated unknown parameters can be accessed with
+after solving the boundary value problem, the estimated unknown parameters can be accessed in the solution
 
 ```@example unknown
 sol.prob.p

From df3990a221e74c86a8139ab6d08db89ceda97d64 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Fri, 8 Aug 2025 23:46:56 +0800
Subject: [PATCH 37/50] Don't set alias_u0 for now

---
 lib/BoundaryValueDiffEqCore/src/utils.jl | 14 ++++----------
 1 file changed, 4 insertions(+), 10 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index dbc7b7cb7..6137a8b53 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -615,17 +615,11 @@ function __split_kwargs(; abstol, adaptive, controller, kwargs...)
     return ((abstol, adaptive, controller), (; abstol, adaptive, kwargs...))
 end
 
-@inline __concrete_kwargs(nlsolve,
-    ::Nothing,
-    nlsolve_kwargs,
-    optimize_kwargs) = (
-    nlsolve_kwargs..., alias = SciMLBase.NonlinearAliasSpecifier(alias_u0 = true))
+@inline __concrete_kwargs(nlsolve, ::Nothing, nlsolve_kwargs, optimize_kwargs) = (;
+    nlsolve_kwargs...)
 @inline __concrete_kwargs(::Nothing, optimize, nlsolve_kwargs, optimize_kwargs) = (;) # Doesn't support for now
-@inline __concrete_kwargs(::Nothing,
-    ::Nothing,
-    nlsolve_kwargs,
-    optimize_kwargs) = (
-    nlsolve_kwargs..., alias = SciMLBase.NonlinearAliasSpecifier(alias_u0 = true))
+@inline __concrete_kwargs(::Nothing, ::Nothing, nlsolve_kwargs, optimize_kwargs) = (;
+    nlsolve_kwargs...)
 
 ## Optimization solver related utils ##
 

From 7a45c7d22fa427a52e2edce404d68d25ffd5ea8b Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 9 Aug 2025 00:19:49 +0800
Subject: [PATCH 38/50] Fix typos

---
 docs/src/basics/error_control.md               | 2 +-
 docs/src/basics/solve.md                       | 4 ++--
 docs/src/solvers/mirk.md                       | 2 +-
 docs/src/solvers/mirkn.md                      | 2 +-
 docs/src/solvers/shooting.md                   | 2 +-
 docs/src/tutorials/solve_nlls_bvp.md           | 2 +-
 lib/BoundaryValueDiffEqCore/src/calc_errors.jl | 2 +-
 7 files changed, 8 insertions(+), 8 deletions(-)

diff --git a/docs/src/basics/error_control.md b/docs/src/basics/error_control.md
index 5d263ed1f..f41a602c2 100644
--- a/docs/src/basics/error_control.md
+++ b/docs/src/basics/error_control.md
@@ -4,7 +4,7 @@ Adaptivity helps ensure the quality of the our numerical solution, and when our
 
 When comes to solving ill-conditioned BVP, for example the singular perturbation problem where the small parameters become extremely small leading to the layers phonemona, the error control adaptivity becomes even more critical, because the minor perturbations can lead to large deviation in the solution. In such cases, adaptivity automatically figure out where to use refined mesh and where to use coarse mesh to achieve the balance of computational efficiency and accuracy.
 
-BoundaryValuDiffEq.jl support error control adaptivity for collocation methods, and the adaptivity is default as defect control adaptivity when using adaptive collocation solvers:
+BoundaryValueDiffEq.jl support error control adaptivity for collocation methods, and the adaptivity is default as defect control adaptivity when using adaptive collocation solvers:
 
 ```julia
 sol = solve(prob, MIRK4(), dt = 0.01, adaptive = true)
diff --git a/docs/src/basics/solve.md b/docs/src/basics/solve.md
index e65b6b7fb..4d48df59a 100644
--- a/docs/src/basics/solve.md
+++ b/docs/src/basics/solve.md
@@ -5,7 +5,7 @@
   - `controller`: Error controller for collocation methods, default as `DefectControl()`, more controller options in [Error Control Adaptivity](@ref error_control).
   - `defect_threshold`: Monitor of the size of defect norm. Defaults to `0.1`.
   - `odesolve_kwargs`: OrdinaryDiffEq.jl solvers kwargs for passing to ODE solving in shooting methods. For more information, see the documentation for OrdinaryDiffEq: [Common Solver Options](https://docs.sciml.ai/DiffEqDocs/latest/basics/common_solver_opts/).
-  - `nlsolve_kwargs`: NonlinearSolve.jl solvers kwargs for passing to nonlinear solving in collocation methods and shooting methods. For more information, see the documentation for NonlinearSolve: [Commom Solver Options](https://docs.sciml.ai/NonlinearSolve/stable/basics/solve/). The default internal nonlinear solver is [customized polyalgorithm](https://github.com/SciML/BoundaryValueDiffEq.jl/blob/master/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl) and the default absolute tolerance of nonlinear solving in collocation and shooting methods is `1e-6`.
-  - `optimize_kwargs`: Optimization.jl solvers kwargs for passing to optimization problem solving in collocation methods and shooting methods. For more information, see the documentation for Optimization: [Commom Solver Options](https://docs.sciml.ai/Optimization/stable/API/solve/). The internal optimization solver should be specified and the default absolute tolerance of optimization problem solving in collocation and shooting methods is `1e-6`.
+  - `nlsolve_kwargs`: NonlinearSolve.jl solvers kwargs for passing to nonlinear solving in collocation methods and shooting methods. For more information, see the documentation for NonlinearSolve: [Common Solver Options](https://docs.sciml.ai/NonlinearSolve/stable/basics/solve/). The default internal nonlinear solver is [customized polyalgorithm](https://github.com/SciML/BoundaryValueDiffEq.jl/blob/master/lib/BoundaryValueDiffEqCore/src/default_nlsolve.jl) and the default absolute tolerance of nonlinear solving in collocation and shooting methods is `1e-6`.
+  - `optimize_kwargs`: Optimization.jl solvers kwargs for passing to optimization problem solving in collocation methods and shooting methods. For more information, see the documentation for Optimization: [Common Solver Options](https://docs.sciml.ai/Optimization/stable/API/solve/). The internal optimization solver should be specified and the default absolute tolerance of optimization problem solving in collocation and shooting methods is `1e-6`.
   - `verbose`:  Toggles whether warnings are thrown when the solver exits early. Defaults to `true`.
   - `ensemblealg`: Whether `MultipleShooting` uses multithreading, default as `EnsembleThreads()`. For more information, see the documentation for OrdinaryDiffEq: [EnsembleAlgorithms](https://docs.sciml.ai/DiffEqDocs/latest/features/ensemble/#EnsembleAlgorithms).
diff --git a/docs/src/solvers/mirk.md b/docs/src/solvers/mirk.md
index 823027c61..7b4d7cabf 100644
--- a/docs/src/solvers/mirk.md
+++ b/docs/src/solvers/mirk.md
@@ -1,6 +1,6 @@
 # [BoundaryValueDiffEqMIRK](@id mirk)
 
-Monotonic Implicit Runge Kutta(MIRK) Methods. To only use the MIRK methods form BoundaryVaueDiffEq.jl, you need to install them use the Julia package manager:
+Monotonic Implicit Runge Kutta(MIRK) Methods. To only use the MIRK methods form BoundaryValueDiffEq.jl, you need to install them use the Julia package manager:
 
 ```julia
 using Pkg
diff --git a/docs/src/solvers/mirkn.md b/docs/src/solvers/mirkn.md
index a646d4780..e6266f34f 100644
--- a/docs/src/solvers/mirkn.md
+++ b/docs/src/solvers/mirkn.md
@@ -1,6 +1,6 @@
 # [BoundaryValueDiffEqMIRKN](@id mirkn)
 
-Monotonic Implicit Runge Kutta Nyström(MIRKN) Methods. To only use the MIRKN methods form BoundaryVaueDiffEq.jl, you need to install them use the Julia package manager:
+Monotonic Implicit Runge Kutta Nyström(MIRKN) Methods. To only use the MIRKN methods form BoundaryValueDiffEq.jl, you need to install them use the Julia package manager:
 
 ```julia
 using Pkg
diff --git a/docs/src/solvers/shooting.md b/docs/src/solvers/shooting.md
index 0da879eed..5200e751e 100644
--- a/docs/src/solvers/shooting.md
+++ b/docs/src/solvers/shooting.md
@@ -1,6 +1,6 @@
 # [BoundaryValueDiffEqShooting](@id shooting)
 
-Single shooting method and multiple shooting method. To only use the Shooting methods form BoundaryVaueDiffEq.jl, you need to install them use the Julia package manager:
+Single shooting method and multiple shooting method. To only use the Shooting methods form BoundaryValueDiffEq.jl, you need to install them use the Julia package manager:
 
 ```julia
 using Pkg
diff --git a/docs/src/tutorials/solve_nlls_bvp.md b/docs/src/tutorials/solve_nlls_bvp.md
index 23e463b39..fa3cf7a6c 100644
--- a/docs/src/tutorials/solve_nlls_bvp.md
+++ b/docs/src/tutorials/solve_nlls_bvp.md
@@ -39,7 +39,7 @@ sol = solve(prob, MIRK4(), dt = 0.01, abstol = 1e-3)
 plot(sol)
 ```
 
-Since this BVP imposes constraints only at the two endpoints, we can use `TwoPointBVProlem` to handle such cases.
+Since this BVP imposes constraints only at the two endpoints, we can use `TwoPointBVProblem` to handle such cases.
 
 ```@example nlls_overdetermined
 function f!(du, u, p, t)
diff --git a/lib/BoundaryValueDiffEqCore/src/calc_errors.jl b/lib/BoundaryValueDiffEqCore/src/calc_errors.jl
index 9c2397624..1271f3bb5 100644
--- a/lib/BoundaryValueDiffEqCore/src/calc_errors.jl
+++ b/lib/BoundaryValueDiffEqCore/src/calc_errors.jl
@@ -35,7 +35,7 @@ struct DefectControl{T} <: AbstractErrorControl
 end
 
 """
-    GlobalErrorControl(; method = HOErorControl())
+    GlobalErrorControl(; method = HOErrorControl())
 
 Global error controller, use high order global error estimation method `HOErrorControl` as default.
 """

From 3a51b8f750ea5a4f1b7b8e027371ef5077ff4673 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 9 Aug 2025 01:55:33 +0800
Subject: [PATCH 39/50] Fix construct issue with FIRK

---
 .../test/expanded/ensemble_tests.jl           |  9 ++--
 .../test/expanded/firk_basic_tests.jl         | 43 ++++++++++---------
 2 files changed, 29 insertions(+), 23 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
index a095ed9c4..68b72b692 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/ensemble_tests.jl
@@ -40,7 +40,8 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested); trajectories = 10, dt = 0.1)
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested);
+                trajectories = 10, dt = 0.1)
             @test sol.converged
         end
     end
@@ -51,7 +52,8 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested); trajectories = 10, dt = 0.1)
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested);
+                trajectories = 10, dt = 0.1)
             @test sol.converged
         end
     end
@@ -62,7 +64,8 @@
                 AutoSparse(AutoFiniteDiff()); bc_diffmode = AutoFiniteDiff(),
                 nonbc_diffmode = AutoSparse(AutoFiniteDiff()))]
         for jac_alg in jac_algs
-            sol = solve(ensemble_prob, solver(nlsolve, jac_alg; nested); trajectories = 10, dt = 0.1)
+            sol = solve(ensemble_prob, solver(; nlsolve, jac_alg, nested_nlsolve = nested);
+                trajectories = 10, dt = 0.1)
             @test sol.converged
         end
     end
diff --git a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
index c85eb9ad7..e6a421cdf 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/expanded/firk_basic_tests.jl
@@ -236,26 +236,29 @@ end
     nested = false
 
     # Using ForwardDiff might lead to Cache expansion warnings
-    @test_nowarn solve(bvp1, LobattoIIIa2(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa5(nl_solve, jac_alg; nested); dt = 0.005)
-
-    @test_nowarn solve(bvp1, LobattoIIIb2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, LobattoIIIb3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIb4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIb5(nl_solve, jac_alg; nested); dt = 0.005)
-
-    @test_nowarn solve(bvp1, LobattoIIIc2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, LobattoIIIc3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIc4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIc5(nl_solve, jac_alg; nested); dt = 0.005)
-
-    @test_nowarn solve(bvp1, RadauIIa1(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, RadauIIa2(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, RadauIIa3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, RadauIIa5(nl_solve, jac_alg; nested); dt = 0.05)
-    @test_nowarn solve(bvp1, RadauIIa7(nl_solve, jac_alg; nested); dt = 0.05)
+    @test_nowarn solve(bvp1, LobattoIIIa2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIa3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIa4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIa5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+
+    @test_nowarn solve(
+        bvp1, LobattoIIIb2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
+    @test_nowarn solve(bvp1, LobattoIIIb3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIb4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIb5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+
+    @test_nowarn solve(
+        bvp1, LobattoIIIc2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
+    @test_nowarn solve(bvp1, LobattoIIIc3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIc4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIc5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+
+    @test_nowarn solve(
+        bvp1, RadauIIa1(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
+    @test_nowarn solve(bvp1, RadauIIa2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, RadauIIa3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, RadauIIa5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.05)
+    @test_nowarn solve(bvp1, RadauIIa7(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.05)
 end
 
 @testitem "Interpolation" setup=[FIRKExpandedConvergenceTests] begin

From 2aaa820c71ff386962d1a02eda08f923625e333d Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 9 Aug 2025 11:51:59 +0800
Subject: [PATCH 40/50] Fix basic tests in nested FIRK

---
 .../test/nested/firk_basic_tests.jl           | 43 ++++++++++---------
 1 file changed, 23 insertions(+), 20 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
index 4a2d96489..a01acde4e 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
@@ -256,26 +256,29 @@ end
     nested = true
 
     # Using ForwardDiff might lead to Cache expansion warnings
-    @test_nowarn solve(bvp1, LobattoIIIa2(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa5(nl_solve, jac_alg; nested); dt = 0.005)
-
-    @test_nowarn solve(bvp1, LobattoIIIb2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, LobattoIIIb3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIb4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIb5(nl_solve, jac_alg; nested); dt = 0.005)
-
-    @test_nowarn solve(bvp1, LobattoIIIc2(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, LobattoIIIc3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIc4(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIc5(nl_solve, jac_alg; nested); dt = 0.005)
-
-    @test_nowarn solve(bvp1, RadauIIa1(nl_solve, jac_alg; nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, RadauIIa2(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, RadauIIa3(nl_solve, jac_alg; nested); dt = 0.005)
-    @test_nowarn solve(bvp1, RadauIIa5(nl_solve, jac_alg; nested); dt = 0.05)
-    @test_nowarn solve(bvp1, RadauIIa7(nl_solve, jac_alg; nested); dt = 0.05)
+    @test_nowarn solve(bvp1, LobattoIIIa2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIa3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIa4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIa5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+
+    @test_nowarn solve(
+        bvp1, LobattoIIIb2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
+    @test_nowarn solve(bvp1, LobattoIIIb3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIb4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIb5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+
+    @test_nowarn solve(
+        bvp1, LobattoIIIc2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
+    @test_nowarn solve(bvp1, LobattoIIIc3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIc4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIc5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+
+    @test_nowarn solve(
+        bvp1, RadauIIa1(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
+    @test_nowarn solve(bvp1, RadauIIa2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, RadauIIa3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, RadauIIa5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.05)
+    @test_nowarn solve(bvp1, RadauIIa7(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.05)
 end
 
 @testitem "Interpolation" setup=[FIRKNestedConvergenceTests] begin

From 1288c889e7144d4619de87408559690e76b91757 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 9 Aug 2025 13:24:02 +0800
Subject: [PATCH 41/50] Fix keyword specification in nested FIRK

---
 .../test/nested/firk_basic_tests.jl           | 36 +++++++++----------
 1 file changed, 18 insertions(+), 18 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
index a01acde4e..118fcf47e 100644
--- a/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqFIRK/test/nested/firk_basic_tests.jl
@@ -252,33 +252,33 @@ end
 
     jac_alg = BVPJacobianAlgorithm(;
         bc_diffmode = AutoFiniteDiff(), nonbc_diffmode = AutoSparse(AutoFiniteDiff()))
-    nl_solve = NewtonRaphson()
+    nlsolve = NewtonRaphson()
     nested = true
 
     # Using ForwardDiff might lead to Cache expansion warnings
-    @test_nowarn solve(bvp1, LobattoIIIa2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIa5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIa2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIa3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIa4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIa5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
     @test_nowarn solve(
-        bvp1, LobattoIIIb2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, LobattoIIIb3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIb4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIb5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIb2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
+    @test_nowarn solve(bvp1, LobattoIIIb3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIb4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIb5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
     @test_nowarn solve(
-        bvp1, LobattoIIIc2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, LobattoIIIc3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIc4(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
-    @test_nowarn solve(bvp1, LobattoIIIc5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+        bvp1, LobattoIIIc2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
+    @test_nowarn solve(bvp1, LobattoIIIc3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIc4(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, LobattoIIIc5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
 
     @test_nowarn solve(
-        bvp1, RadauIIa1(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
-    @test_nowarn solve(bvp1, RadauIIa2(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
-    @test_nowarn solve(bvp1, RadauIIa3(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.005)
-    @test_nowarn solve(bvp1, RadauIIa5(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.05)
-    @test_nowarn solve(bvp1, RadauIIa7(; nl_solve, jac_alg, nested_nlsolve = nested); dt = 0.05)
+        bvp1, RadauIIa1(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005, adaptive = false)
+    @test_nowarn solve(bvp1, RadauIIa2(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, RadauIIa3(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.005)
+    @test_nowarn solve(bvp1, RadauIIa5(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.05)
+    @test_nowarn solve(bvp1, RadauIIa7(; nlsolve, jac_alg, nested_nlsolve = nested); dt = 0.05)
 end
 
 @testitem "Interpolation" setup=[FIRKNestedConvergenceTests] begin

From 9118781b370ba321c7cf253e9c7932f1180670d6 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 9 Aug 2025 14:58:01 +0800
Subject: [PATCH 42/50] Fix internal solve

---
 lib/BoundaryValueDiffEqFIRK/src/firk.jl   | 6 ++++--
 lib/BoundaryValueDiffEqMIRK/src/mirk.jl   | 2 +-
 lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl | 2 +-
 3 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/lib/BoundaryValueDiffEqFIRK/src/firk.jl b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
index 50a5845a7..a6d7abb57 100644
--- a/lib/BoundaryValueDiffEqFIRK/src/firk.jl
+++ b/lib/BoundaryValueDiffEqFIRK/src/firk.jl
@@ -406,8 +406,10 @@ end
 
 function __perform_firk_iteration(cache::Union{FIRKCacheExpand, FIRKCacheNested}, abstol, adaptive::Bool)
     nlprob = __construct_problem(cache, vec(cache.y₀), copy(cache.y₀))
-    nlsolve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve)
-    sol_nlprob = __solve(nlprob, nlsolve_alg; cache.nlsolve_kwargs..., alias_u0 = true)
+    solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize)
+    kwargs = __concrete_kwargs(
+        cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs)
+    sol_nlprob = __solve(nlprob, solve_alg; kwargs...)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
     defect_norm = 2 * abstol
diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
index def010411..4a9851a03 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
@@ -192,7 +192,7 @@ function __perform_mirk_iteration(cache::MIRKCache, abstol, adaptive::Bool, cont
     solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize)
     kwargs = __concrete_kwargs(
         cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs)
-    sol_nlprob = solve(nlprob, solve_alg, kwargs...)
+    sol_nlprob = solve(nlprob, solve_alg; kwargs...)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
     error_norm = 2 * abstol
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
index c21107976..a0bacc04d 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
@@ -119,7 +119,7 @@ function __perform_mirkn_iteration(cache::MIRKNCache)
     solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize)
     kwargs = __concrete_kwargs(
         cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs)
-    sol_nlprob = solve(nlprob, solve_alg, kwargs...)
+    sol_nlprob = solve(nlprob, solve_alg; kwargs...)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
     return sol_nlprob, sol_nlprob.retcode

From c477425aa5678f1c868151421fd01053406fccc3 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Wed, 13 Aug 2025 23:54:15 +0800
Subject: [PATCH 43/50] Use correct solve

---
 lib/BoundaryValueDiffEqAscher/src/ascher.jl            |  2 +-
 lib/BoundaryValueDiffEqMIRK/src/mirk.jl                |  2 +-
 lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl              |  2 +-
 .../src/multiple_shooting.jl                           |  4 ++--
 lib/BoundaryValueDiffEqShooting/src/single_shooting.jl |  2 +-
 test/misc/default_solvers.jl                           | 10 +++++-----
 test/misc/qa_tests.jl                                  |  6 ++----
 7 files changed, 13 insertions(+), 15 deletions(-)

diff --git a/lib/BoundaryValueDiffEqAscher/src/ascher.jl b/lib/BoundaryValueDiffEqAscher/src/ascher.jl
index e385b6947..41804f463 100644
--- a/lib/BoundaryValueDiffEqAscher/src/ascher.jl
+++ b/lib/BoundaryValueDiffEqAscher/src/ascher.jl
@@ -181,7 +181,7 @@ function __perform_ascher_iteration(cache::AscherCache{iip, T}, abstol, adaptive
     solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize)
     kwargs = __concrete_kwargs(
         cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs)
-    nlsol = solve(nlprob, solve_alg; kwargs...)
+    nlsol = __solve(nlprob, solve_alg; kwargs...)
     error_norm = 2 * abstol
     info = nlsol.retcode
 
diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
index 4a9851a03..3a45d2474 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
@@ -192,7 +192,7 @@ function __perform_mirk_iteration(cache::MIRKCache, abstol, adaptive::Bool, cont
     solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize)
     kwargs = __concrete_kwargs(
         cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs)
-    sol_nlprob = solve(nlprob, solve_alg; kwargs...)
+    sol_nlprob = __solve(nlprob, solve_alg; kwargs...)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
     error_norm = 2 * abstol
diff --git a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
index a0bacc04d..05bdb2cd2 100644
--- a/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
+++ b/lib/BoundaryValueDiffEqMIRKN/src/mirkn.jl
@@ -119,7 +119,7 @@ function __perform_mirkn_iteration(cache::MIRKNCache)
     solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize)
     kwargs = __concrete_kwargs(
         cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs)
-    sol_nlprob = solve(nlprob, solve_alg; kwargs...)
+    sol_nlprob = __solve(nlprob, solve_alg; kwargs...)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
     return sol_nlprob, sol_nlprob.retcode
diff --git a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
index 8b15099e8..a384f4eda 100644
--- a/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/multiple_shooting.jl
@@ -148,7 +148,7 @@ function __solve_nlproblem!(
         u_at_nodes, prob.p, M, length(nodes), nothing)
 
     nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize)
-    solve(nlprob, nlsolve_alg; kwargs...)
+    __solve(nlprob, nlsolve_alg; kwargs...)
 
     return nothing
 end
@@ -223,7 +223,7 @@ function __solve_nlproblem!(::StandardBVProblem, alg::MultipleShooting, bcresid_
         prob, alg, loss_fn, jac_fn, jac_prototype, resid_prototype,
         u_at_nodes, prob.p, M, length(nodes), nothing)
     nlsolve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize)
-    solve(nlprob, nlsolve_alg; kwargs...)
+    __solve(nlprob, nlsolve_alg; kwargs...)
 
     return nothing
 end
diff --git a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
index 231260977..c8e9a567d 100644
--- a/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
+++ b/lib/BoundaryValueDiffEqShooting/src/single_shooting.jl
@@ -80,7 +80,7 @@ function SciMLBase.__solve(prob::BVProblem, alg_::Shooting; abstol = 1e-6,
     solve_alg = __concrete_solve_algorithm(nlprob, alg.nlsolve, alg.optimize)
     kwargs = __concrete_kwargs(alg.nlsolve, alg.optimize, nlsolve_kwargs, optimize_kwargs)
     #TODO: add verbose kwarg
-    nlsol = solve(nlprob, solve_alg; kwargs...)
+    nlsol = __solve(nlprob, solve_alg; kwargs...)
 
     # There is no way to reinit with the same cache with different cache. But not saving
     # the internal values gives a significant speedup. So we just create a new cache
diff --git a/test/misc/default_solvers.jl b/test/misc/default_solvers.jl
index 91cf9543a..7c936b88c 100644
--- a/test/misc/default_solvers.jl
+++ b/test/misc/default_solvers.jl
@@ -1,24 +1,24 @@
 @testitem "Default Solvers" begin
     using BoundaryValueDiffEq, Test
-    
+
     function f(du, u, p, t)
         (x, v) = u
         du[1] = v
         du[2] = -x
     end
-    
+
     function bc!(resid, sol, p, t)
         resid[1] = sol[1][1]
         resid[2] = sol[end][1] - 1
     end
-    
+
     tspan = (0.0, 100.0)
     u0 = [0.0, 1.0]
     bvp = BVProblem(f, bc!, u0, tspan)
     resid_f = Array{Float64}(undef, 2)
     sol = solve(bvp, Shooting(Tsit5()))
     sol2 = solve(bvp)
-    
+
     @test sol2.alg == Tsit5()
     @test all(sol.u .== sol2.u)
-end
\ No newline at end of file
+end
diff --git a/test/misc/qa_tests.jl b/test/misc/qa_tests.jl
index e406bc0ca..78dc32757 100644
--- a/test/misc/qa_tests.jl
+++ b/test/misc/qa_tests.jl
@@ -1,10 +1,8 @@
 @testitem "Quality Assurance" begin
     using Aqua
 
-    Aqua.test_all(BoundaryValueDiffEq; 
-        ambiguities = false,
-        piracies = (broken = false, 
-                   treat_as_own = [SciMLBase.BVProblem]))
+    Aqua.test_all(BoundaryValueDiffEq; ambiguities = false,
+        piracies = (broken = false, treat_as_own = [SciMLBase.BVProblem]))
 end
 
 @testitem "JET Package Test" begin

From 16c129f1f69ece05e3222c2c9584ceb738bf0ea3 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Thu, 14 Aug 2025 17:47:20 +0800
Subject: [PATCH 44/50] Finally done

---
 lib/BoundaryValueDiffEqCore/src/utils.jl                   | 4 ++++
 lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl | 2 +-
 lib/BoundaryValueDiffEqMIRK/src/mirk.jl                    | 2 +-
 3 files changed, 6 insertions(+), 2 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 6137a8b53..3743f9ce5 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -711,3 +711,7 @@ function __construct_internal_problem(
             prob, optf, y, p; lcons = lcons, ucons = ucons)
     end
 end
+
+# Some optimization algorithms (solvers from interfacing packages) don't support the __solve(prob) interface
+@inline __internal_solve(prob::SciMLBase.NonlinearProblem, alg; kwargs...) = __solve(prob, alg; kwargs...)
+@inline __internal_solve(prob::SciMLBase.OptimizationProblem, alg; kwargs...) = solve(prob, alg; kwargs...)
diff --git a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
index b14341097..9801f9b31 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/BoundaryValueDiffEqMIRK.jl
@@ -23,7 +23,7 @@ using BoundaryValueDiffEqCore: AbstractBoundaryValueDiffEqAlgorithm,
                                __use_both_error_control, __default_coloring_algorithm,
                                DiffCacheNeeded, NoDiffCacheNeeded, __split_kwargs,
                                __concrete_kwargs, __FastShortcutNonlinearPolyalg,
-                               __construct_internal_problem
+                               __construct_internal_problem, __internal_solve
 
 using ConcreteStructs: @concrete
 using DiffEqBase: DiffEqBase
diff --git a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
index 3a45d2474..98d89646c 100644
--- a/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
+++ b/lib/BoundaryValueDiffEqMIRK/src/mirk.jl
@@ -192,7 +192,7 @@ function __perform_mirk_iteration(cache::MIRKCache, abstol, adaptive::Bool, cont
     solve_alg = __concrete_solve_algorithm(nlprob, cache.alg.nlsolve, cache.alg.optimize)
     kwargs = __concrete_kwargs(
         cache.alg.nlsolve, cache.alg.optimize, cache.nlsolve_kwargs, cache.optimize_kwargs)
-    sol_nlprob = __solve(nlprob, solve_alg; kwargs...)
+    sol_nlprob = __internal_solve(nlprob, solve_alg; kwargs...)
     recursive_unflatten!(cache.y₀, sol_nlprob.u)
 
     error_norm = 2 * abstol

From c935e0fc5123ae934aad47fac49f26eac1fa0437 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Thu, 14 Aug 2025 19:06:54 +0800
Subject: [PATCH 45/50] Fix NLLS internal solve dispatch

---
 lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl | 6 ++++++
 lib/BoundaryValueDiffEqCore/src/utils.jl                  | 4 ----
 2 files changed, 6 insertions(+), 4 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl b/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl
index c3de4a03e..800b6e952 100644
--- a/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl
+++ b/lib/BoundaryValueDiffEqCore/src/default_internal_solve.jl
@@ -74,3 +74,9 @@ end
         return __FastShortcutBVPCompatibleNonlinearPolyalg(eltype(prob.u0))
     end
 end
+
+# Some optimization algorithms (solvers from interfacing packages) don't support the __solve(prob) interface
+@inline __internal_solve(
+    prob::Union{SciMLBase.NonlinearProblem, SciMLBase.NonlinearLeastSquaresProblem},
+    alg; kwargs...) = __solve(prob, alg; kwargs...)
+@inline __internal_solve(prob::SciMLBase.OptimizationProblem, alg; kwargs...) = solve(prob, alg; kwargs...)
diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 3743f9ce5..6137a8b53 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -711,7 +711,3 @@ function __construct_internal_problem(
             prob, optf, y, p; lcons = lcons, ucons = ucons)
     end
 end
-
-# Some optimization algorithms (solvers from interfacing packages) don't support the __solve(prob) interface
-@inline __internal_solve(prob::SciMLBase.NonlinearProblem, alg; kwargs...) = __solve(prob, alg; kwargs...)
-@inline __internal_solve(prob::SciMLBase.OptimizationProblem, alg; kwargs...) = solve(prob, alg; kwargs...)

From f6bd8420d388ffffa35a104ee5f3b31a1186d19e Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Thu, 14 Aug 2025 23:44:00 +0800
Subject: [PATCH 46/50] Docs should all be fine

---
 docs/make.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/make.jl b/docs/make.jl
index f0ca989e8..07383973e 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -27,7 +27,7 @@ makedocs(; sitename = "BoundaryValueDiffEq.jl",
     clean = true,
     doctest = false,
     checkdocs = :exports,
-    warnonly = [:missing_docs],
+    warnonly = [:missing_docs, :cross_references],
     plugins = [bib, interlinks],
     format = Documenter.HTML(assets = ["assets/favicon.ico"],
         canonical = "https://docs.sciml.ai/BoundaryValueDiffEq/stable/"),

From 3bfe1ac35dc39f1b2e0b75ea6ded4d6d6898048a Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Fri, 15 Aug 2025 23:39:27 +0800
Subject: [PATCH 47/50] Better support for inequality constraints

---
 README.md                                | 15 +++++++++++++
 lib/BoundaryValueDiffEqCore/src/utils.jl | 27 +++++++++++++++++-------
 2 files changed, 34 insertions(+), 8 deletions(-)

diff --git a/README.md b/README.md
index 6c015b4c9..f61d47fe9 100644
--- a/README.md
+++ b/README.md
@@ -37,6 +37,21 @@ prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan)
 sol = solve(prob, MIRK4(), dt = 0.05)
 ```
 
+```julia
+tspan = (0.0, 1.0)
+function lotka_volterra!(du, u, p, t)
+    du[1] = 1.5 * u[1] - 1.0 * u[1] * u[2]
+    du[2] = -3.0 * u[2] + 1.0 * u[1] * u[2]
+end
+function bc!(residual, u, p, t)
+    residual[1] = u[1] - 1.0
+    residual[2] = u[2] - 1.0
+end
+prob = BVProblem(
+    lotka_volterra!, bc!, [4.0, 2.0], tspan, lcons = [0.0, 0.0], ucons = [Inf, Inf])
+sol = solve(prob, MIRK4(; optimizer = Ipopt.Optimizer()), dt = 0.05)
+```
+
 ## Available Solvers
 
 For the list of available solvers, please refer to the [DifferentialEquations.jl BVP Solvers page](https://docs.sciml.ai/DiffEqDocs/stable/solvers/bvp_solve/). For options for the `solve` command, see the [common solver options page](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/).
diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 6137a8b53..338dd94e6 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -627,6 +627,21 @@ end
 @inline __default_cost(f) = f
 @inline __default_cost(fun::BVPFunction) = __default_cost(fun.cost)
 
+@inline function __extract_lcons_ucons(prob::AbstractBVProblem, ::Type{T}, M, N) where {T}
+    inequality_length = length(prob.f.lcons)
+    lcons = if isnothing(prob.f.lcons)
+        zeros(T, N*M)
+    else
+        vcat(prob.f.lcons, zeros(T, N*M - inequality_length))
+    end
+    ucons = if isnothing(prob.f.ucons)
+        zeros(T, N*M)
+    else
+        vcat(prob.f.ucons, zeros(T, N*M - inequality_length))
+    end
+    return lcons, ucons
+end
+
 """
     __construct_internal_problem
 
@@ -646,8 +661,7 @@ function __construct_internal_problem(prob::AbstractBVProblem, alg, loss, jac,
         optf = OptimizationFunction{true}(__default_cost(prob.f), AutoFiniteDiff(), # Need to investigate the ForwardDiff dual problem
             cons = loss,
             cons_j = jac, cons_jac_prototype = jac_prototype)
-        lcons = zeros(T, N*M)
-        ucons = zeros(T, N*M)
+        lcons, ucons = __extract_lcons_ucons(prob, T, M, N)
         return __internal_optimization_problem(
             prob, optf, y, p; lcons = lcons, ucons = ucons)
     end
@@ -665,8 +679,7 @@ function __construct_internal_problem(prob::TwoPointBVProblem, alg, loss, jac,
         optf = OptimizationFunction{true}(
             __default_cost(prob.f), get_dense_ad(alg.jac_alg.diffmode),
             cons = loss, cons_j = jac, cons_jac_prototype = jac_prototype)
-        lcons = zeros(T, N*M)
-        ucons = zeros(T, N*M)
+        lcons, ucons = __extract_lcons_ucons(prob, T, M, N)
 
         return __internal_optimization_problem(
             prob, optf, y, p; lcons = lcons, ucons = ucons)
@@ -684,8 +697,7 @@ function __construct_internal_problem(prob, alg, loss, jac, jac_prototype,
         optf = OptimizationFunction{true}(
             __default_cost(prob.f), get_dense_ad(alg.jac_alg.diffmode),
             cons = loss, cons_j = jac, cons_jac_prototype = jac_prototype)
-        lcons = zeros(T, N*M)
-        ucons = zeros(T, N*M)
+        lcons, ucons = __extract_lcons_ucons(prob, T, M, N)
 
         return __internal_optimization_problem(
             prob, optf, y, p; lcons = lcons, ucons = ucons)
@@ -704,8 +716,7 @@ function __construct_internal_problem(
         optf = OptimizationFunction{true}(
             __default_cost(prob.f), get_dense_ad(alg.jac_alg.nonbc_diffmode),
             cons = loss, cons_j = jac, cons_jac_prototype = jac_prototype)
-        lcons = zeros(T, N*M)
-        ucons = zeros(T, N*M)
+        lcons, ucons = __extract_lcons_ucons(prob, T, M, N)
 
         return __internal_optimization_problem(
             prob, optf, y, p; lcons = lcons, ucons = ucons)

From f38194d786eca4abc255c3a240610efaccafade5 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 16 Aug 2025 01:12:55 +0800
Subject: [PATCH 48/50] Better inequality specification

---
 lib/BoundaryValueDiffEqCore/src/utils.jl | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 338dd94e6..7969ff83e 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -628,16 +628,17 @@ end
 @inline __default_cost(fun::BVPFunction) = __default_cost(fun.cost)
 
 @inline function __extract_lcons_ucons(prob::AbstractBVProblem, ::Type{T}, M, N) where {T}
-    inequality_length = length(prob.f.lcons)
-    lcons = if isnothing(prob.f.lcons)
+    lcons = if isnothing(prob.lcons)
         zeros(T, N*M)
     else
-        vcat(prob.f.lcons, zeros(T, N*M - inequality_length))
+        lcons_length = length(prob.lcons)
+        vcat(prob.lcons, zeros(T, N*M - lcons_length))
     end
-    ucons = if isnothing(prob.f.ucons)
+    ucons = if isnothing(prob.ucons)
         zeros(T, N*M)
     else
-        vcat(prob.f.ucons, zeros(T, N*M - inequality_length))
+        ucons_length = length(prob.ucons)
+        vcat(prob.ucons, zeros(T, N*M - ucons_length))
     end
     return lcons, ucons
 end

From 532392cf1ec9a6c88100e3c6e95fac1b990deabc Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 16 Aug 2025 23:35:47 +0800
Subject: [PATCH 49/50] Add test case for inequality constraints

---
 .../test/mirk_basic_tests.jl                  | 19 +++++++++++++++++++
 1 file changed, 19 insertions(+)

diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
index 7efa04a3b..a91039765 100644
--- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
@@ -478,3 +478,22 @@ end
         end
     end
 end
+
+@testitem "BVP with inequality constraints" setup=[MIRKConvergenceTests] begin
+    using LinearAlgebra, DiffEqDevTools, OptimizationMOI, Ipopt
+
+    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] = u(pi / 4)[1] + pi / 2
+        residual[2] = u(pi / 2)[1] - pi / 2
+    end
+    prob = BVProblem(simplependulum!, bc!, [pi / 2, pi / 2], tspan,
+        lcons = [0.0, 0.0], ucons = [Inf, Inf])
+    @test_nowarn sol = solve(prob, MIRK4(; optimize = Ipopt.Optimizer()), dt = 0.05)
+end

From ab8ea4afb5e1fba35e80fd3c4899072c12f6a647 Mon Sep 17 00:00:00 2001
From: Qingyu Qu <2283984853@qq.com>
Date: Sat, 16 Aug 2025 23:45:33 +0800
Subject: [PATCH 50/50] Ad edge case handling

---
 lib/BoundaryValueDiffEqCore/src/utils.jl             | 12 ++++++++++++
 lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl |  4 ++--
 2 files changed, 14 insertions(+), 2 deletions(-)

diff --git a/lib/BoundaryValueDiffEqCore/src/utils.jl b/lib/BoundaryValueDiffEqCore/src/utils.jl
index 7969ff83e..b89cfff3f 100644
--- a/lib/BoundaryValueDiffEqCore/src/utils.jl
+++ b/lib/BoundaryValueDiffEqCore/src/utils.jl
@@ -655,6 +655,9 @@ function __construct_internal_problem(prob::AbstractBVProblem, alg, loss, jac,
     # multiple shooting always use iip
     iip = SciMLBase.isinplace(prob)
     if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize))
+        !isnothing(prob.lcons) ||
+            !isnothing(prob.ucons) &&
+            throw(ArgumentError("Inequality constraints can only be handled by optimization solvers, please specify e.g. `MIRK4(; optimize = Ipopt.Optimizer())`."))
         nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype,
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
@@ -673,6 +676,9 @@ function __construct_internal_problem(prob::TwoPointBVProblem, alg, loss, jac,
     T = eltype(y)
     iip = SciMLBase.isinplace(prob)
     if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize))
+        !isnothing(prob.lcons) ||
+            !isnothing(prob.ucons) &&
+            throw(ArgumentError("Inequality constraints can only be handled by optimization solvers, please specify e.g. `MIRK4(; optimize = Ipopt.Optimizer())`."))
         nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype,
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
@@ -691,6 +697,9 @@ function __construct_internal_problem(prob, alg, loss, jac, jac_prototype,
         resid_prototype, y, p, M::Int, N::Int, ::Nothing)
     T = eltype(y)
     if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize))
+        !isnothing(prob.lcons) ||
+            !isnothing(prob.ucons) &&
+            throw(ArgumentError("Inequality constraints can only be handled by optimization solvers, please specify e.g. `MIRK4(; optimize = Ipopt.Optimizer())`."))
         nlf = NonlinearFunction{true}(loss; jac = jac, resid_prototype = resid_prototype,
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
@@ -710,6 +719,9 @@ function __construct_internal_problem(
     T = eltype(y)
     iip = SciMLBase.isinplace(prob)
     if !isnothing(alg.nlsolve) || (isnothing(alg.nlsolve) && isnothing(alg.optimize))
+        !isnothing(prob.lcons) ||
+            !isnothing(prob.ucons) &&
+            throw(ArgumentError("Inequality constraints can only be handled by optimization solvers, please specify e.g. `MIRK4(; optimize = Ipopt.Optimizer())`."))
         nlf = NonlinearFunction{iip}(loss; jac = jac, resid_prototype = resid_prototype,
             jac_prototype = jac_prototype)
         return __internal_nlsolve_problem(prob, resid_prototype, y, nlf, y, p)
diff --git a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
index a91039765..5c0a72474 100644
--- a/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
+++ b/lib/BoundaryValueDiffEqMIRK/test/mirk_basic_tests.jl
@@ -479,8 +479,8 @@ end
     end
 end
 
-@testitem "BVP with inequality constraints" setup=[MIRKConvergenceTests] begin
-    using LinearAlgebra, DiffEqDevTools, OptimizationMOI, Ipopt
+@testitem "BVP with inequality constraints" begin
+    using BoundaryValueDiffEqMIRK, OptimizationMOI, Ipopt
 
     tspan = (0.0, pi / 2)
     function simplependulum!(du, u, p, t)