From fec8df9539e15a1c1725acb885a517f9d5523352 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Thu, 17 Apr 2025 21:35:26 +0200
Subject: [PATCH 1/9] add SciMLSensitivity

---
 docs/Project.toml            | 9 +++++++--
 docs/src/tutorial-goddard.md | 1 +
 test/Project.toml            | 8 +++++++-
 3 files changed, 15 insertions(+), 3 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 022ff0565..d14b28531 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -17,9 +17,11 @@ MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
+SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 
 [compat]
@@ -27,6 +29,8 @@ BenchmarkTools = "1.6"
 CTBase = "0.15"
 CTDirect = "0.14"
 CTFlows = "0.8"
+CTModels = "0.2"
+CTParser = "0.2"
 CommonSolve = "0.2"
 DifferentiationInterface = "0.6"
 Documenter = "1.8"
@@ -34,11 +38,12 @@ DocumenterMermaid = "0.2"
 ForwardDiff = "0.10"
 JLD2 = "0.5"
 JSON3 = "1.14"
+LinearAlgebra = "1"
 MINPACK = "1.3"
 MadNLP = "0.8"
 NLPModelsIpopt = "0.10"
-NonlinearSolve = "4.4"
-OrdinaryDiffEq = "6"
+NonlinearSolve = "4.6"
+OrdinaryDiffEq = "6.93"
 Percival = "0.7"
 Plots = "1.40"
 Suppressor = "0.2"
diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index a301d82dd..1bd977b3b 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -303,6 +303,7 @@ nothing # hide
 Let us do some benchmarking.
 
 ```@example main
+using SciMLSensitivity
 using BenchmarkTools
 @benchmark solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false))
 ```
diff --git a/test/Project.toml b/test/Project.toml
index 199f6cf56..a092c2c2b 100644
--- a/test/Project.toml
+++ b/test/Project.toml
@@ -10,7 +10,13 @@ SplitApplyCombine = "03a91e81-4c3e-53e1-a0a4-9c0c8f19dd66"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [compat]
+CTModels = "0.2"
+DifferentiationInterface = "0.6"
+ForwardDiff = "0.10"
+LinearAlgebra = "1"
+MINPACK = "1.3"
 NLPModelsIpopt = "0.10"
-OrdinaryDiffEq = "6"
+OrdinaryDiffEq = "6.93"
 SplitApplyCombine = "1.2"
+Test = "1"
 julia = "1.10"

From 93d62504cb64c9bc16c47c71caa7d0708ed78d63 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Thu, 17 Apr 2025 22:58:53 +0200
Subject: [PATCH 2/9] fixing goddard tuto

---
 docs/src/tutorial-goddard.md | 53 +++++++++++++++++++-----------------
 1 file changed, 28 insertions(+), 25 deletions(-)

diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index 1bd977b3b..a7f2d6bd1 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -56,12 +56,12 @@ using Plots           # to plot the solution
 We define the problem
 
 ```@example main
-t0 = 0      # initial time
-r0 = 1      # initial altitude
-v0 = 0      # initial speed
-m0 = 1      # initial mass
-vmax = 0.1  # maximal authorized speed
-mf = 0.6    # final mass to target
+const t0 = 0      # initial time
+const r0 = 1      # initial altitude
+const v0 = 0      # initial speed
+const m0 = 1      # initial mass
+const vmax = 0.1  # maximal authorized speed
+const mf = 0.6    # final mass to target
 
 ocp = @def begin # definition of the optimal control problem
 
@@ -70,7 +70,7 @@ ocp = @def begin # definition of the optimal control problem
     x = (r, v, m) ∈ R³, state
     u ∈ R, control
 
-    x(t0) == [ r0, v0, m0 ]
+    x(t0) == [r0, v0, m0]
     m(tf) == mf,         (1)
     0 ≤ u(t) ≤ 1
     r(t) ≥ r0
@@ -80,7 +80,7 @@ ocp = @def begin # definition of the optimal control problem
 
     r(tf) → max
 
-end;
+end
 
 # Dynamics
 const Cd = 310
@@ -91,12 +91,12 @@ const b = 2
 F0(x) = begin
     r, v, m = x
     D = Cd * v^2 * exp(-β*(r - 1)) # Drag force
-    return [ v, -D/m - 1/r^2, 0 ]
+    return [v, -D/m - 1/r^2, 0]
 end
 
 F1(x) = begin
     r, v, m = x
-    return [ 0, Tmax/m, -b*Tmax ]
+    return [0, Tmax/m, -b*Tmax]
 end
 nothing # hide
 ```
@@ -203,13 +203,13 @@ these expressions are straightforwardly translated into Julia code:
 
 ```@example main
 # Controls
-u0 = 0                                  # off control
-u1 = 1                                  # bang control
+const u0 = 0                            # off control
+const u1 = 1                            # bang control
 
 H0 = Lift(F0)                           # H0(x, p) = p' * F0(x)
-H01  = @Lie { H0, H1 }
-H001 = @Lie { H0, H01 }
-H101 = @Lie { H1, H01 }
+H01  = @Lie {H0, H1}
+H001 = @Lie {H0, H01}
+H101 = @Lie {H1, H01}
 us(x, p) = -H001(x, p) / H101(x, p)     # singular control
 
 ub(x) = -(F0⋅g)(x) / (F1⋅g)(x)          # boundary control
@@ -239,7 +239,7 @@ function shoot!(s, p0, t1, t2, t3, tf)
     xf, pf = f0(t3, x3, p3, tf)
 
     s[1] = xf[3] - mf                             # final mass constraint
-    s[2:3] = pf[1:2] - [ 1, 0 ]                   # transversality conditions
+    s[2:3] = pf[1:2] - [1, 0]                     # transversality conditions
     s[4] = H1(x1, p1)                             # H1 = H01 = 0
     s[5] = H01(x1, p1)                            # at the entrance of the singular arc
     s[6] = g(x2)                                  # g = 0 when entering the boundary arc
@@ -283,7 +283,7 @@ println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
 We aggregate the data to define the initial guess vector.
 
 ```@example main
-ξ = [ p0 ; t1 ; t2 ; t3 ; tf ] # initial guess
+ξ = [p0, t1, t2, t3, tf] # initial guess
 ```
 
 ### NonlinearSolve.jl
@@ -343,13 +343,22 @@ println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
 
 ### MINPACK.jl
 
-```@setup main
+Instead of the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package we can use the 
+[MINPACK.jl](https://github.com/sglyon/MINPACK.jl) package to solve 
+the shooting equation. To compute the Jacobian of the shooting function we use the 
+[DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with 
+[ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) backend.
+
+```@example main
 using MINPACK
+```
+
+```@setup main
 function fsolve(f, j, x; kwargs...)
     try
         MINPACK.fsolve(f, j, x; kwargs...)
     catch e
-        println("Erreur using MINPACK")
+        println("Error using MINPACK")
         println(e)
         println("hybrj not supported. Replaced by hybrd even if it is not visible on the doc.")
         MINPACK.fsolve(f, x; kwargs...)
@@ -357,12 +366,6 @@ function fsolve(f, j, x; kwargs...)
 end
 ```
 
-Instead of the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package we can use the 
-[MINPACK.jl](https://github.com/sglyon/MINPACK.jl) package to solve 
-the shooting equation. To compute the Jacobian of the shooting function we use the 
-[DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with 
-[ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) backend.
-
 ```@example main
 using DifferentiationInterface
 import ForwardDiff

From f6ff527466c1bc685fedb2c661c269c1e0021062 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Thu, 17 Apr 2025 23:04:56 +0200
Subject: [PATCH 3/9] add SciMLSensitivity

---
 docs/Project.toml        | 1 +
 docs/src/tutorial-iss.md | 1 +
 2 files changed, 2 insertions(+)

diff --git a/docs/Project.toml b/docs/Project.toml
index d14b28531..2221a7c64 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -46,5 +46,6 @@ NonlinearSolve = "4.6"
 OrdinaryDiffEq = "6.93"
 Percival = "0.7"
 Plots = "1.40"
+SciMLSensitivity = "7.76"
 Suppressor = "0.2"
 julia = "1.10"
diff --git a/docs/src/tutorial-iss.md b/docs/src/tutorial-iss.md
index 251a0c814..99c616b08 100644
--- a/docs/src/tutorial-iss.md
+++ b/docs/src/tutorial-iss.md
@@ -167,6 +167,7 @@ nothing # hide
 Let us do some benchmarking.
 
 ```@example main
+using SciMLSensitivity
 using BenchmarkTools
 @benchmark solve(prob; show_trace=Val(false))
 ```

From c3674fdea6c66cfdd079cd86024bec8d571799dc Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Thu, 17 Apr 2025 23:44:57 +0200
Subject: [PATCH 4/9] remove SciMLSensitivity, use SimpleNewtonRaphson

---
 docs/Project.toml            |  3 ---
 docs/src/tutorial-goddard.md | 16 ++++++----------
 docs/src/tutorial-iss.md     | 31 +++++++++++--------------------
 3 files changed, 17 insertions(+), 33 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index 2221a7c64..fa2e8368c 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -17,11 +17,9 @@ MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
-SciMLSensitivity = "1ed8b502-d754-442c-8d5d-10ac956f44a1"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
 
 [compat]
@@ -46,6 +44,5 @@ NonlinearSolve = "4.6"
 OrdinaryDiffEq = "6.93"
 Percival = "0.7"
 Plots = "1.40"
-SciMLSensitivity = "7.76"
 Suppressor = "0.2"
 julia = "1.10"
diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index a7f2d6bd1..2f5e724c1 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -300,26 +300,22 @@ prob = NonlinearProblem(nle!, ξ)
 nothing # hide
 ```
 
-Let us do some benchmarking.
+Let us do some benchmarking. This will be useful to compare the performance with the MINPACK.jl package below.
 
 ```@example main
-using SciMLSensitivity
 using BenchmarkTools
-@benchmark solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false))
-```
-
-For small nonlinear systems, it could be faster to use the 
-[`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/).
-
-```@example main
 @benchmark solve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(false))
 ```
 
+!!! note
+
+    For small nonlinear systems, the [`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/) may be faster.
+    
 Now, let us solve the problem and retrieve the initial costate solution.
 
 ```@example main
 # resolution of S(ξ) = 0
-indirect_sol = solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(true))
+indirect_sol = solve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(true))
 
 # we retrieve the costate solution together with the times
 p0 = indirect_sol.u[1:3]
diff --git a/docs/src/tutorial-iss.md b/docs/src/tutorial-iss.md
index 99c616b08..7f5b27d2d 100644
--- a/docs/src/tutorial-iss.md
+++ b/docs/src/tutorial-iss.md
@@ -145,8 +145,7 @@ nothing # hide
 
 ## Resolution of the shooting equation
 
-At the end, solving (BVP) is equivalent to solve $S(p_0) = 0$. This is what we call the 
-**indirect simple shooting method**. We define an initial guess.
+At the end, solving (BVP) is equivalent to solve $S(p_0) = 0$. This is what we call the **indirect simple shooting method**. We define an initial guess.
 
 ```@example main
 ξ = [ 0.1 ]    # initial guess
@@ -155,8 +154,7 @@ nothing # hide
 
 ### NonlinearSolve.jl
 
-We first use the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package to solve the shooting
-equation. Let us define the problem.
+We first use the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package to solve the shooting equation. Let us define the problem.
 
 ```@example main
 nle! = (s, ξ, λ) -> s[1] = S(ξ[1])    # auxiliary function
@@ -164,16 +162,17 @@ prob = NonlinearProblem(nle!, ξ)      # NLE problem with initial guess
 nothing # hide
 ```
 
-Let us do some benchmarking.
+Let us do some benchmarking. This will be useful to compare the performance with the MINPACK.jl package below.
 
 ```@example main
 using SciMLSensitivity
 using BenchmarkTools
-@benchmark solve(prob; show_trace=Val(false))
+@benchmark solve(prob, SimpleNewtonRaphson(); show_trace=Val(false))
 ```
 
-For small nonlinear systems, it could be faster to use the 
-[`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/).
+!!! note
+
+    For small nonlinear systems, the [`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/) may be faster.
 
 ```@example main
 @benchmark solve(prob, SimpleNewtonRaphson(); show_trace=Val(false))
@@ -190,6 +189,8 @@ println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
 
 ### MINPACK.jl
 
+Instead of the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package we can use [MINPACK.jl](https://github.com/sglyon/MINPACK.jl) to solve the shooting equation. To compute the Jacobian of the shooting function we use [DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) with [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) backend.
+
 ```@setup main
 using MINPACK
 function fsolve(f, j, x; kwargs...)
@@ -204,12 +205,6 @@ function fsolve(f, j, x; kwargs...)
 end
 ```
 
-Instead of the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package we can use the 
-[MINPACK.jl](https://github.com/sglyon/MINPACK.jl) package to solve 
-the shooting equation. To compute the Jacobian of the shooting function we use the 
-[DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with 
-[ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) backend.
-
 ```@example main
 using DifferentiationInterface
 import ForwardDiff
@@ -225,8 +220,7 @@ jnle! = (js, ξ) -> jacobian!(nle!, similar(ξ), js, backend, ξ)    # Jacobian
 nothing # hide
 ```
 
-We are now in position to solve the problem with the `hybrj` solver from MINPACK.jl through the `fsolve` 
-function, providing the Jacobian. Let us do some benchmarking.
+We are now in position to solve the problem with the `hybrj` solver from MINPACK.jl through the `fsolve`  function, providing the Jacobian. Let us do some benchmarking to compare the performance with the NonlinearSolve.jl package above.
 
 ```@example main
 @benchmark fsolve(nle!, jnle!, ξ; show_trace=false)    # initial guess given to the solver
@@ -258,10 +252,7 @@ sol = φ((t0, tf), x0, p0_sol)
 plot(sol)
 ```
 
-In the indirect shooting method, the research of the optimal control is replaced by the computation
-of its associated extremal. This computation is equivalent to finding the initial covector solution
-to the shooting function. Let us plot the extremal in the phase space and the shooting function with 
-the solution.
+In the indirect shooting method, the research of the optimal control is replaced by the computation of its associated extremal. This computation is equivalent to finding the initial covector solution to the shooting function. Let us plot the extremal in the phase space and the shooting function with  the solution.
 
 ```@raw html
 <article class="docstring">

From f2f476f76c3f19c2ba73ace65275f44fcfa39128 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Thu, 17 Apr 2025 23:45:27 +0200
Subject: [PATCH 5/9] goddard tuto fix

---
 docs/src/tutorial-goddard.md | 62 ++++++++++++------------------------
 1 file changed, 20 insertions(+), 42 deletions(-)

diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index a7f2d6bd1..c94f33a9b 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -229,7 +229,7 @@ Then, we define the shooting function according to the optimal structure we have
 that is a concatenation of four arcs.
 
 ```@example main
-x0 = [ r0, v0, m0 ] # initial state
+x0 = [r0, v0, m0] # initial state
 
 function shoot!(s, p0, t1, t2, t3, tf)
 
@@ -283,7 +283,7 @@ println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
 We aggregate the data to define the initial guess vector.
 
 ```@example main
-ξ = [p0, t1, t2, t3, tf] # initial guess
+ξ = [p0; t1; t2; t3; tf] # initial guess
 ```
 
 ### NonlinearSolve.jl
@@ -305,6 +305,14 @@ Let us do some benchmarking.
 ```@example main
 using SciMLSensitivity
 using BenchmarkTools
+function nlsolve(prob; kwargs...)
+    try
+        NonlinearSolve.solve(prob; kwargs...)
+    catch e
+        println("Error using NonlinearSolve")
+        println(e)
+    end
+end
 @benchmark solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false))
 ```
 
@@ -312,33 +320,15 @@ For small nonlinear systems, it could be faster to use the
 [`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/).
 
 ```@example main
-@benchmark solve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(false))
-```
-
-Now, let us solve the problem and retrieve the initial costate solution.
-
-```@example main
-# resolution of S(ξ) = 0
-indirect_sol = solve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(true))
-
-# we retrieve the costate solution together with the times
-p0 = indirect_sol.u[1:3]
-t1 = indirect_sol.u[4]
-t2 = indirect_sol.u[5]
-t3 = indirect_sol.u[6]
-tf = indirect_sol.u[7]
-
-println("")
-println("p0 = ", p0)
-println("t1 = ", t1)
-println("t2 = ", t2)
-println("t3 = ", t3)
-println("tf = ", tf)
-
-# Norm of the shooting function at solution
-s = similar(p0, 7)
-shoot!(s, p0, t1, t2, t3, tf)
-println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
+function nlsolve(prob, meth; kwargs...)
+    try
+        NonlinearSolve.solve(prob, meth; kwargs...)
+    catch e
+        println("Error using NonlinearSolve")
+        println(e)
+    end
+end
+@benchmark nlsolve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(false))
 ```
 
 ### MINPACK.jl
@@ -349,10 +339,6 @@ the shooting equation. To compute the Jacobian of the shooting function we use t
 [DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with 
 [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) backend.
 
-```@example main
-using MINPACK
-```
-
 ```@setup main
 function fsolve(f, j, x; kwargs...)
     try
@@ -388,15 +374,7 @@ We are now in position to solve the problem with the `hybrj` solver from MINPACK
 function, providing the Jacobian. Let us do some benchmarking.
 
 ```@example main
-@benchmark fsolve(nle!, jnle!, ξ; show_trace=false) # initial guess given to the solver
-```
-
-We can also use the [preparation step](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/stable/tutorial1/#Preparing-for-multiple-gradients) of DifferentiationInterface.jl.
-
-```@example main
-extras = prepare_jacobian(nle!, similar(ξ), backend, ξ)
-jnle_prepared!(js, ξ) = jacobian!(nle!, similar(ξ), js, backend, ξ, extras)
-@benchmark fsolve(nle!, jnle_prepared!, ξ; show_trace=false)
+@benchmark fsolve(nle!, jnle!, ξ; show_trace=true) # initial guess given to the solver
 ```
 
 Now, let us solve the problem and retrieve the initial costate solution.

From a2e28c3de82c877c6371ddf3f3c2c98994d6b201 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Fri, 18 Apr 2025 00:12:51 +0200
Subject: [PATCH 6/9] goddard tuto fix

---
 docs/src/tutorial-goddard.md | 6 +++---
 1 file changed, 3 insertions(+), 3 deletions(-)

diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index c3b7fd742..a00a5b272 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -49,6 +49,7 @@ using OrdinaryDiffEq  # to get the Flow function from OptimalControl
 using NonlinearSolve  # interface to NLE solvers
 using MINPACK         # NLE solver: use to solve the shooting equation
 using Plots           # to plot the solution
+using BenchmarkTools  # do some benchmarking
 ```
 
 ## Optimal control problem
@@ -302,8 +303,7 @@ nothing # hide
 
 Let us do some benchmarking. This will be useful to compare the performance with the MINPACK.jl package below.
 
-```@example main
-using BenchmarkTools
+```@setup main
 function nlsolve(prob; kwargs...)
     try
         NonlinearSolve.solve(prob; kwargs...)
@@ -318,7 +318,7 @@ end
 For small nonlinear systems, it could be faster to use the 
 [`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/).
 
-```@example main
+```@setup main
 function nlsolve(prob, meth; kwargs...)
     try
         NonlinearSolve.solve(prob, meth; kwargs...)

From 5a0213cf8196a0540ff139d805bb07825b800112 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Fri, 18 Apr 2025 00:48:56 +0200
Subject: [PATCH 7/9] fix goddard

---
 docs/Project.toml            | 1 +
 docs/src/tutorial-goddard.md | 2 +-
 docs/src/tutorial-iss.md     | 1 -
 3 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/Project.toml b/docs/Project.toml
index fa2e8368c..4fab10d5f 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -17,6 +17,7 @@ MINPACK = "4854310b-de5a-5eb6-a2a5-c1dee2bd17f9"
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
 NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
 NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OptimalControl = "5f98b655-cc9a-415a-b60e-744165666948"
 OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Percival = "01435c0c-c90d-11e9-3788-63660f8fbccc"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index a00a5b272..7453f5be4 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -284,7 +284,7 @@ println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
 We aggregate the data to define the initial guess vector.
 
 ```@example main
-ξ = [p0; t1; t2; t3; tf] # initial guess
+ξ = [p0..., t1, t2, t3, tf] # initial guess
 ```
 
 ### NonlinearSolve.jl
diff --git a/docs/src/tutorial-iss.md b/docs/src/tutorial-iss.md
index 7f5b27d2d..b63678d52 100644
--- a/docs/src/tutorial-iss.md
+++ b/docs/src/tutorial-iss.md
@@ -165,7 +165,6 @@ nothing # hide
 Let us do some benchmarking. This will be useful to compare the performance with the MINPACK.jl package below.
 
 ```@example main
-using SciMLSensitivity
 using BenchmarkTools
 @benchmark solve(prob, SimpleNewtonRaphson(); show_trace=Val(false))
 ```

From 010ef443ced662bef29951e139a12bf036ff9faf Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Fri, 18 Apr 2025 00:53:12 +0200
Subject: [PATCH 8/9] just do SimpleNewtonRaphson()

---
 docs/src/tutorial-goddard.md | 15 ++-------------
 1 file changed, 2 insertions(+), 13 deletions(-)

diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index 7453f5be4..781e06ee2 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -303,20 +303,9 @@ nothing # hide
 
 Let us do some benchmarking. This will be useful to compare the performance with the MINPACK.jl package below.
 
-```@setup main
-function nlsolve(prob; kwargs...)
-    try
-        NonlinearSolve.solve(prob; kwargs...)
-    catch e
-        println("Error using NonlinearSolve")
-        println(e)
-    end
-end
-@benchmark nlsolve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false))
-```
+!!! note
 
-For small nonlinear systems, it could be faster to use the 
-[`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/).
+    For small nonlinear systems, the [`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/) may be faster.
 
 ```@setup main
 function nlsolve(prob, meth; kwargs...)

From 29f2790b15263775521cb89760ec60e03f2567e0 Mon Sep 17 00:00:00 2001
From: Jean-Baptiste Caillau <jean-baptiste.caillau@univ-cotedazur.fr>
Date: Fri, 18 Apr 2025 01:03:39 +0200
Subject: [PATCH 9/9] goddard + iss tuto fix

---
 docs/src/tutorial-goddard.md | 55 ++---------------------------
 docs/src/tutorial-iss.md     | 68 +++++-------------------------------
 2 files changed, 12 insertions(+), 111 deletions(-)

diff --git a/docs/src/tutorial-goddard.md b/docs/src/tutorial-goddard.md
index a00a5b272..4686e49ab 100644
--- a/docs/src/tutorial-goddard.md
+++ b/docs/src/tutorial-goddard.md
@@ -49,7 +49,6 @@ using OrdinaryDiffEq  # to get the Flow function from OptimalControl
 using NonlinearSolve  # interface to NLE solvers
 using MINPACK         # NLE solver: use to solve the shooting equation
 using Plots           # to plot the solution
-using BenchmarkTools  # do some benchmarking
 ```
 
 ## Optimal control problem
@@ -287,52 +286,9 @@ We aggregate the data to define the initial guess vector.
 ξ = [p0; t1; t2; t3; tf] # initial guess
 ```
 
-### NonlinearSolve.jl
-
-We first use the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package to solve the shooting
-equation. Let us define the problem.
-
-```@example main
-# auxiliary function with aggregated inputs
-nle! = (s, ξ, λ) -> shoot!(s, ξ[1:3], ξ[4], ξ[5], ξ[6], ξ[7])
-
-# NLE problem with initial guess
-prob = NonlinearProblem(nle!, ξ)
-nothing # hide
-```
-
-Let us do some benchmarking. This will be useful to compare the performance with the MINPACK.jl package below.
-
-```@setup main
-function nlsolve(prob; kwargs...)
-    try
-        NonlinearSolve.solve(prob; kwargs...)
-    catch e
-        println("Error using NonlinearSolve")
-        println(e)
-    end
-end
-@benchmark nlsolve(prob; abstol=1e-8, reltol=1e-8, show_trace=Val(false))
-```
-
-For small nonlinear systems, it could be faster to use the 
-[`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/).
-
-```@setup main
-function nlsolve(prob, meth; kwargs...)
-    try
-        NonlinearSolve.solve(prob, meth; kwargs...)
-    catch e
-        println("Error using NonlinearSolve")
-        println(e)
-    end
-end
-@benchmark nlsolve(prob, SimpleNewtonRaphson(); abstol=1e-8, reltol=1e-8, show_trace=Val(false))
-```
-
 ### MINPACK.jl
 
-Instead of the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package we can use the 
+We can use [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package or, instead, the 
 [MINPACK.jl](https://github.com/sglyon/MINPACK.jl) package to solve 
 the shooting equation. To compute the Jacobian of the shooting function we use the 
 [DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) package with 
@@ -370,13 +326,8 @@ nothing # hide
 ```
 
 We are now in position to solve the problem with the `hybrj` solver from MINPACK.jl through the `fsolve` 
-function, providing the Jacobian. Let us do some benchmarking.
-
-```@example main
-@benchmark fsolve(nle!, jnle!, ξ; show_trace=true) # initial guess given to the solver
-```
-
-Now, let us solve the problem and retrieve the initial costate solution.
+function, providing the Jacobian.
+Let us solve the problem and retrieve the initial costate solution.
 
 ```@example main
 # resolution of S(ξ) = 0
diff --git a/docs/src/tutorial-iss.md b/docs/src/tutorial-iss.md
index 7f5b27d2d..d67bcc3f9 100644
--- a/docs/src/tutorial-iss.md
+++ b/docs/src/tutorial-iss.md
@@ -30,11 +30,11 @@ Let us consider the following optimal control problem:
 with $t_0 = 0$, $t_f = 1$, $x_0 = -1$, $x_f = 0$, $\alpha=1.5$ and $\forall\, t \in [t_0, t_f]$, $x(t) \in \R$.
 
 ```@example main
-t0 = 0
-tf = 1
-x0 = -1
-xf = 0
-α  = 1.5
+const t0 = 0
+const tf = 1
+const x0 = -1
+const xf = 0
+const α  = 1.5
 ocp = @def begin
 
     t ∈ [t0, tf], time
@@ -148,48 +148,13 @@ nothing # hide
 At the end, solving (BVP) is equivalent to solve $S(p_0) = 0$. This is what we call the **indirect simple shooting method**. We define an initial guess.
 
 ```@example main
-ξ = [ 0.1 ]    # initial guess
+ξ = [0.1]    # initial guess
 nothing # hide
 ```
 
-### NonlinearSolve.jl
-
-We first use the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package to solve the shooting equation. Let us define the problem.
-
-```@example main
-nle! = (s, ξ, λ) -> s[1] = S(ξ[1])    # auxiliary function
-prob = NonlinearProblem(nle!, ξ)      # NLE problem with initial guess
-nothing # hide
-```
-
-Let us do some benchmarking. This will be useful to compare the performance with the MINPACK.jl package below.
-
-```@example main
-using SciMLSensitivity
-using BenchmarkTools
-@benchmark solve(prob, SimpleNewtonRaphson(); show_trace=Val(false))
-```
-
-!!! note
-
-    For small nonlinear systems, the [`SimpleNewtonRaphson()` descent algorithm](https://docs.sciml.ai/NonlinearSolve/stable/tutorials/code_optimization/) may be faster.
-
-```@example main
-@benchmark solve(prob, SimpleNewtonRaphson(); show_trace=Val(false))
-```
-
-Now, let us solve the problem and retrieve the initial costate solution.
-
-```@example main
-indirect_sol = solve(prob; show_trace=Val(true))      # resolution of S(p0) = 0  
-p0_sol = indirect_sol.u[1]                            # costate solution
-println("\ncostate:    p0 = ", p0_sol)
-println("shoot: |S(p0)| = ", abs(S(p0_sol)), "\n")
-```
-
 ### MINPACK.jl
 
-Instead of the [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package we can use [MINPACK.jl](https://github.com/sglyon/MINPACK.jl) to solve the shooting equation. To compute the Jacobian of the shooting function we use [DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) with [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) backend.
+We can use [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl) package or, instead, [MINPACK.jl](https://github.com/sglyon/MINPACK.jl) to solve the shooting equation. To compute the Jacobian of the shooting function we use [DifferentiationInterface.jl](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface) with [ForwardDiff.jl](https://github.com/JuliaDiff/ForwardDiff.jl) backend.
 
 ```@setup main
 using MINPACK
@@ -220,21 +185,8 @@ jnle! = (js, ξ) -> jacobian!(nle!, similar(ξ), js, backend, ξ)    # Jacobian
 nothing # hide
 ```
 
-We are now in position to solve the problem with the `hybrj` solver from MINPACK.jl through the `fsolve`  function, providing the Jacobian. Let us do some benchmarking to compare the performance with the NonlinearSolve.jl package above.
-
-```@example main
-@benchmark fsolve(nle!, jnle!, ξ; show_trace=false)    # initial guess given to the solver
-```
-
-We can also use the [preparation step](https://gdalle.github.io/DifferentiationInterface.jl/DifferentiationInterface/stable/tutorial1/#Preparing-for-multiple-gradients) of DifferentiationInterface.jl.
-
-```@example main
-extras = prepare_jacobian(nle!, similar(ξ), backend, ξ)
-jnle_prepared!(js, ξ) = jacobian!(nle!, similar(ξ), js, backend, ξ, extras)
-@benchmark fsolve(nle!, jnle_prepared!, ξ; show_trace=false)
-```
-
-Now, let us solve the problem and retrieve the initial costate solution.
+We are now in position to solve the problem with the `hybrj` solver from MINPACK.jl through the `fsolve` function, providing the Jacobian.
+Let us solve the problem and retrieve the initial costate solution.
 
 ```@example main
 indirect_sol = fsolve(nle!, jnle!, ξ; show_trace=true)    # resolution of S(p0) = 0
@@ -400,5 +352,3 @@ nothing # hide
 ```@example main
 pretty_plot(S, p0_sol; size=(800, 450))
 ```
-
-