SciML · ChrisRackauckas · Jun 1, 2025 · Jun 1, 2025
diff --git a/docs/pages.jl b/docs/pages.jl
@@ -2,8 +2,8 @@
 
 pages = ["index.md",
     "Getting Started with BVP solving in Julia" => "tutorials/getting_started.md",
-    "Tutorials" => Any[
-        "tutorials/continuation.md", "tutorials/solve_nlls_bvp.md", "tutorials/extremum.md"],
+    "Tutorials" => Any["tutorials/continuation.md",
+        "tutorials/solve_nlls_bvp.md", "tutorials/extremum.md"],
     "Basics" => Any["basics/bvp_problem.md", "basics/bvp_functions.md",
         "basics/solve.md", "basics/autodiff.md", "basics/error_control.md"],
     "Solver Summaries and Recommendations" => Any[

diff --git a/ext/BoundaryValueDiffEqODEInterfaceExt.jl b/ext/BoundaryValueDiffEqODEInterfaceExt.jl
@@ -61,13 +61,17 @@ function SciMLBase.__solve(prob::BVProblem, alg::BVPM2; dt = 0.0, reltol = 1e-3,
         @closure (t, u, du) -> du .= vec(prob.f(reshape(u, u0_size), prob.p, t))
     end
     bvp2m_bc = if SciMLBase.isinplace(prob)
-        @closure (ya, yb, bca, bcb) -> begin
+        @closure (ya,
+            yb,
+            bca,
+            bcb) -> begin
             prob.f.bc[1](reshape(bca, left_bc_size), reshape(ya, u0_size), prob.p)
             prob.f.bc[2](reshape(bcb, right_bc_size), reshape(yb, u0_size), prob.p)
             return nothing
         end
     else
-        @closure (ya, yb, bca, bcb) -> begin
+        @closure (
+            ya, yb, bca, bcb) -> begin
             bca .= vec(prob.f.bc[1](reshape(ya, u0_size), prob.p))
             bcb .= vec(prob.f.bc[2](reshape(yb, u0_size), prob.p))
             return nothing
@@ -140,15 +144,18 @@ function SciMLBase.__solve(prob::BVProblem, alg::BVPSOL; maxiters = 1000,
     end
 
     bvpsol_bc = if SciMLBase.isinplace(prob)
-        @closure (ya, yb, r) -> begin
+        @closure (ya, yb,
+            r) -> begin
             left_bc = reshape(@view(r[1:no_left_bc]), left_bc_size)
             right_bc = reshape(@view(r[(no_left_bc + 1):end]), right_bc_size)
             prob.f.bc[1](left_bc, reshape(ya, u0_size), prob.p)
             prob.f.bc[2](right_bc, reshape(yb, u0_size), prob.p)
             return nothing
         end
     else
-        @closure (ya, yb, r) -> begin
+        @closure (ya,
+            yb,
+            r) -> begin
             r[1:no_left_bc] .= vec(prob.f.bc[1](reshape(ya, u0_size), prob.p))
             r[(no_left_bc + 1):end] .= vec(prob.f.bc[2](reshape(yb, u0_size), prob.p))
             return nothing

diff --git a/lib/BoundaryValueDiffEqAscher/src/adaptivity.jl b/lib/BoundaryValueDiffEqAscher/src/adaptivity.jl
@@ -188,15 +188,13 @@ function error_estimate!(cache::AscherCache)
         # in valstr in case they prove to be needed later for an error estimate.
         x = mesh[i] + (mesh_dt[i]) * 2.0 / 3.0
         @views approx(cache, x, valstr[i][3])
-        error[i] .= wgterr .*
-                    abs.(valstr[i][3] .-
+        error[i] .= wgterr .* abs.(valstr[i][3] .-
                          (isodd(i) ? valstr[Int((i + 1) / 2)][2] : valstr[Int(i / 2)][4]))
 
         x = mesh[i] + (mesh_dt[i]) / 3.0
         @views approx(cache, x, valstr[i][2])
         error[i] .= error[i] .+
-                    wgterr .*
-                    abs.(valstr[i][2] .-
+                    wgterr .* abs.(valstr[i][2] .-
                          (isodd(i) ? valstr[Int((i + 1) / 2)][1] : valstr[Int(i / 2)][3]))
     end
     return maximum(reduce(hcat, error), dims = 2)

diff --git a/lib/BoundaryValueDiffEqAscher/src/ascher.jl b/lib/BoundaryValueDiffEqAscher/src/ascher.jl
@@ -102,7 +102,8 @@ function SciMLBase.__init(
 
     iip = isinplace(prob)
 
-    f, bc = if prob.u0 isa AbstractVector
+    f,
+    bc = if prob.u0 isa AbstractVector
         prob.f, prob.f.bc
     elseif iip
         vecf! = @closure (du, u, p, t) -> __vec_f!(du, u, p, t, prob.f, size(u0))
@@ -337,10 +338,11 @@ function __construct_nlproblem(cache::AscherCache{iip, T}) where {iip, T}
     end
 
     jac = if iip
-        @closure (J, u, p) -> __ascher_mpoint_jacobian!(
-            J, u, diffmode, jac_cache, loss, lz, cache.p)
+        @closure (J, u,
+            p) -> __ascher_mpoint_jacobian!(J, u, diffmode, jac_cache, loss, lz, cache.p)
     else
-        @closure (u, p) -> __ascher_mpoint_jacobian(
+        @closure (u,
+            p) -> __ascher_mpoint_jacobian(
             jac_prototype, u, diffmode, jac_cache, loss, cache.p)
     end
 

diff --git a/lib/BoundaryValueDiffEqAscher/src/collocation.jl b/lib/BoundaryValueDiffEqAscher/src/collocation.jl
@@ -1,5 +1,6 @@
 function Φ!(cache::AscherCache{iip, T}, z, res, pt::StandardBVProblem) where {iip, T}
-    (; f, mesh, mesh_dt, ncomp, ny, bc, k, p, zeta, residual, zval, yval, gval, delz, dmz, deldmz, g, w, v, ipvtg, ipvtw, TU) = cache
+    (; f, mesh, mesh_dt, ncomp, ny, bc, k, p, zeta, residual, zval,
+        yval, gval, delz, dmz, deldmz, g, w, v, ipvtg, ipvtw, TU) = cache
     (; acol, rho) = TU
     ncy = ncomp + ny
     n = length(mesh) - 1
@@ -156,7 +157,8 @@ function Φ!(cache::AscherCache{iip, T}, z, res, pt::StandardBVProblem) where {i
 end
 
 function Φ!(cache::AscherCache{iip, T}, z, res, pt::TwoPointBVProblem) where {iip, T}
-    (; f, mesh, mesh_dt, ncomp, ny, bc, k, p, zeta, bcresid_prototype, residual, zval, yval, gval, delz, dmz, deldmz, g, w, v, dmzo, ipvtg, ipvtw, TU) = cache
+    (; f, mesh, mesh_dt, ncomp, ny, bc, k, p, zeta, bcresid_prototype, residual,
+        zval, yval, gval, delz, dmz, deldmz, g, w, v, dmzo, ipvtg, ipvtw, TU) = cache
     (; acol, rho) = TU
     ncy = ncomp + ny
     n = length(mesh) - 1
@@ -319,7 +321,8 @@ end
 @inline __get_value(z) = isa(z, ForwardDiff.Dual) ? z.value : z
 
 function Φ(cache::AscherCache{iip, T}, z, pt::StandardBVProblem) where {iip, T}
-    (; f, mesh, mesh_dt, ncomp, ny, bc, k, p, zeta, residual, zval, yval, gval, delz, dmz, deldmz, g, w, v, dmzo, ipvtg, ipvtw, TU) = cache
+    (; f, mesh, mesh_dt, ncomp, ny, bc, k, p, zeta, residual, zval, yval,
+        gval, delz, dmz, deldmz, g, w, v, dmzo, ipvtg, ipvtw, TU) = cache
     (; acol, rho) = TU
     ncy = ncomp + ny
     n = length(mesh) - 1
@@ -476,7 +479,8 @@ function Φ(cache::AscherCache{iip, T}, z, pt::StandardBVProblem) where {iip, T}
 end
 
 function Φ(cache::AscherCache{iip, T}, z, pt::TwoPointBVProblem) where {iip, T}
-    (; f, mesh, mesh_dt, ncomp, ny, bc, k, p, zeta, residual, zval, yval, gval, delz, dmz, deldmz, g, w, v, dmzo, ipvtg, ipvtw, TU) = cache
+    (; f, mesh, mesh_dt, ncomp, ny, bc, k, p, zeta, residual, zval, yval,
+        gval, delz, dmz, deldmz, g, w, v, dmzo, ipvtg, ipvtw, TU) = cache
     (; acol, rho) = TU
     ncy = ncomp + ny
     n = length(mesh) - 1

diff --git a/lib/BoundaryValueDiffEqCore/src/types.jl b/lib/BoundaryValueDiffEqCore/src/types.jl
@@ -6,7 +6,8 @@
 end
 
 @inline __materialize_jacobian_algorithm(_, alg::BVPJacobianAlgorithm) = alg
-@inline __materialize_jacobian_algorithm(_, alg::ADTypes.AbstractADType) = BVPJacobianAlgorithm(alg)
+@inline __materialize_jacobian_algorithm(
+    _, alg::ADTypes.AbstractADType) = BVPJacobianAlgorithm(alg)
 @inline __materialize_jacobian_algorithm(::Nothing, ::Nothing) = BVPJacobianAlgorithm()
 @inline function __materialize_jacobian_algorithm(nlsolve::N, ::Nothing) where {N}
     ad = hasfield(N, :jacobian_ad) ? nlsolve.jacobian_ad : missing

diff --git a/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl b/lib/BoundaryValueDiffEqFIRK/src/BoundaryValueDiffEqFIRK.jl
@@ -100,6 +100,7 @@ include("sparse_jacobians.jl")
 
     @compile_workload begin
         @sync for prob in probs, alg in algs
+
             Threads.@spawn solve(prob, alg; dt = 0.2)
         end
     end
@@ -112,6 +113,7 @@ include("sparse_jacobians.jl")
 
     @compile_workload begin
         @sync for prob in probs, alg in algs
+
             Threads.@spawn solve(prob, alg; dt = 0.2)
         end
     end
@@ -124,6 +126,7 @@ include("sparse_jacobians.jl")
 
     @compile_workload begin
         @sync for prob in probs, alg in algs
+
             Threads.@spawn solve(prob, alg; dt = 0.2)
         end
     end
@@ -136,6 +139,7 @@ include("sparse_jacobians.jl")
 
     @compile_workload begin
         @sync for prob in probs, alg in algs
+
             Threads.@spawn solve(prob, alg; dt = 0.2)
         end
     end
@@ -194,6 +198,7 @@ include("sparse_jacobians.jl")
 
     @compile_workload begin
         @sync for prob in probs, alg in algs
+
             Threads.@spawn solve(prob, alg; dt = 0.2, abstol = 1e-2)
         end
     end
@@ -210,6 +215,7 @@ include("sparse_jacobians.jl")
 
     @compile_workload begin
         @sync for prob in probs, alg in algs
+
             Threads.@spawn solve(prob, alg; dt = 0.2, abstol = 1e-2)
         end
     end
@@ -226,6 +232,7 @@ include("sparse_jacobians.jl")
 
     @compile_workload begin
         @sync for prob in probs, alg in algs
+
             Threads.@spawn solve(prob, alg; dt = 0.2, abstol = 1e-2)
         end
     end
@@ -242,6 +249,7 @@ include("sparse_jacobians.jl")
 
     @compile_workload begin
         @sync for prob in probs, alg in algs
+
             Threads.@spawn solve(prob, alg; dt = 0.2, abstol = 1e-2)
         end
     end
@@ -258,6 +266,7 @@ include("sparse_jacobians.jl")
 
     @compile_workload begin
         @sync for prob in probs, alg in algs
+
             Threads.@spawn solve(prob, alg; dt = 0.2, abstol = 1e-2)
         end
     end

diff --git a/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl b/lib/BoundaryValueDiffEqFIRK/src/adaptivity.jl
@@ -226,8 +226,8 @@ function s_constraints(M, h)
         row_start = (i - 1) * M + 1
         for k in 0:(M - 1)
             for j in 1:6
-                A[row_start + k, j + k * 6] = j == 1.0 ? 0.0 :
-                                              (j - 1) * t[i + k * 6]^(j - 2)
+                A[row_start + k,
+                    j + k * 6] = j == 1.0 ? 0.0 : (j - 1) * t[i + k * 6]^(j - 2)
             end
         end
     end

diff --git a/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl b/lib/BoundaryValueDiffEqFIRK/src/algorithms.jl
@@ -95,8 +95,14 @@ for stage in (1, 2, 3, 5, 7)
             defect_threshold::T = 0.1
             max_num_subintervals::Int = 3000
         end
-        $(alg)(nlsolve::N, jac_alg::J; nested = false, nested_nlsolve_kwargs::NamedTuple = (;), defect_threshold::T = 0.1, max_num_subintervals::Int = 3000) where {N, J, T} = $(alg){
-            N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
+        $(alg)(nlsolve::N,
+            jac_alg::J;
+            nested = false,
+            nested_nlsolve_kwargs::NamedTuple = (;),
+            defect_threshold::T = 0.1,
+            max_num_subintervals::Int = 3000) where {N,
+            J,
+            T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
             defect_threshold, max_num_subintervals)
     end
 end
@@ -194,8 +200,14 @@ for stage in (2, 3, 4, 5)
             defect_threshold::T = 0.1
             max_num_subintervals::Int = 3000
         end
-        $(alg)(nlsolve::N, jac_alg::J; nested = false, nested_nlsolve_kwargs::NamedTuple = (;), defect_threshold::T = 0.1, max_num_subintervals::Int = 3000) where {N, J, T} = $(alg){
-            N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
+        $(alg)(nlsolve::N,
+            jac_alg::J;
+            nested = false,
+            nested_nlsolve_kwargs::NamedTuple = (;),
+            defect_threshold::T = 0.1,
+            max_num_subintervals::Int = 3000) where {N,
+            J,
+            T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
             defect_threshold, max_num_subintervals)
     end
 end
@@ -295,8 +307,14 @@ for stage in (2, 3, 4, 5)
             defect_threshold::T = 0.1
             max_num_subintervals::Int = 3000
         end
-        $(alg)(nlsolve::N, jac_alg::J; nested = false, nested_nlsolve_kwargs::NamedTuple = (;), defect_threshold::T = 0.1, max_num_subintervals::Int = 3000) where {N, J, T} = $(alg){
-            N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
+        $(alg)(nlsolve::N,
+            jac_alg::J;
+            nested = false,
+            nested_nlsolve_kwargs::NamedTuple = (;),
+            defect_threshold::T = 0.1,
+            max_num_subintervals::Int = 3000) where {N,
+            J,
+            T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
             defect_threshold, max_num_subintervals)
     end
 end
@@ -396,8 +414,14 @@ for stage in (2, 3, 4, 5)
             defect_threshold::T = 0.1
             max_num_subintervals::Int = 3000
         end
-        $(alg)(nlsolve::N, jac_alg::J; nested = false, nested_nlsolve_kwargs::NamedTuple = (;), defect_threshold::T = 0.1, max_num_subintervals::Int = 3000) where {N, J, T} = $(alg){
-            N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
+        $(alg)(nlsolve::N,
+            jac_alg::J;
+            nested = false,
+            nested_nlsolve_kwargs::NamedTuple = (;),
+            defect_threshold::T = 0.1,
+            max_num_subintervals::Int = 3000) where {N,
+            J,
+            T} = $(alg){N, J, T}(nlsolve, jac_alg, nested, nested_nlsolve_kwargs,
             defect_threshold, max_num_subintervals)
     end
 end