Update README

Author: ErikQQY

README.md

@@ -37,21 +37,6 @@ 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/).

docs/src/tutorials/optimal_control.md

@@ -0,0 +1 @@
+

lib/BoundaryValueDiffEqCore/src/utils.jl

@@ -652,7 +652,6 @@ algorithm used. It returns either a `NonlinearProblem` or an `OptimizationProble
 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,

test/misc/dynamic_optimization.jl

@@ -0,0 +1,71 @@
+@testitem "Linear systems" begin
+    using BoundaryValueDiffEq
+    using OptimizationMOI, Ipopt
+end
+
+@testitem "Rocket launch" begin
+    using BoundaryValueDiffEq
+    using 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] + big(pi / 2)
+        residual[2] = u(pi / 2)[1] - big(pi / 2)
+    end
+    u0 = BigFloat.([pi / 2, pi / 2])
+    multi_point_bvp = BVProblem(simplependulum!, bc!, u0, tspan)
+
+    @testset "BigFloat compatibility with Multi-point BVP" begin
+        for solver in [MIRK4(), RadauIIa5(), LobattoIIIa4(nested_nlsolve = true)]
+            sol = solve(multi_point_bvp, solver, dt = 0.05)
+            @test SciMLBase.successful_retcode(sol.retcode)
+        end
+    end
+
+    function f!(du, u, p, t)
+        du[1] = u[2]
+        du[2] = u[1]
+    end
+    function bca!(resid_a, u_a, p)
+        resid_a[1] = u_a[1] - 1
+    end
+    function bcb!(resid_b, u_b, p)
+        resid_b[1] = u_b[1]
+    end
+    bvp_function = BVPFunction(f!, (bca!, bcb!), bcresid_prototype = (zeros(1), zeros(1)), twopoint = Val(true))
+    tspan = (0.0, 1.0)
+    two_point_bvp = BVProblem(bvp_function, BigFloat.([1.0, 0.0]), tspan)
+
+    @testset "BigFloat compatibility with Two-point BVP" begin
+        for solver in [MIRK4(), RadauIIa5(), LobattoIIIa4(nested_nlsolve = true)]
+            sol = solve(two_point_bvp, solver, dt = 0.05)
+            @test SciMLBase.successful_retcode(sol.retcode)
+        end
+    end
+
+    function second_f!(ddu, du, u, p, t)
+        ϵ = 0.1
+        ddu[1] = u[2]
+        ddu[2] = (-u[1] * du[2] - u[3] * du[3]) / ϵ
+        ddu[3] = (du[1] * u[3] - u[1] * du[3]) / ϵ
+    end
+    function second_bc!(res, du, u, p, t)
+        res[1] = u(0.0)[1]
+        res[2] = u(1.0)[1]
+        res[3] = u(0.0)[3] + 1
+        res[4] = u(1.0)[3] - 1
+        res[5] = du(0.0)[1]
+        res[6] = du(1.0)[1]
+    end
+    u0 = BigFloat.([1.0, 1.0, 1.0])
+    tspan = (0.0, 1.0)
+    prob = SecondOrderBVProblem(second_f!, second_bc!, u0, tspan)
+    @test_broken sol4 = solve(prob, MIRKN4(), dt = 0.01)
+    @test_broken SciMLBase.successful_retcode(sol4.retcode)
+end