diff --git a/docs/make.jl b/docs/make.jl
index 3154259..37e6bf6 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -17,7 +17,7 @@ cp(
repo_url = "github.com/control-toolbox/Tutorials.jl"
makedocs(;
- draft=false, # if draft is true, then the julia code from .md is not executed
+ draft=true, # if draft is true, then the julia code from .md is not executed
# to disable the draft mode in a specific markdown file, use the following:
# ```@meta
# Draft = false
@@ -38,6 +38,12 @@ makedocs(;
"Tutorials and Advanced Features" => [
"Discrete continuation" => "tutorial-continuation.md",
"Discretisation methods" => "tutorial-discretisation.md",
+ "Free times" => [
+ "Final time" => "tutorial-free-times-final.md",
+ "Initial time" => "tutorial-free-times-initial.md",
+ "Final and initial times" => "tutorial-free-times-final-initial.md",
+ "Orbital transfer min time" => "tutorial-free-times-orbital.md",
+ ],
"NLP manipulations" => "tutorial-nlp.md",
"Indirect simple shooting" => "tutorial-iss.md",
"Goddard: direct, indirect" => "tutorial-goddard.md",
diff --git a/docs/src/tutorial-continuation.md b/docs/src/tutorial-continuation.md
index eefd682..5286ba6 100644
--- a/docs/src/tutorial-continuation.md
+++ b/docs/src/tutorial-continuation.md
@@ -1,3 +1,4 @@
+
# Discrete continuation
By using the warm start option, it is easy to implement a basic discrete continuation method, in which a sequence of problems is solved by using each solution as the initial guess for the next problem.
@@ -53,7 +54,7 @@ Then we perform the continuation with a simple *for* loop, using each solution t
init = ()
data = DataFrame(T=Float64[], Objective=Float64[], Iterations=Int[])
for T ∈ range(1, 2, length=5)
- ocp = problem(T)
+ ocp = problem(T)
sol = solve(ocp; init=init, display=false)
global init = sol
push!(data, (T=T, Objective=objective(sol), Iterations=iterations(sol)))
diff --git a/docs/src/tutorial-free-times-final-initial.md b/docs/src/tutorial-free-times-final-initial.md
new file mode 100644
index 0000000..e326e92
--- /dev/null
+++ b/docs/src/tutorial-free-times-final-initial.md
@@ -0,0 +1,54 @@
+```@meta
+Draft = false
+```
+
+# [Optimal control problem with free initial and free final times](@id tutorial-free-times-final-initial)
+
+In this tutorial, we explore optimal control problems with free initial time `t₀` and final time `t_f`.
+
+
+## Here are the required packages for the tutorial:
+
+```@example both_time
+using OptimalControl
+using NLPModelsIpopt
+using BenchmarkTools
+using DataFrames
+using Plots
+using Printf
+```
+## Definition of the problem
+
+```@example both_time
+
+@def ocp begin
+ v=(t0, tf) ∈ R², variable
+ t ∈ [t0, tf], time
+ x ∈ R², state
+ u ∈ R, control
+ -1 ≤ u(t) ≤ 1
+ x(t0) == [0, 0]
+ x(tf) == [1, 0]
+ 0.05 ≤ t0 ≤ 10
+ 0.05 ≤ tf ≤ 10
+ 0.01 ≤ tf - t0 ≤ Inf
+ ẋ(t) == [x₂(t), u(t)]
+ t0 → max
+end
+
+nothing # hide
+```
+
+# Direct resolution with both final and inital times being free:
+
+
+We now solve the problem using a direct method, with automatic treatment of the free initial time.
+
+```@example both_time
+sol = solve(ocp; grid_size=100)
+```
+And plot the solution.
+
+```@example both_time
+plot(sol; label="direct", size=(800, 800))
+```
diff --git a/docs/src/tutorial-free-times-final.md b/docs/src/tutorial-free-times-final.md
new file mode 100644
index 0000000..f047e6a
--- /dev/null
+++ b/docs/src/tutorial-free-times-final.md
@@ -0,0 +1,351 @@
+```@meta
+Draft = false
+```
+
+# [Optimal control problem with free final time](@id tutorial-free-times-final)
+
+In this tutorial, we explore an optimal control problem with free final time `tf`.
+
+Here is the required packages for the tutorial:
+
+```@example main-free-final
+using LinearAlgebra: norm
+using NLPModelsIpopt
+using NonlinearSolve
+using OptimalControl
+using OrdinaryDiffEq # to get the Flow function from OptimalControl
+using Plots
+using Printf
+```
+
+## Definition of the problem
+
+```@example main-free-final
+t0 = 0 # initial time
+x0 = [0, 0] # initial point
+xf_target = [1, 0] # final point
+
+@def ocp begin
+ tf ∈ R, variable
+ t ∈ [t0, tf], time
+ x ∈ R², state
+ u ∈ R, control
+ -1 ≤ u(t) ≤ 1
+ x(t0) == x0
+ x(tf) == xf_target
+ 0.05 ≤ tf ≤ Inf
+ ẋ(t) == [x₂(t), u(t)]
+ tf → min
+end
+nothing # hide
+```
+
+In the definition, the line that shows that the problem has a **free** final time is the following:
+
+```julia
+tf ∈ R, variable
+```
+
+The variable `tf` is indeed the final time according to:
+
+```julia
+t ∈ [t0, tf], time
+```
+
+## Direct resolution:
+
+You may need to add this type of constrain in your problem definition :
+
+To help the direct solver to converge, we have added a constraint on the variable `tf`:
+
+```julia
+0.05 ≤ tf ≤ Inf
+```
+
+We now solve the problem using a direct method, with automatic treatment of the free final time.
+
+```@example main-free-final
+sol = solve(ocp; grid_size=100)
+```
+
+And plot the solution.
+
+```@example main-free-final
+plt = plot(sol; label="direct", size=(800, 800))
+```
+
+## Verification of results
+
+!!! note "Mathematical computations"
+
+ ```@raw html
+
+ Click to show/hide mathematical computations.
+ ```
+
+ Here is the theorical part. The pseudo-Hamiltonian is:
+
+ ```math
+ H(x, p, u) = p_1 x_2 + p_2 u.
+ ```
+
+ Pontryagin's theorem gives:
+
+ ```math
+ \begin{aligned}
+ \dot{p}_1(t) &= 0, \\
+ \dot{p}_2(t) &= -p_1(t). \\
+ \end{aligned}
+ ```
+
+ Hence, $p_1(t) = \mathrm{cst} = \alpha$ and $p_2(t) = -\alpha t + \beta$, $\beta = p_2(0)$. According to Pontryagin's theorem, we have:
+
+ ```math
+ \begin{aligned}
+ p_1(t) &= 1, \\
+ p_2(t) &= 1-t. \\
+ \end{aligned}
+ ```
+
+ Besides, the pseudo-Hamiltonian satisfies:
+ ```math
+ H(x(t_f), p(t_f), u(t_f)) = -p° = 1.
+ ```
+
+ For this problem, we can prove that the optimal control $u$ satisfies:
+
+ ```math
+ u(t) = \left\{
+ \begin{aligned}
+ 1 & \quad\text{if}\quad t \in [0, t_1], \\
+ -1 & \quad\text{if}\quad t \in (t_1, t_f], \\
+ \end{aligned}
+ \right.
+ ```
+
+ where $t_1$ is a constant to determined, as $t_f$. To do so, we have to compute the state $x$:
+
+ On t ∈ [0,t1] :
+ ```math
+ \begin{aligned}
+ \dot{x}_1(t) &= x_2(t) \\
+ \dot{x}_2(t) &= 1 \\
+ \end{aligned}
+ ```
+
+ ```math
+ \begin{aligned}
+ x_2(t) &= t \\
+ x_1(t) &= (1/2)*t^2
+ \end{aligned}
+
+ ```
+
+ When t = t1 :
+ ```math
+ \begin{aligned}
+ x_2(t1) &= t_1 \\
+ x_1(t1) &= (1/2)*t_1^2
+ \end{aligned}
+ ```
+
+ On $t \in [t_1,t_f]$
+ ```math
+ \begin{aligned}
+ x_2(t) &= -t + 2t_1 \\
+ x_1(t) &= -\frac{1}{2}t^2 + 2t_1 t + C_2 \\
+ \end{aligned}
+ ```
+ ```math
+ \begin{aligned}
+ x_1(t_1) &= -\frac{1}{2}t_1^2 + 2t_1^2 + C_2 &= \frac{1}{2}t_1^2 \\
+ C_2 &= -t_1^2 \\
+ \end{aligned}
+ ```
+
+ ```math
+ \begin{aligned}
+ x_1(t) &= -\frac{1}{2}t^2 + 2t_1 t - t_1^2
+ \end{aligned}
+ ```
+ Finally you can solve the terminal conditions :
+ ```math
+ \begin{aligned}
+ x_1(t_f) &= 1, \quad x_2(t_f) &= 0 \\
+ \end{aligned}
+ ```
+ ```math
+ \begin{aligned}
+ x_2(t_f) &= -t_f + 2t_1 &= 0 \\
+ t_f &= 2t_1 \\
+ \end{aligned}
+ ```
+ ```math
+ \begin{aligned}
+ x_1(t_f)& = -\frac{1}{2}t_f^2 + 2t_1 t_f - t_1^2 &= 1 \\
+ \end{aligned}
+ ```
+ and
+ ```math
+ \begin{aligned}
+ -2t_1^2 + 4t_1^2 - t_1^2 &= t_1^2 &= 1 \\
+ \end{aligned}
+ ```
+ gets us
+ ```math
+ \begin{aligned}
+ t_1 &= 1, \quad t_f &= 2
+ \end{aligned}
+ ```
+ To sum up we find the following solutions :
+ ```math
+
+ x_1(t) = \begin{cases}
+ \frac{1}{2}t^2 & \text{si } t \in [0, 1) \\
+ -\frac{1}{2}t^2 + 2t - 1 & \text{si } t \in [1, 2]
+ \end{cases}
+ \qquad
+ x_2(t) = \begin{cases}
+ t & \text{si } t \in [0, 1) \\
+ 2 - t & \text{si } t \in [1, 2]
+ \end{cases}
+ ```
+ ```math
+ p_1(t) = 1, \quad p_2(t) = 1-t, \quad p^0=-1
+ ```
+ ```math
+ u(t) = \begin{cases}
+ 1 & \text{si } t \in [0, 1) \\
+ -1 & \text{si } t \in [1, 2]
+ \end{cases}
+ ```
+ ```@raw html
+
+ ```
+
+Now we can compare the results found with the direct method whith the theoritical analysis:
+```@example main-free-final
+tf = variable(sol)
+u = control(sol)
+p = costate(sol)
+x = state(sol)
+p° = -1
+H(t) = p(t)[1]*x(t)[2] + p(t)[2]*u(t)
+
+@printf("H(tf) = %.3f\n", H(tf))
+@printf("x(tf) = [%.3f, %.3f]\n", x(tf)[1], x(tf)[2])
+@printf("p(tf) = [%.3f, %.3f]\n", p(tf)[1], p(tf)[2])
+```
+
+The numerical results closely match the theoretical predictions: the final state $x(t_f)=[1,0]$ is exactly satisfied.
+The costate and Hamiltonian values at final time show a small deviation (≈ 0.01), likely due to numerical precision.
+Overall, the direct method confirms the theoretical analysis with excellent accuracy.
+
+We can analyse the influence of using different discretization sizes (`grid_size`), and observed the following results for the optimal $t_f$:
+
+```@example main-free-final
+for N in [20, 50, 100, 200]
+ solN = solve(ocp; grid_size=N, display=false)
+ @printf("grid_size = %3d → tf = %.5f\n", N, objective(solN))
+end
+```
+
+This example shows that problems with a free final time can be sensitive to discretization. A small grid may lead to suboptimal or slightly inaccurate results.
+
+## Indirect method
+
+We first define the pseudo-Hamiltonian and the switching function needed to define the shooting function.
+
+```@example main-free-final
+# Pseudo-Hamitlonian
+H(x, p, u) = p[1]*x[2] + p[2]*u
+
+# Hamiltonian lift
+F1(x) = [0, 1]
+H1 = Lift(F1)
+nothing # hide
+```
+
+Then, we define the different flows, associated to control laws $u(t)=+1$ and $u(t)=-1$.
+
+```@example main-free-final
+const u_pos = 1
+const u_neg = -1
+
+f_pos = Flow(ocp, (x, p, tf) -> u_pos) # it depends on tf since tf is the variable
+f_neg = Flow(ocp, (x, p, tf) -> u_neg)
+nothing # hide
+```
+
+We can now define the shooting function following the structure given by the direct method.
+
+```@example main-free-final
+function shoot!(s, p0, t1, tf)
+
+ # the flows
+ x1, p1 = f_pos(t0, x0, p0, t1)
+ xf, pf = f_neg(t1, x1, p1, tf)
+
+ # final control
+ uf = -1
+
+ # shooting conditions
+ s[1:2] = xf .- xf_target # reach the target
+ s[3] = H(xf, pf, uf) - 1 # final condition on the pseudo-Hamiltonian
+ s[4] = H1(x1, p1) # switching condition
+
+end
+nothing # hide
+```
+
+Before solving our problem we must find a good initial guess to help the convergence of the algorithm.
+
+```@example main-free-final
+t = time_grid(sol)
+x = state(sol)
+p = costate(sol)
+φ(t) = H1(x(t), p(t)) # switching function
+
+p0 = p(t0)
+t1 = t[argmin(abs.(φ.(t)))]
+tf = t[end]
+
+println("p0 = ", p0)
+println("t1 = ", t1)
+println("tf = ", tf)
+
+# Norm of the shooting function at initial guess
+s = zeros(4)
+shoot!(s, p0, t1, tf)
+println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
+```
+
+And we finally can solve the problem using an indirect method.
+
+```@example main-free-final
+# Aggregated function
+nle! = (s, ξ, λ) -> shoot!(s, ξ[1:2], ξ[3], ξ[4])
+
+# Nonlinear problem
+ξ_guess = [p0..., t1, tf] # Initial guess
+prob = NonlinearProblem(nle!, ξ_guess)
+
+# Resolution
+indirect_sol = solve(prob; show_trace=Val(true), abstol=1e-8, reltol=1e-8)
+```
+
+We can now plot and compare with the direct method.
+
+```@example main-free-final
+# Data
+p0 = indirect_sol.u[1:2]
+t1 = indirect_sol.u[3]
+tf = indirect_sol.u[4]
+
+# Compute the optimal solution
+f = f_pos * (t1, f_neg) # concatenate the flows
+flow_sol = f((t0, tf), x0, p0; saveat=range(t0, tf, 100))
+
+# Plot the solution
+plot!(plt, flow_sol; label="indirect", color=:red)
+```
\ No newline at end of file
diff --git a/docs/src/tutorial-free-times-initial.md b/docs/src/tutorial-free-times-initial.md
new file mode 100644
index 0000000..9b9a70d
--- /dev/null
+++ b/docs/src/tutorial-free-times-initial.md
@@ -0,0 +1,239 @@
+```@meta
+Draft = false
+```
+
+# [Optimal control problem with free initial time](@id tutorial-free-times-initial)
+
+
+
+In this tutorial, we explore an optimal control problem with free initial time `t0`.
+
+Here are the required packages for the tutorial:
+
+```@example initial_time
+using LinearAlgebra: norm
+using OptimalControl
+using NLPModelsIpopt
+using BenchmarkTools
+using DataFrames
+using Plots
+using Printf
+```
+## Definition of the problem
+```@example initial_time
+
+@def ocp begin
+ t0 ∈ R, variable
+ tf=0
+ t ∈ [t0, tf], time
+ x ∈ R², state
+ u ∈ R, control
+ -1 ≤ u(t) ≤ 1
+ x(t0) == [0, 0]
+ x(tf) == [1, 0]
+ 0.05 ≤ -t0 ≤ Inf
+ ẋ(t) == [x₂(t), u(t)]
+ -t0 → min
+end
+nothing # hide
+```
+
+# Direct resolution with free initial time :
+
+We now solve the problem using a direct method, with automatic treatment of the free initial time.
+
+```@example initial_time
+sol = solve(ocp; grid_size=100)
+```
+And plot the solution.
+
+```@example initial_time
+plot(sol; label="direct", size=(800, 800))
+```
+
+
+
+## Verification of results
+!!! note "Mathematical computations"
+ ```@raw html
+
+ Click to show/hide mathematical verification.
+ ```
+ Here is the theoretical part using Pontryagin's Maximum Principle:
+ ```math
+ H = p_1 x_2 + p_2 u + 1
+ ```
+
+ Conditions from Pontryagin’s theorem:
+ ```math
+ \begin{aligned}
+
+ \dot{p}_1(t) &= 0, \quad \Rightarrow \quad p_1 = c_1 \quad (\text{constant}) \\
+
+ \dot{p}_2(t) &= -p_1 \quad \Rightarrow \quad p_2 = -c_1 t + c_2
+ \end{aligned}
+ ```
+
+ Switching condition:
+ ```math
+ p_2(t_s) = 0 \quad \Rightarrow \quad c_2 = c_1 t_s
+ ```
+
+ Optimal control:
+ ```math
+ u(t) = 1 \quad \text{on} \quad [t_0, t_s] \\
+ u(t) = -1 \quad \text{on} \quad [t_s, 0]
+ ```
+
+ Now we integrate the system:
+
+ On ( t in [t_0, t_s] ) :
+ ```math
+ x_2' = u = 1 \quad \Rightarrow \quad x_2(t) = t - t_0 \\
+ x_1' = x_2 \quad \Rightarrow \quad x_1(t) = \frac{(t - t_0)^2}{2}
+ ```
+
+ At switching time ( t = t_s ) :
+ ```math
+ \dot{x}_2(t_s) = t_s - t_0 \\
+ \dot{x}_1(t_s) = \frac{(t_s - t_0)^2}{2}
+ ```
+
+ when ( t in [t_s, 0] ) :
+ ```math
+ \dot{x}_2(t) = u(t) = -1 \quad \Rightarrow \quad x_2(t) = x_2(t_s) - (t - t_s) \\
+ \dot{x}_1(t) = x_2(t) \quad \Rightarrow \quad x_1(t) = x_1(t_s) + \int_{t_s}^t x_2(s) ds
+ ```
+
+ Final velocity condition:
+ ```math
+ x_2(0) = 0 \quad \Rightarrow \quad t_s - t_0 + t_s = 0 \quad \Rightarrow \quad t_0 = 2 t_s
+ ```
+
+ Final position:
+ ```math
+ x_1(0) = x_1(t_s) + \frac{t_s^2}{2} \quad \Rightarrow \quad x_1(0) = t_s^2 = 1 \quad \Rightarrow \quad t_s = -1
+ ```
+
+ We deduce:
+ ```math
+ t_0 = 2 * t_s = -2
+ ```
+
+ ### Final solution:
+ - Switching time:
+ ```math
+ t_s = -1
+ ```
+ - Initial time:
+ ```math
+ t_0 = -2
+ ```
+
+ Control:
+ ```math
+ u(t) = 1 \quad \text{on} \quad [-2, -1] \\
+ u(t) = -1 \quad \text{on} \quad [-1, 0]
+ ```
+ ```@raw html
+
+ ```
+## Indirect method
+
+
+Define the pseudo-Hamiltonian and switching function
+```@example initial_time
+
+H(x, p, u) = p[1]*x[2] + p[2]*u
+
+# Hamiltonian lift for switching condition
+F1(x) = [0, 1]
+H1 = Lift(F1)
+
+# Define flows for u = +1 and u = -1
+const u_pos = 1
+const u_neg = -1
+
+
+f_pos = Flow(ocp, (x, p, t0) -> u_pos)
+f_neg = Flow(ocp, (x, p, t0) -> u_neg)
+```
+
+Shooting function adapted for free initial time
+```@example initial_time
+
+function shoot!(s, p0, t1, t0)
+ # Initial conditions
+ x0 = [0, 0] # fixed initial state
+ tf = 0 # fixed final time
+
+ # The flows (time runs from t0 to tf = 0)
+ x1, p1 = f_pos(t0, x0, p0, t1) # from t0 to t1
+ xf, pf = f_neg(t1, x1, p1, tf) # from t1 to tf = 0
+
+ # Final control
+ uf = -1
+
+ # Target final state
+ xf_target = [1, 0]
+
+ # Shooting conditions
+ s[1:2] = xf .- xf_target # reach the target at tf = 0
+ s[3] = H(xf, pf, uf) - 1 # final condition on pseudo-Hamiltonian
+ s[4] = H1(x1, p1) # switching condition at t1
+end
+```
+
+Get initial guess from direct solution
+```@example initial_time
+
+t = time_grid(sol)
+x = state(sol)
+p = costate(sol)
+φ(t) = H1(x(t), p(t)) # switching function
+
+p0 = p(t[1]) # initial costate (at t0)
+t1 = t[argmin(abs.(φ.(t)))] # switching time
+t0 = t[1] # initial time (negative value)
+
+println("p0 = ", p0)
+println("t1 = ", t1)
+println("t0 = ", t0)
+
+#Test shooting function at initial guess
+s = zeros(4)
+shoot!(s, p0, t1, t0)
+println("\nNorm of the shooting function: ‖s‖ = ", norm(s), "\n")
+
+# Solve the nonlinear system
+using NonlinearSolve
+
+# Aggregated function
+nle! = (s, ξ, λ) -> shoot!(s, ξ[1:2], ξ[3], ξ[4])
+
+# Initial guess: [p0[1], p0[2], t1, t0]
+ξ_guess = [p0..., t1, t0]
+prob = NonlinearProblem(nle!, ξ_guess)
+
+# Solve the problem
+indirect_sol = solve(prob; show_trace=Val(true), abstol=1e-8, reltol=1e-8)
+# Extract solution
+p0_opt = indirect_sol.u[1:2]
+t1_opt = indirect_sol.u[3]
+t0_opt = indirect_sol.u[4]
+```
+We can now plot and compare with the direct method.# Compute and plot the optimal trajectory
+
+```@example initial_time
+
+x0 = [0, 0]
+tf = 0
+
+# Concatenate flows: u=+1 from t0 to t1, then u=-1 from t1 to tf
+f = f_pos * (t1_opt, f_neg)
+flow_sol = f((t0_opt, tf), x0, p0_opt; saveat=range(t0_opt, tf, 100))
+
+# Plot comparison
+plt = plot(sol; label="direct", size=(800, 800))
+plot!(plt, flow_sol; label="indirect", color=:red)
+```
\ No newline at end of file
diff --git a/docs/src/tutorial-free-times-orbital.md b/docs/src/tutorial-free-times-orbital.md
new file mode 100644
index 0000000..4a59643
--- /dev/null
+++ b/docs/src/tutorial-free-times-orbital.md
@@ -0,0 +1,280 @@
+```@meta
+Draft = false
+```
+
+# [Optimal control problem with free final orbital time](@id tutorial-free-times-orbital)
+
+## A more concrete example about the change of orbit of a satellite:
+
+## Here are the required packages for the tutorial:
+
+```@example orbit
+using OptimalControl
+using NLPModelsIpopt
+using Plots
+using Printf
+using LinearAlgebra
+using NLsolve
+```
+## Definition of the problem
+
+```@example orbit
+
+const x₁₀ = -42272.67 # initial position x
+const x₂₀ = 0 # initial position y
+const x₃₀ = 0 # initial velocity in x
+const x₄₀ = -5696.72 # initial velocity in y
+const μ = 5.1658620912*1e12 # gravitational parameter
+const γ_max = 0.05 # maximal thrust norm
+const r_f = 42165 # target orbit radius (final distance to origin)
+const rf3 = r_f^3
+const α = sqrt(μ/rf3);
+
+function min_orbit_tf()
+ @def ocp begin
+ tf ∈ R, variable
+ t ∈ [0, tf], time
+ x ∈ R⁴, state
+ u ∈ R², control
+ u₁(t)^2 + u₂(t)^2 ≤ γ_max^2
+ x(0) == [x₁₀, x₂₀, x₃₀, x₄₀]
+ x₁(tf)^2 + x₂(tf)^2 == r_f^2
+ 0.05 ≤ tf ≤ Inf
+ ẋ(t) == [
+ x₃(t),
+ x₄(t),
+ -μ * x₁(t) / ((x₁(t)^2 + x₂(t)^2)^(3/2)) + u₁(t),
+ -μ * x₂(t) / ((x₁(t)^2 + x₂(t)^2)^(3/2)) + u₂(t)
+ ]
+ tf → min
+ end
+
+ return ocp
+end
+nothing # hide
+```
+
+# Direct resolution with minimal orbital time :
+
+We now solve the problem using a direct method, with automatic treatment of the free initial time.
+
+```@example orbit
+ocp = min_orbit_tf()
+sol = solve(ocp;init=(variable=13.4,), grid_size=100)
+```
+And plot the solution.
+
+```@example orbit
+plot(sol; label="direct", size=(800, 800))
+```
+# Indirect resolution with minimal orbital time :
+
+
+## Hamiltonian and optimal control
+
+From the Pontryagin Maximum Principle, the optimal control is bang-bang:
+
+```@example orbit
+ocp = min_orbit_tf()
+
+function optimal_control(x, p, tf)
+ p₃, p₄ = p[3], p[4]
+ p_norm = sqrt(p₃^2 + p₄^2)
+
+ if p_norm > 1e-10
+ u₁ = γ_max * p₃ / p_norm
+ u₂ = γ_max * p₄ / p_norm
+ return [u₁, u₂]
+ else
+ return [0.0, 0.0]
+ end
+end
+
+function hamiltonian(x, p, u)
+ x₁, x₂, x₃, x₄ = x
+ p₁, p₂, p₃, p₄ = p
+ u₁, u₂ = u
+
+ r = sqrt(x₁^2 + x₂^2)
+ r³ = r^3
+
+ return (p₁*x₃ + p₂*x₄ +
+ p₃*(-μ*x₁/r³ + u₁) +
+ p₄*(-μ*x₂/r³ + u₂))
+end
+# Create flow using OptimalControl.jl
+flow = Flow(ocp, optimal_control)
+```
+
+
+## Shooting function
+
+The shooting function encodes the boundary conditions:
+
+```@example orbit
+function shooting_function!(s, ξ)
+ p0 = ξ[1:4]
+ tf = ξ[5]
+
+ x0 = [x₁₀, x₂₀, x₃₀, x₄₀]
+
+ try
+ # Use OptimalControl.jl Flow interface with variable
+ x_final, p_final = flow(0.0, x0, p0, tf, tf)
+
+ x₁f, x₂f, x₃f, x₄f = x_final
+ p₁f, p₂f, p₃f, p₄f = p_final
+
+ s[1] = x₁f^2 + x₂f^2 - r_f^2
+ s[2] = p₁f
+ s[3] = p₂f
+
+ u_final = optimal_control(x_final, p_final, tf)
+ H_final = hamiltonian(x_final, p_final, u_final)
+ s[4] = H_final
+
+ s[5] = p₁f^2 + p₂f^2 + p₃f^2 + p₄f^2 - 1.0
+
+ catch e
+ fill!(s, 1e6)
+ end
+
+ return nothing
+end
+```
+
+## Initial guess and solver setup
+
+```@example orbit
+# Initial guess for the shooting variables [p₀, tf]
+
+function get_initial_guess()
+ tf_guess = 13.4
+
+ p0_guess = [1.0323e-4, 4.915e-5, 3.568e-4, -1.554e-4]
+
+ return vcat(p0_guess, tf_guess)
+end
+
+ξ_init = get_initial_guess()
+println("Initial guess: p₀ = $(ξ_init[1:4]), tf = $(ξ_init[5])")
+```
+
+## Solve the shooting problem
+
+```@example orbit
+# Solve the shooting problem using NLsolve
+function solve_indirect_method()
+ println("Solving indirect method...")
+
+ # Set up the shooting problem
+ result = nlsolve(shooting_function!, ξ_init,
+ method=:trust_region,
+ ftol=1e-10,
+ xtol=1e-10,
+ iterations=1000,
+ show_trace=false)
+
+ if result.f_converged || result.x_converged
+ println("Shooting method converged!")
+ println("Final residual norm: $(norm(result.zero))")
+ return result.zero
+ else
+ println("Shooting method failed to converge")
+ println("Final residual norm: $(norm(result.zero))")
+ return nothing
+ end
+end
+
+ξ_solution = solve_indirect_method()
+
+if ξ_solution !== nothing
+ p0_opt = ξ_solution[1:4]
+ tf_opt = ξ_solution[5]
+
+ println("\nOptimal solution:")
+ println("p₀ = $(p0_opt)")
+ println("tf = $(tf_opt)")
+
+ # Verify the solution
+ s_check = zeros(5)
+ shooting_function!(s_check, ξ_solution)
+ println("Shooting function residual: $(s_check)")
+end
+```
+
+## Generate the indirect solution trajectory
+
+```@example orbit
+if ξ_solution !== nothing
+ # Generate the optimal trajectory
+ x0 = [x₁₀, x₂₀, x₃₀, x₄₀]
+ p0_opt = ξ_solution[1:4]
+ tf_opt = ξ_solution[5]
+
+ # Define time points for plotting FIRST
+ t_indirect = range(0, tf_opt, length=1000)
+
+ # Solve the optimal trajectory
+ sol_indirect = flow((0.0, tf_opt), x0, p0_opt)
+
+ # Extract trajectories by evaluating at time points
+ x_traj = zeros(length(t_indirect), 4)
+ p_traj = zeros(length(t_indirect), 4)
+ u_traj = zeros(length(t_indirect), 2)
+
+ for (i, t) in enumerate(t_indirect)
+ x_traj[i, :], p_traj[i, :] = flow(0.0, x0, p0_opt, t)
+ u_traj[i, :] = optimal_control(x_traj[i, :], p_traj[i, :], tf_opt)
+ end
+
+ println("Indirect solution generated successfully!")
+ println("Final time: $(tf_opt)")
+ println("Final position: $(x_traj[end, 1:2])")
+ println("Final radius: $(sqrt(x_traj[end, 1]^2 + x_traj[end, 2]^2))")
+end
+```
+# Visualisation of results
+```@raw html
+
+Click to show/hide indirect method visualization code
+```
+
+```@example orbit
+# Simple visualization - just the basic plots
+if ξ_solution !== nothing
+ p0_opt = ξ_solution[1:4]
+ tf_opt = ξ_solution[5]
+
+ x0 = [x₁₀, x₂₀, x₃₀, x₄₀]
+
+ t_indirect = range(0, tf_opt, length=1000)
+
+ x_traj = zeros(length(t_indirect), 4)
+ p_traj = zeros(length(t_indirect), 4)
+ u_traj = zeros(length(t_indirect), 2)
+
+ # Use flow function directly to evaluate at each time point
+ for (i, t) in enumerate(t_indirect)
+ x_traj[i, :], p_traj[i, :] = flow(0.0, x0, p0_opt, t)
+ u_traj[i, :] = optimal_control(x_traj[i, :], p_traj[i, :], tf_opt)
+ end
+
+ # Combined plot with 3 rows
+ plt = plot(layout=(3,1), size=(800, 900))
+
+ plot!(plt[1], t_indirect, [x_traj[:,1] x_traj[:,2] x_traj[:,3] x_traj[:,4]],
+ label=["x₁" "x₂" "x₃" "x₄"], title="States")
+
+ plot!(plt[2], t_indirect, [p_traj[:,1] p_traj[:,2] p_traj[:,3] p_traj[:,4]],
+ label=["p₁" "p₂" "p₃" "p₄"], title="Costates")
+
+ plot!(plt[3], t_indirect, [u_traj[:,1] u_traj[:,2]],
+ label=["u₁" "u₂"], title="Control")
+
+ plt # This returns the plot for documentation
+end
+```
+```@raw html
+
+```
\ No newline at end of file