+```
+
+We assume that the mass is constant and equal to one, and that there is no friction. The dynamics are given by
+
+```math
+ \dot x_1(t) = x_2(t), \quad \dot x_2(t) = u(t),\quad u(t) \in \R,
+```
+
+which is simply the [double integrator](https://en.wikipedia.org/w/index.php?title=Double_integrator&oldid=1071399674) system. Let us consider a transfer starting at time $t_0 = 0$ and ending at time $t_f = 1$, for which we want to minimise the transfer energy
+
+```math
+ \frac{1}{2}\int_{0}^{1} u^2(t) \, \mathrm{d}t
+```
+
+starting from $x(0) = (-1, 0)$ and aiming to reach the target $x(1) = (0, 0)$.
+
+First, we need to import the [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl) package to define the optimal control problem, [NLPModelsIpopt.jl](https://jso.dev/NLPModelsIpopt.jl) to solve it, and [Plots.jl](https://docs.juliaplots.org) to visualise the solution.
+
+```@example main
+using OptimalControl
+using NLPModelsIpopt
+using Plots
+```
+
+## Optimal control problem
+
+Let us define the problem with the [`@def`](@ref) macro:
+
+```@raw html
+
+```
+
+First, we need to import the [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl) package to define the
+optimal control problem and [NLPModelsIpopt.jl](https://jso.dev/NLPModelsIpopt.jl) to solve it.
+We also need to import the [Plots.jl](https://docs.juliaplots.org) package to plot the solution.
+
+```@example main
+using OptimalControl
+using NLPModelsIpopt
+using Plots
+```
+
+## Optimal control problem
+
+Let us define the problem:
+
+```@raw html
+
+
+```
+
+This prompt generates the syntactically correct code:
+
+```julia
+# Parameters
+h₀ = 1
+v₀ = 0
+m₀ = 1
+g₀ = 1
+Tc = 3.5
+hc = 500
+vc = 620
+mc = 0.6
+
+# Derived constants
+Dc = 0.5 * vc * m₀ / g₀
+mf = mc * m₀
+c = 0.5 * sqrt(g₀ * h₀)
+Tmax = Tc * m₀ * g₀
+
+# Auxiliary functions
+D(h, v) = Dc * v^2 * exp(-hc * (h - h₀) / h₀)
+g(h) = g₀ * (h₀ / h)^2
+
+rocket = @def begin
+ # Variable (free final time)
+ tf ∈ R, variable
+
+ # Time
+ t ∈ [0, tf], time
+
+ # State: (altitude, velocity, mass)
+ x = (h, v, m) ∈ R³, state
+
+ # Control: thrust
+ T ∈ R, control
+
+ # Dynamics
+ ∂(h)(t) == v(t)
+ ∂(v)(t) == (T(t) - D(h(t), v(t)) - m(t) * g(h(t))) / m(t)
+ ∂(m)(t) == -T(t) / c
+
+ # Initial conditions
+ h(0) == h₀
+ v(0) == v₀
+ m(0) == m₀
+
+ # Final condition
+ m(tf) == mf
+
+ # State constraints
+ h(t) ≥ h₀
+ v(t) ≥ v₀
+ mf ≤ m(t) ≤ m₀
+
+ # Control constraint
+ 0 ≤ T(t) ≤ Tmax
+
+ # Variable constraint
+ tf ≥ 0
+
+ # Objective: maximize final altitude h(tf)
+ -h(tf) → min
+end
+```
\ No newline at end of file
diff --git a/docs/src/manual-flow-api.md b/.save/docs/src/manual-flow-api.md
similarity index 100%
rename from docs/src/manual-flow-api.md
rename to .save/docs/src/manual-flow-api.md
diff --git a/.save/docs/src/manual-flow-ocp.md b/.save/docs/src/manual-flow-ocp.md
new file mode 100644
index 000000000..9a6f34d1d
--- /dev/null
+++ b/.save/docs/src/manual-flow-ocp.md
@@ -0,0 +1,587 @@
+# [How to compute flows from optimal control problems](@id manual-flow-ocp)
+
+```@meta
+CollapsedDocStrings = false
+```
+
+In this tutorial, we explain the `Flow` function, in particular to compute flows from an optimal control problem.
+
+## Basic usage
+
+Les us define a basic optimal control problem.
+
+```@example main
+using OptimalControl
+
+t0 = 0
+tf = 1
+x0 = [-1, 0]
+
+ocp = @def begin
+
+ t ∈ [ t0, tf ], time
+ x = (q, v) ∈ R², state
+ u ∈ R, control
+
+ x(t0) == x0
+ x(tf) == [0, 0]
+ ẋ(t) == [v(t), u(t)]
+
+ ∫( 0.5u(t)^2 ) → min
+
+end
+nothing # hide
+```
+
+The **pseudo-Hamiltonian** of this problem is
+
+```math
+ H(x, p, u) = p_q\, v + p_v\, u + p^0 u^2 /2,
+```
+
+where $p^0 = -1$ since we are in the normal case. From the Pontryagin maximum principle, the maximising control is given in feedback form by
+
+```math
+u(x, p) = p_v
+```
+
+since $\partial^2_{uu} H = p^0 = - 1 < 0$.
+
+```@example main
+u(x, p) = p[2]
+nothing # hide
+```
+
+Actually, if $(x, u)$ is a solution of the optimal control problem,
+then, the Pontryagin maximum principle tells us that there exists a costate $p$ such that $u(t) = u(x(t), p(t))$
+and such that the pair $(x, p)$ satisfies:
+
+```math
+\begin{array}{l}
+ \dot{x}(t) = \displaystyle\phantom{-}\nabla_p H(x(t), p(t), u(x(t), p(t))), \\[0.5em]
+ \dot{p}(t) = \displaystyle - \nabla_x H(x(t), p(t), u(x(t), p(t))).
+\end{array}
+```
+
+The `Flow` function aims to compute $t \mapsto (x(t), p(t))$ from the optimal control problem `ocp` and the control in feedback form `u(x, p)`.
+
+!!! note "Nota bene"
+
+ Actually, writing $z = (x, p)$, then the pair $(x, p)$ is also solution of
+
+ ```math
+ \dot{z}(t) = \vec{\mathbf{H}}(z(t)),
+ ```
+ where $\mathbf{H}(z) = H(z, u(z))$ and $\vec{\mathbf{H}} = (\nabla_p \mathbf{H}, -\nabla_x \mathbf{H})$. This is what is actually computed by `Flow`.
+
+Let us try to get the associated flow:
+
+```julia
+julia> f = Flow(ocp, u)
+ERROR: ExtensionError. Please make: julia> using OrdinaryDiffEq
+```
+
+As you can see, an error occurred since we need the package [OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs).
+This package provides numerical integrators to compute solutions of the ordinary differential equation
+$\dot{z}(t) = \vec{\mathbf{H}}(z(t))$.
+
+!!! note "OrdinaryDiffEq.jl"
+
+ The package OrdinaryDiffEq.jl is part of [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs). You can either use one or the other.
+
+```@example main
+using OrdinaryDiffEq
+f = Flow(ocp, u)
+nothing # hide
+```
+
+Now we have the flow of the associated Hamiltonian vector field, we can use it. Some simple calculations shows that the initial covector $p(0)$ solution of the Pontryagin maximum principle is $[12, 6]$. Let us check that integrating the flow from $(t_0, x_0, p_0) = (0, [-1, 0], [12, 6])$ to the final time $t_f$ we reach the target $x_f = [0, 0]$.
+
+```@example main
+p0 = [12, 6]
+xf, pf = f(t0, x0, p0, tf)
+xf
+```
+
+If you prefer to get the state, costate and control trajectories at any time, you can call the flow like this:
+
+```@example main
+sol = f((t0, tf), x0, p0)
+nothing # hide
+```
+
+In this case, you obtain a data that you can plot exactly like when solving the optimal control problem
+with the function [`solve`](@ref). See for instance the [basic example](@ref example-double-integrator-energy-solve-plot) or the
+[plot tutorial](@ref manual-plot).
+
+```@example main
+using Plots
+plot(sol)
+```
+
+You can notice from the graph of `v` that the integrator has made very few steps:
+
+```@example main
+time_grid(sol)
+```
+
+!!! note "Time grid"
+
+ The function [`time_grid`](@ref) returns the discretised time grid returned by the solver. In this case, the solution has been computed by numerical integration with an adaptive step-length Runge-Kutta scheme.
+
+To have a better visualisation (the accuracy won't change), you can provide a fine grid.
+
+```@example main
+sol = f((t0, tf), x0, p0; saveat=range(t0, tf, 100))
+plot(sol)
+```
+
+The argument `saveat` is an option from OrdinaryDiffEq.jl. Please check the
+[list of common options](https://docs.sciml.ai/DiffEqDocs/stable/basics/common_solver_opts/#solver_options).
+For instance, one can change the integrator with the keyword argument `alg` or the absolute tolerance with
+`abstol`. Note that you can set an option when declaring the flow or set an option in a particular call of the flow.
+In the following example, the integrator will be `BS5()` and the absolute tolerance will be `abstol=1e-8`.
+
+```@example main
+f = Flow(ocp, u; alg=BS5(), abstol=1) # alg=BS5(), abstol=1
+xf, pf = f(t0, x0, p0, tf; abstol=1e-8) # alg=BS5(), abstol=1e-8
+```
+
+## Non-autonomous case
+
+Let us consider the following optimal control problem:
+
+```@example main
+t0 = 0
+tf = π/4
+x0 = 0
+xf = tan(π/4) - 2log(√(2)/2)
+
+ocp = @def begin
+
+ t ∈ [t0, tf], time
+ x ∈ R, state
+ u ∈ R, control
+
+ x(t0) == x0
+ x(tf) == xf
+ ẋ(t) == u(t) * (1 + tan(t)) # The dynamics depend explicitly on t
+
+ 0.5∫( u(t)^2 ) → min
+
+end
+nothing # hide
+```
+
+The pseudo-Hamiltonian of this problem is
+
+```math
+ H(t, x, p, u) = p\, u\, (1+\tan\, t) + p^0 u^2 /2,
+```
+
+where $p^0 = -1$ since we are in the normal case. We can notice that the pseudo-Hamiltonian is non-autonomous since it explicitly depends on the time $t$.
+
+```@example main
+is_autonomous(ocp)
+```
+
+From the Pontryagin maximum principle, the maximising control is given in feedback form by
+
+```math
+u(t, x, p) = p\, (1+\tan\, t)
+```
+
+since $\partial^2_{uu} H = p^0 = - 1 < 0$.
+
+```@example main
+u(t, x, p) = p * (1 + tan(t))
+nothing # hide
+```
+
+As before, the `Flow` function aims to compute $(x, p)$ from the optimal control problem `ocp` and the control in feedback form `u(t, x, p)`.
+Since the problem is non-autonomous, we must provide a control law that depends on time.
+
+```@example main
+f = Flow(ocp, u)
+nothing # hide
+```
+
+Now we have the flow of the associated Hamiltonian vector field, we can use it. Some simple calculations shows that the initial covector $p(0)$ solution of the Pontryagin maximum principle is $1$. Let us check that integrating the flow from $(t_0, x_0) = (0, 0)$ to the final time $t_f = \pi/4$ we reach the target $x_f = \tan(\pi/4) - 2 \log(\sqrt{2}/2)$.
+
+```@example main
+p0 = 1
+xf, pf = f(t0, x0, p0, tf)
+xf - (tan(π/4) - 2log(√(2)/2))
+```
+
+## Variable
+
+Let us consider an optimal control problem with a (decision / optimisation) variable.
+
+```@example main
+t0 = 0
+x0 = 0
+
+ocp = @def begin
+
+ tf ∈ R, variable # the optimisation variable is tf
+ t ∈ [t0, tf], time
+ x ∈ R, state
+ u ∈ R, control
+
+ x(t0) == x0
+ x(tf) == 1
+ ẋ(t) == tf * u(t)
+
+ tf + 0.5∫(u(t)^2) → min
+
+end
+nothing # hide
+```
+
+As you can see, the variable is the final time `tf`. Note that the dynamics depends on `tf`.
+From the Pontryagin maximum principle, the solution is given by:
+
+```@example main
+tf = (3/2)^(1/4)
+p0 = 2tf/3
+nothing # hide
+```
+
+The input arguments of the maximising control are now the state `x`, the costate `p` and the variable `tf`.
+
+```@example main
+u(x, p, tf) = tf * p
+nothing # hide
+```
+
+Let us check that the final condition `x(tf) = 1` is satisfied.
+
+```@example main
+f = Flow(ocp, u)
+xf, pf = f(t0, x0, p0, tf, tf)
+```
+
+The usage of the flow `f` is the following: `f(t0, x0, p0, tf, v)` where `v` is the variable. If one wants
+to compute the state at time `t1 = 0.5`, then, one must write:
+
+```@example main
+t1 = 0.5
+x1, p1 = f(t0, x0, p0, t1, tf)
+```
+
+!!! note "Free times"
+
+ In the particular cases: the initial time `t0` is the only variable, the final time `tf` is the only variable, or the initial and final times `t0` and `tf` are the only variables and are in order `v=(t0, tf)`, the times do not need to be repeated in the call of the flow:
+
+ ```@example main
+ xf, pf = f(t0, x0, p0, tf)
+ ```
+
+Since the variable is the final time, we can make the time-reparameterisation $t = s\, t_f$ to normalise
+the time $s$ in $[0, 1]$.
+
+```@example main
+ocp = @def begin
+
+ tf ∈ R, variable
+ s ∈ [0, 1], time
+ x ∈ R, state
+ u ∈ R, control
+
+ x(0) == 0
+ x(1) == 1
+ ẋ(s) == tf^2 * u(s)
+
+ tf + (0.5*tf)*∫(u(s)^2) → min
+
+end
+
+f = Flow(ocp, u)
+xf, pf = f(0, x0, p0, 1, tf)
+```
+
+Another possibility is to add a new state variable $t_f(s)$. The problem has no variable anymore.
+
+```@example main
+ocp = @def begin
+
+ s ∈ [0, 1], time
+ y = (x, tf) ∈ R², state
+ u ∈ R, control
+
+ x(0) == 0
+ x(1) == 1
+ dx = tf(s)^2 * u(s)
+ dtf = 0 * u(s) # 0
+ ẏ(s) == [dx, dtf]
+
+ tf(1) + 0.5∫(tf(s) * u(s)^2) → min
+
+end
+
+u(y, q) = y[2] * q[1]
+
+f = Flow(ocp, u)
+yf, pf = f(0, [x0, tf], [p0, 0], 1)
+```
+
+!!! danger "Bug"
+
+ Note that in the previous optimal control problem, we have `dtf = 0 * u(s)` instead of `dtf = 0`. The latter does not work.
+
+!!! note "Goddard problem"
+
+ In the [Goddard problem](https://control-toolbox.org/Tutorials.jl/stable/tutorial-goddard.html#tutorial-goddard-structure), you may find other constructions of flows, especially for singular and boundary arcs.
+
+
+## Concatenation of arcs
+
+In this part, we present how to concatenate several flows. Let us consider the following problem.
+
+```@example main
+t0 = 0
+tf = 1
+x0 = -1
+xf = 0
+
+@def ocp begin
+
+ t ∈ [ t0, tf ], time
+ x ∈ R, state
+ u ∈ R, control
+
+ x(t0) == x0
+ x(tf) == xf
+ -1 ≤ u(t) ≤ 1
+ ẋ(t) == -x(t) + u(t)
+
+ ∫( abs(u(t)) ) → min
+
+end
+nothing # hide
+```
+
+From the Pontryagin maximum principle, the optimal control is a concatenation of an off arc ($u=0$) followed by a
+positive bang arc ($u=1$). The initial costate is
+
+```math
+p_0 = \frac{1}{x_0 - (x_f-1) e^{t_f}}
+```
+
+and the switching time is $t_1 = -\ln(p_0)$.
+
+```@example main
+p0 = 1/( x0 - (xf-1) * exp(tf) )
+t1 = -log(p0)
+nothing # hide
+```
+
+Let us define the two flows and the concatenation. Note that the concatenation of two flows is a flow.
+
+```@example main
+f0 = Flow(ocp, (x, p) -> 0) # off arc: u = 0
+f1 = Flow(ocp, (x, p) -> 1) # positive bang arc: u = 1
+
+f = f0 * (t1, f1) # f0 followed by f1 whenever t ≥ t1
+nothing # hide
+```
+
+Now, we can check that the state reach the target.
+
+```@example main
+sol = f((t0, tf), x0, p0)
+plot(sol)
+```
+
+!!! note "Goddard problem"
+
+ In the [Goddard problem](https://control-toolbox.org/Tutorials.jl/stable/tutorial-goddard.html#tutorial-goddard-plot), you may find more complex concatenations.
+
+For the moment, this concatenation is not equivalent to an exact concatenation.
+
+```@example main
+f = Flow(x -> x)
+g = Flow(x -> -x)
+
+x0 = 1
+φ(t) = (f * (t/2, g))(0, x0, t)
+ψ(t) = g(t/2, f(0, x0, t/2), t)
+
+println("φ(t) = ", abs(φ(1)-x0))
+println("ψ(t) = ", abs(ψ(1)-x0))
+
+t = range(1, 5e2, 201)
+
+plt = plot(yaxis=:log, legend=:bottomright, title="Comparison of concatenations", xlabel="t")
+plot!(plt, t, t->abs(φ(t)-x0), label="OptimalControl")
+plot!(plt, t, t->abs(ψ(t)-x0), label="Classical")
+```
+
+## State constraints
+
+We consider an optimal control problem with a state constraints of order 1.[^1]
+
+[^1]: B. Bonnard, L. Faubourg, G. Launay & E. Trélat, Optimal Control With State Constraints And The Space Shuttle Re-entry Problem, J. Dyn. Control Syst., 9 (2003), no. 2, 155–199.
+
+```@example main
+t0 = 0
+tf = 2
+x0 = 1
+xf = 1/2
+lb = 0.1
+
+ocp = @def begin
+
+ t ∈ [t0, tf], time
+ x ∈ R, state
+ u ∈ R, control
+
+ -1 ≤ u(t) ≤ 1
+ x(t0) == x0
+ x(tf) == xf
+ x(t) - lb ≥ 0 # state constraint
+ ẋ(t) == u(t)
+
+ ∫( x(t)^2 ) → min
+
+end
+nothing # hide
+```
+
+The pseudo-Hamiltonian of this problem is
+
+```math
+ H(x, p, u, \mu) = p\, u + p^0 x^2 + \mu\, c(x),
+```
+
+where $ p^0 = -1 $ since we are in the normal case, and where $c(x) = x - l_b$. Along a boundary arc, when $c(x(t)) = 0$, we have $x(t) = l_b$, so $ x(\cdot) $ is constant. Differentiating, we obtain $\dot{x}(t) = u(t) = 0$. Hence, along a boundary arc, the control in feedback form is:
+
+
+```math
+u(x) = 0.
+```
+
+From the maximisation condition, along a boundary arc, we have $p(t) = 0$. Differentiating, we obtain $\dot{p}(t) = 2 x(t) - \mu(t) = 0$. Hence, along a boundary arc, the dual variable $\mu$ is given in feedback form by:
+
+```math
+\mu(x) = 2x.
+```
+
+!!! note
+
+ Within OptimalControl.jl, the constraint must be given in the form:
+ ```julia
+ c([t, ]x, u[, v])
+ ```
+ the control law in feedback form must be given as:
+ ```julia
+ u([t, ]x, p[, v])
+ ```
+ and the dual variable:
+ ```julia
+ μ([t, ]x, p[, v])
+ ```
+ The time `t` must be provided when the problem is [non-autonomous](@ref manual-model-time-dependence) and the variable `v` must be given when the optimal control problem contains a [variable](@ref manual-abstract-variable) to optimise.
+
+The optimal control is a concatenation of 3 arcs: a negative bang arc followed by a boundary arc, followed by a positive bang arc. The initial covector is approximately $p(0)=-0.982237546583301$, the first switching time is $t_1 = 0.9$, and the exit time of the boundary is $t_2 = 1.6$. Let us check this by concatenating the three flows.
+
+```@example main
+u(x) = 0 # boundary control
+c(x) = x-lb # constraint
+μ(x) = 2x # dual variable
+
+f1 = Flow(ocp, (x, p) -> -1)
+f2 = Flow(ocp, (x, p) -> u(x), (x, u) -> c(x), (x, p) -> μ(x))
+f3 = Flow(ocp, (x, p) -> +1)
+
+t1 = 0.9
+t2 = 1.6
+f = f1 * (t1, f2) * (t2, f3)
+
+p0 = -0.982237546583301
+xf, pf = f(t0, x0, p0, tf)
+xf
+```
+
+## Jump on the costate
+
+Let consider the following problem:
+
+```@example main
+t0=0
+tf=1
+x0=[0, 1]
+l = 1/9
+@def ocp begin
+ t ∈ [ t0, tf ], time
+ x ∈ R², state
+ u ∈ R, control
+ x(t0) == x0
+ x(tf) == [0, -1]
+ x₁(t) ≤ l, (x_con)
+ ẋ(t) == [x₂(t), u(t)]
+ 0.5∫(u(t)^2) → min
+end
+nothing # hide
+```
+
+The pseudo-Hamiltonian of this problem is
+
+```math
+ H(x, p, u, \mu) = p_1\, x_2 + p_2\, u + 0.5\, p^0 u^2 + \mu\, c(x),
+```
+
+where $ p^0 = -1 $ since we are in the normal case, and where the constraint is $c(x) = l - x_1 \ge 0$. Along a boundary arc, when $c(x(t)) = 0$, we have $x_1(t) = l$, so $\dot{x}_1(t) = x_2(t) = 0$. Differentiating again, we obtain $\dot{x}_2(t) = u(t) = 0$ (the constraint is of order 2). Hence, along a boundary arc, the control in feedback form is:
+
+
+```math
+u(x, p) = 0.
+```
+
+From the maximisation condition, along a boundary arc, we have $p_2(t) = 0$. Differentiating, we obtain $\dot{p}_2(t) = -p_1(t) = 0$. Differentiating again, we obtain $\dot{p}_1(t) = \mu(t) = 0$. Hence, along a boundary arc, the Lagrange multiplier $\mu$ is given in feedback form by:
+
+```math
+\mu(x, p) = 0.
+```
+
+Outside a boundary arc, the maximisation condition gives $u(x, p) = p_2$. A deeper analysis of the problem shows that the optimal solution has 3 arcs, the first and the third ones are interior to the constraint. The second arc is a boundary arc, that is $x_1(t) = l$ along the second arc. We denote by $t_1$ and $t_2$ the two switching times. We have $t_1 = 3l = 1/3$ and $t_2 = 1 - 3l = 2/3$, since $l=1/9$. The initial costate solution is $p(0) = [-18, -6]$.
+
+!!! danger "Important"
+
+ The costate is discontinuous at $t_1$ and $t_2$ with a jump of $18$.
+
+Let us compute the solution concatenating the flows with the jumps.
+
+```@example main
+t1 = 3l
+t2 = 1 - 3l
+p0 = [-18, -6]
+
+fs = Flow(ocp,
+ (x, p) -> p[2] # control along regular arc
+ )
+fc = Flow(ocp,
+ (x, p) -> 0, # control along boundary arc
+ (x, u) -> l-x[1], # state constraint
+ (x, p) -> 0 # Lagrange multiplier
+ )
+
+ν = 18 # jump value of p1 at t1 and t2
+
+f = fs * (t1, [ν, 0], fc) * (t2, [ν, 0], fs)
+
+xf, pf = f(t0, x0, p0, tf) # xf should be [0, -1]
+```
+
+Let us solve the problem with a direct method to compare with the solution from the flow.
+
+```@example main
+using NLPModelsIpopt
+
+direct_sol = solve(ocp)
+plot(direct_sol; label="direct", size=(800, 700))
+
+flow_sol = f((t0, tf), x0, p0; saveat=range(t0, tf, 100))
+plot!(flow_sol; label="flow", state_style=(color=3,), linestyle=:dash)
+```
\ No newline at end of file
diff --git a/.save/docs/src/manual-flow-others.md b/.save/docs/src/manual-flow-others.md
new file mode 100644
index 000000000..e7d5960ae
--- /dev/null
+++ b/.save/docs/src/manual-flow-others.md
@@ -0,0 +1,116 @@
+# [How to compute Hamiltonian flows and trajectories](@id manual-flow-others)
+
+```@meta
+CollapsedDocStrings = false
+```
+
+In this tutorial, we explain the `Flow` function, in particular to compute flows from a Hamiltonian vector fields, but also from general vector fields.
+
+## Introduction
+
+Consider the simple optimal control problem from the [basic example page](@ref example-double-integrator-energy). The **pseudo-Hamiltonian** is
+
+```math
+ H(x, p, u) = p_q\, v + p_v\, u + p^0 u^2 /2,
+```
+
+where $x=(q,v)$, $p=(p_q,p_v)$, $p^0 = -1$ since we are in the normal case. From the Pontryagin maximum principle, the maximising control is given in feedback form by
+
+```math
+u(x, p) = p_v
+```
+
+since $\partial^2_{uu} H = p^0 = - 1 < 0$.
+
+```@example main
+u(x, p) = p[2]
+nothing # hide
+```
+
+Actually, if $(x, u)$ is a solution of the optimal control problem,
+then, the Pontryagin maximum principle tells us that there exists a costate $p$ such that $u(t) = u(x(t), p(t))$
+and such that the pair $(x, p)$ satisfies:
+
+```math
+\begin{array}{l}
+ \dot{x}(t) = \displaystyle\phantom{-}\nabla_p H(x(t), p(t), u(x(t), p(t))), \\[0.5em]
+ \dot{p}(t) = \displaystyle - \nabla_x H(x(t), p(t), u(x(t), p(t))).
+\end{array}
+```
+
+!!! note "Nota bene"
+
+ Actually, writing $z = (x, p)$, then the pair $(x, p)$ is also solution of
+
+ ```math
+ \dot{z}(t) = \vec{\mathbf{H}}(z(t)),
+ ```
+ where $\mathbf{H}(z) = H(z, u(z))$ and $\vec{\mathbf{H}} = (\nabla_p \mathbf{H}, -\nabla_x \mathbf{H})$.
+
+Let us import the necessary packages.
+
+```@example main
+using OptimalControl
+using OrdinaryDiffEq
+```
+
+The package [OrdinaryDiffEq.jl](https://docs.sciml.ai/DiffEqDocs) provides numerical integrators to compute solutions of ordinary differential equations.
+
+!!! note "OrdinaryDiffEq.jl"
+
+ The package OrdinaryDiffEq.jl is part of [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs). You can either use one or the other.
+
+## Extremals from the Hamiltonian
+
+The pairs $(x, p)$ solution of the Hamitonian vector field are called *extremals*. We can compute some constructing the flow from the optimal control problem and the control in feedback form. Another way to compute extremals is to define explicitly the Hamiltonian.
+
+```@example main
+H(x, p, u) = p[1] * x[2] + p[2] * u - 0.5 * u^2 # pseudo-Hamiltonian
+H(x, p) = H(x, p, u(x, p)) # Hamiltonian
+
+z = Flow(Hamiltonian(H))
+
+t0 = 0
+tf = 1
+x0 = [-1, 0]
+p0 = [12, 6]
+xf, pf = z(t0, x0, p0, tf)
+```
+
+## Extremals from the Hamiltonian vector field
+
+You can also provide the Hamiltonian vector field.
+
+```@example main
+Hv(x, p) = [x[2], p[2]], [0.0, -p[1]] # Hamiltonian vector field
+
+z = Flow(HamiltonianVectorField(Hv))
+xf, pf = z(t0, x0, p0, tf)
+```
+
+Note that if you call the flow on `tspan=(t0, tf)`, then you obtain the output solution
+from OrdinaryDiffEq.jl.
+
+```@example main
+sol = z((t0, tf), x0, p0)
+xf, pf = sol(tf)[1:2], sol(tf)[3:4]
+```
+
+## Trajectories
+
+You can also compute trajectories from the control dynamics $(x, u) \mapsto (v, u)$ and a control law
+$t \mapsto u(t)$.
+
+```@example main
+u(t) = 6-12t
+x = Flow((t, x) -> [x[2], u(t)]; autonomous=false) # the vector field depends on t
+x(t0, x0, tf)
+```
+
+Again, giving a `tspan` you get an output solution from OrdinaryDiffEq.jl.
+
+```@example main
+using Plots
+sol = x((t0, tf), x0)
+plot(sol)
+```
diff --git a/.save/docs/src/manual-initial-guess.md b/.save/docs/src/manual-initial-guess.md
new file mode 100644
index 000000000..06db0c3e6
--- /dev/null
+++ b/.save/docs/src/manual-initial-guess.md
@@ -0,0 +1,190 @@
+# [Initial guess (or iterate) for the resolution](@id manual-initial-guess)
+
+```@meta
+CurrentModule = OptimalControl
+```
+
+We present the different possibilities to provide an initial guess to solve an
+optimal control problem with the [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl) package.
+
+First, we need to import OptimalControl.jl to define the
+optimal control problem and [NLPModelsIpopt.jl](https://jso.dev/NLPModelsIpopt.jl) to solve it.
+We also need to import [Plots.jl](https://docs.juliaplots.org) to plot solutions.
+
+```@example main
+using OptimalControl
+using NLPModelsIpopt
+using Plots
+```
+
+For the illustrations, we define the following optimal control problem.
+
+```@example main
+t0 = 0
+tf = 10
+α = 5
+
+ocp = @def begin
+ t ∈ [t0, tf], time
+ v ∈ R, variable
+ x ∈ R², state
+ u ∈ R, control
+ x(t0) == [ -1, 0 ]
+ x₁(tf) == 0
+ ẋ(t) == [ x₂(t), x₁(t) + α*x₁(t)^2 + u(t) ]
+ x₂(tf)^2 + ∫( 0.5u(t)^2 ) → min
+end
+nothing # hide
+```
+
+## Default initial guess
+
+We first solve the problem without giving an initial guess.
+This will default to initialize all variables to 0.1.
+
+```@example main
+# solve the optimal control problem without initial guess
+sol = solve(ocp; display=false)
+println("Number of iterations: ", iterations(sol))
+nothing # hide
+```
+
+Let us plot the solution of the optimal control problem.
+
+```@example main
+plot(sol; size=(600, 450))
+```
+
+Note that the following formulations are equivalent to not giving an initial guess.
+
+```@example main
+sol = solve(ocp; init=nothing, display=false)
+println("Number of iterations: ", iterations(sol))
+
+sol = solve(ocp; init=(), display=false)
+println("Number of iterations: ", iterations(sol))
+nothing # hide
+```
+!!! tip "Interactions with an optimal control solution"
+
+ To get the number of iterations of the solver, check the [`iterations`](@ref) function.
+
+To reduce the number of iterations and improve the convergence, we can give an initial guess to the solver.
+This initial guess can be built from constant values, interpolated vectors, functions, or existing solutions.
+Except when initializing from a solution, the arguments are to be passed as a named tuple ```init=(state=..., control=..., variable=...)``` whose fields are optional. Missing fields will revert to default initialization (ie constant 0.1).
+
+## Constant initial guess
+
+We first illustrate the constant initial guess, using vectors or scalars according to the dimension.
+
+```@example main
+# solve the optimal control problem with initial guess with constant values
+sol = solve(ocp; init=(state=[-0.2, 0.1], control=-0.2, variable=0.05), display=false)
+println("Number of iterations: ", iterations(sol))
+nothing # hide
+```
+
+Partial initializations are also valid, as shown below. Note the ending comma when a single argument is passed, since it must be a tuple.
+```@example main
+# initialisation only on the state
+sol = solve(ocp; init=(state=[-0.2, 0.1],), display=false)
+println("Number of iterations: ", iterations(sol))
+
+# initialisation only on the control
+sol = solve(ocp; init=(control=-0.2,), display=false)
+println("Number of iterations: ", iterations(sol))
+
+# initialisation only on the state and the variable
+sol = solve(ocp; init=(state=[-0.2, 0.1], variable=0.05), display=false)
+println("Number of iterations: ", iterations(sol))
+nothing # hide
+```
+
+## Functional initial guess
+For the state and control, we can also provide functions of time as initial guess.
+
+```@example main
+# initial guess as functions of time
+x(t) = [ -0.2t, 0.1t ]
+u(t) = -0.2t
+
+sol = solve(ocp; init=(state=x, control=u, variable=0.05), display=false)
+println("Number of iterations: ", iterations(sol))
+nothing # hide
+```
+
+## Vector initial guess (interpolated)
+Initialization can also be provided with vectors / matrices to be interpolated along a given time grid.
+In this case the time steps must be given through an additional argument ```time```, which can be a vector or line/column matrix.
+For the values to be interpolated both matrices and vectors of vectors are allowed, but the shape should be *number of time steps x variable dimension*.
+Simple vectors are also allowed for variables of dimension 1.
+
+```@example main
+# initial guess as vector of points
+t_vec = LinRange(t0,tf,4)
+x_vec = [[0, 0], [-0.1, 0.3], [-0.15,0.4], [-0.3, 0.5]]
+u_vec = [0, -0.8, -0.3, 0]
+
+sol = solve(ocp; init=(time=t_vec, state=x_vec, control=u_vec, variable=0.05), display=false)
+println("Number of iterations: ", iterations(sol))
+nothing # hide
+```
+
+Note: in the free final time case, the given time grid should be consistent with the initial guess provided for the final time (in the optimization variables).
+
+## Mixed initial guess
+
+The constant, functional and vector initializations can be mixed, for instance as
+
+```@example main
+# we can mix constant values with functions of time
+sol = solve(ocp; init=(state=[-0.2, 0.1], control=u, variable=0.05), display=false)
+println("Number of iterations: ", iterations(sol))
+
+# wa can mix every possibility
+sol = solve(ocp; init=(time=t_vec, state=x_vec, control=u, variable=0.05), display=false)
+println("Number of iterations: ", iterations(sol))
+nothing # hide
+```
+
+## Solution as initial guess (warm start)
+
+Finally, we can use an existing solution to provide the initial guess.
+The dimensions of the state, control and optimization variable must coincide.
+This particular feature allows an easy implementation of discrete continuations.
+
+```@example main
+# generate the initial solution
+sol_init = solve(ocp; display=false)
+
+# solve the problem using solution as initial guess
+sol = solve(ocp; init=sol_init, display=false)
+println("Number of iterations: ", iterations(sol))
+nothing # hide
+```
+
+Note that you can also manually pick and choose which data to reuse from a solution, by recovering the
+functions ```state(sol)```, ```control(sol)``` and the values ```variable(sol)```.
+For instance the following formulation is equivalent to the ```init=sol``` one.
+
+```@example main
+# use a previous solution to initialise picking data
+sol = solve(ocp;
+ init = (
+ state = state(sol),
+ control = control(sol),
+ variable = variable(sol)
+ ),
+ display=false)
+println("Number of iterations: ", iterations(sol))
+nothing # hide
+```
+
+!!! tip "Interactions with an optimal control solution"
+
+ Please check [`state`](@ref), [`costate`](@ref), [`control`](@ref) and [`variable`](@ref variable(::Solution)) to get data from the solution. The functions `state`, `costate` and `control` return functions of time and `variable` returns a vector.
+
+## Costate / multipliers
+
+For the moment there is no option to provide an initial guess for the costate / multipliers.
+
diff --git a/.save/docs/src/manual-model.md b/.save/docs/src/manual-model.md
new file mode 100644
index 000000000..d09dbffdd
--- /dev/null
+++ b/.save/docs/src/manual-model.md
@@ -0,0 +1,390 @@
+# [The optimal control problem object: structure and usage](@id manual-model)
+
+```@meta
+CollapsedDocStrings = false
+```
+
+In this manual, we'll first recall the **main functionalities** you can use when working with an optimal control problem (OCP). This includes essential operations like:
+
+* **Solving an OCP**: How to find the optimal solution for your defined problem.
+* **Computing flows from an OCP**: Understanding the dynamics and trajectories derived from the optimal solution.
+* **Printing an OCP**: How to display a summary of your problem's definition.
+
+After covering these core functionalities, we'll delve into the **structure of an OCP**. Since an OCP is structured as a [`Model`](@ref) struct, we'll first explain how to **access its underlying attributes**, such as the problem's dynamics, costs, and constraints. Following this, we'll shift our focus to the **simple properties** inherent to an OCP, learning how to determine aspects like whether the problem:
+
+* **Is autonomous**: Does its dynamics depend explicitly on time?
+* **Has a fixed or free initial/final time**: Is the duration of the control problem predetermined or not?
+
+---
+
+**Content**
+
+- [Main functionalities](@ref manual-model-main-functionalities)
+- [Model struct](@ref manual-model-struct)
+- [Attributes and properties](@ref manual-model-attributes)
+
+---
+
+## [Main functionalities](@id manual-model-main-functionalities)
+
+Let's define a basic optimal control problem.
+
+```@example main
+using OptimalControl
+
+t0 = 0
+tf = 1
+x0 = [-1, 0]
+
+ocp = @def begin
+ t ∈ [ t0, tf ], time
+ x = (q, v) ∈ R², state
+ u ∈ R, control
+ x(t0) == x0
+ x(tf) == [0, 0]
+ ẋ(t) == [v(t), u(t)]
+ 0.5∫( u(t)^2 ) → min
+end
+nothing # hide
+```
+
+To print it, simply:
+
+```@example main
+ocp
+```
+
+We can now solve the problem (for more details, visit the [solve manual](@ref manual-solve)):
+
+```@example main
+using NLPModelsIpopt
+solve(ocp)
+nothing # hide
+```
+
+You can also compute flows (for more details, see the [flow manual](@ref manual-flow-ocp)) from the optimal control problem, providing a control law in feedback form. The **pseudo-Hamiltonian** of this problem is
+
+```math
+ H(x, p, u) = p_q\, v + p_v\, u + p^0 u^2 /2,
+```
+
+where $p^0 = -1$ since we are in the normal case. From the Pontryagin maximum principle, the maximising control is given in feedback form by
+
+```math
+u(x, p) = p_v
+```
+
+since $\partial^2_{uu} H = p^0 = - 1 < 0$.
+
+```@example main
+u = (x, p) -> p[2] # control law in feedback form
+
+using OrdinaryDiffEq # needed to import numerical integrators
+f = Flow(ocp, u) # compute the Hamiltonian flow function
+
+p0 = [12, 6] # initial covector solution
+xf, pf = f(t0, x0, p0, tf) # flow from (x0, p0) at time t0 to tf
+xf # should be (0, 0)
+```
+
+!!! note
+
+ A more advanced feature allows for the discretization of the optimal control problem. From the discretized version, you can obtain a Nonlinear Programming problem (or optimization problem) and solve it using any appropriate NLP solver. For more details, visit the [NLP manipulation tutorial](https://control-toolbox.org/Tutorials.jl/stable/tutorial-nlp.html).
+
+## [Model struct](@id manual-model-struct)
+
+The optimal control problem `ocp` is a [`Model`](@ref) struct.
+
+```@docs; canonical=false
+Model
+```
+
+Each field can be accessed directly (`ocp.times`, etc) or by a getter:
+
+- [`times`](@ref)
+- [`state`](@ref)
+- [`control`](@ref)
+- [`variable`](@ref)
+- [`dynamics`](@ref)
+- `objective`
+- [`constraints`](@ref)
+- [`definition`](@ref)
+- [`get_build_examodel`](@ref)
+
+For instance, we can retrieve the `times` and `definition` values.
+
+```@example main
+times(ocp)
+```
+
+```@example main
+definition(ocp)
+```
+
+!!! note
+
+ We refer to the CTModels documentation for more details about this struct and its fields.
+
+## [Attributes and properties](@id manual-model-attributes)
+
+Numerous attributes can be retrieved. To illustrate this, a more complex optimal control problem is defined.
+
+```@example main
+ocp = @def begin
+ v = (w, tf) ∈ R², variable
+ s ∈ [0, tf], time
+ q = (x, y) ∈ R², state
+ u ∈ R, control
+ 0 ≤ tf ≤ 2, (1)
+ u(s) ≥ 0, (cons_u)
+ x(s) + u(s) ≤ 10, (cons_mixed)
+ w == 0
+ x(0) == -1
+ y(0) - tf == 0, (cons_bound)
+ q(tf) == [0, 0]
+ q̇(s) == [y(s)+w, u(s)]
+ 0.5∫( u(s)^2 ) → min
+end
+nothing # hide
+```
+
+### Control, state and variable
+
+You can access the name of the control, state, and variable, along with the names of their components and their dimensions.
+
+```@example main
+using DataFrames
+data = DataFrame(
+ Data=Vector{Symbol}(),
+ Name=Vector{String}(),
+ Components=Vector{Vector{String}}(),
+ Dimension=Vector{Int}(),
+)
+
+# control
+push!(data,(
+ :control,
+ control_name(ocp),
+ control_components(ocp),
+ control_dimension(ocp),
+))
+
+# state
+push!(data,(
+ :state,
+ state_name(ocp),
+ state_components(ocp),
+ state_dimension(ocp),
+))
+
+# variable
+push!(data,(
+ :variable,
+ variable_name(ocp),
+ variable_components(ocp),
+ variable_dimension(ocp),
+))
+```
+
+!!! note
+
+ The names of the components are used for instance when plotting the solution. See the [plot manual](@ref manual-plot).
+
+### Constraints
+
+You can retrieve labelled constraints with the [`constraint`](@ref) function. The `constraint(ocp, label)` method returns a tuple of the form `(type, f, lb, ub)`.
+The signature of the function `f` depends on the symbol `type`. For `:boundary` and `:variable` constraints, the signature is `f(x0, xf, v)` where `x0` is the initial state, `xf` the final state and `v` the variable. For other constraints, the signature is `f(t, x, u, v)`. Here, `t` represents time, `x` the state, `u` the control, and `v` the variable.
+
+```@example main
+(type, f, lb, ub) = constraint(ocp, :eq1)
+println("type: ", type)
+x0 = [0, 1]
+xf = [2, 3]
+v = [1, 4]
+println("val: ", f(x0, xf, v))
+println("lb: ", lb)
+println("ub: ", ub)
+```
+
+```@example main
+(type, f, lb, ub) = constraint(ocp, :cons_bound)
+println("type: ", type)
+println("val: ", f(x0, xf, v))
+println("lb: ", lb)
+println("ub: ", ub)
+```
+
+```@example main
+(type, f, lb, ub) = constraint(ocp, :cons_u)
+println("type: ", type)
+t = 0
+x = [1, 2]
+u = 3
+println("val: ", f(t, x, u, v))
+println("lb: ", lb)
+println("ub: ", ub)
+```
+
+```@example main
+(type, f, lb, ub) = constraint(ocp, :cons_mixed)
+println("type: ", type)
+println("val: ", f(t, x, u, v))
+println("lb: ", lb)
+println("ub: ", ub)
+```
+
+!!! note
+
+ To get the dual variable (or Lagrange multiplier) associated to the constraint, use the [`dual`](@ref) method.
+
+### Dynamics
+
+The dynamics stored in `ocp` are an [in-place function](https://docs.julialang.org/en/v1/manual/functions/#man-argument-passing) (the first argument is mutated upon call) of the form `f!(dx, t, x, u, v)`. Here, `t` represents time, `x` the state, `u` the control, and `v` the variable, with `dx` being the output value.
+
+```@example main
+f! = dynamics(ocp)
+t = 0
+x = [0., 1]
+u = 2
+v = [1, 4]
+dx = similar(x)
+f!(dx, t, x, u, v)
+dx
+```
+
+### Criterion and objective
+
+The criterion can be `:min` or `:max`.
+
+```@example main
+criterion(ocp)
+```
+
+The objective function is either in Mayer, Lagrange or Bolza form.
+
+- Mayer:
+```math
+g(x(t_0), x(t_f), v) \to \min
+```
+- Lagrange:
+```math
+\int_{t_0}^{t_f} f^0(t, x(t), u(t), v)\, \mathrm{d}t \to \min
+```
+- Bolza:
+```math
+g(x(t_0), x(t_f), v) + \int_{t_0}^{t_f} f^0(t, x(t), u(t), v)\, \mathrm{d}t \to \min
+```
+
+The objective of problem `ocp` is `0.5∫( u(t)^2 ) → min`, hence, in Lagrange form. The signature of the Mayer part of the objective is `g(x0, xf, v)` but in our case, the method `mayer` will return an error.
+
+```@repl main
+g = mayer(ocp)
+```
+
+The signature of the Lagrange part of the objective is `f⁰(t, x, u, v)`.
+
+```@example main
+f⁰ = lagrange(ocp)
+f⁰(t, x, u, v)
+```
+
+To avoid having to capture exceptions, you can check the form of the objective:
+
+```@example main
+println("Mayer: ", has_mayer_cost(ocp))
+println("Lagrange: ", has_lagrange_cost(ocp))
+```
+
+### Times
+
+The time variable is not named `t` but `s` in `ocp`.
+
+```@example main
+time_name(ocp)
+```
+
+The initial time is `0`.
+
+```@example main
+initial_time(ocp)
+```
+
+Since the initial time has the value `0`, its name is `string(0)`.
+
+```@example main
+initial_time_name(ocp)
+```
+
+In contrast, the final time is `tf`, since in `ocp` we have `s ∈ [0, tf]`.
+
+```@example main
+final_time_name(ocp)
+```
+
+To get the value of the final time, since it is part of the variable `v = (w, tf)` of `ocp`, we need to provide a variable to the function `final_time`.
+
+```@example main
+v = [1, 2]
+tf = final_time(ocp, v)
+```
+
+```@repl main
+final_time(ocp)
+```
+
+To check whether the initial or final time is fixed or free (i.e., part of the variable), you can use the following functions:
+
+```@example main
+println("Fixed initial time: ", has_fixed_initial_time(ocp))
+println("Fixed final time: ", has_fixed_final_time(ocp))
+```
+
+Or, similarly:
+
+```@example main
+println("Free initial time: ", has_free_initial_time(ocp))
+println("Free final time: ", has_free_final_time(ocp))
+```
+
+### [Time dependence](@id manual-model-time-dependence)
+
+Optimal control problems can be **autonomous** or **non-autonomous**. In an autonomous problem, neither the dynamics nor the Lagrange cost explicitly depends on the time variable.
+
+The following problem is autonomous.
+
+```@example main
+ocp = @def begin
+ t ∈ [ 0, 1 ], time
+ x ∈ R, state
+ u ∈ R, control
+ ẋ(t) == u(t) # no explicit dependence on t
+ x(1) + 0.5∫( u(t)^2 ) → min # no explicit dependence on t
+end
+is_autonomous(ocp)
+```
+
+The following problem is non-autonomous since the dynamics depends on `t`.
+
+```@example main
+ocp = @def begin
+ t ∈ [ 0, 1 ], time
+ x ∈ R, state
+ u ∈ R, control
+ ẋ(t) == u(t) + t # explicit dependence on t
+ x(1) + 0.5∫( u(t)^2 ) → min
+end
+is_autonomous(ocp)
+```
+
+Finally, this last problem is non-autonomous because the Lagrange part of the cost depends on `t`.
+
+```@example main
+ocp = @def begin
+ t ∈ [ 0, 1 ], time
+ x ∈ R, state
+ u ∈ R, control
+ ẋ(t) == u(t)
+ x(1) + 0.5∫( t + u(t)^2 ) → min # explicit dependence on t
+end
+is_autonomous(ocp)
+```
diff --git a/.save/docs/src/manual-plot.md b/.save/docs/src/manual-plot.md
new file mode 100644
index 000000000..1238496e8
--- /dev/null
+++ b/.save/docs/src/manual-plot.md
@@ -0,0 +1,450 @@
+# [How to plot a solution](@id manual-plot)
+
+```@meta
+CollapsedDocStrings = false
+```
+
+In this tutorial, we explain the different options for plotting the solution of an optimal control problem using the `plot` and `plot!` functions, which are extensions of the [Plots.jl](https://docs.juliaplots.org) package. Use `plot` to create a new plot object, and `plot!` to add to an existing one:
+
+```julia
+plot(args...; kw...) # creates a new Plot, and set it to be the `current`
+plot!(args...; kw...) # modifies Plot `current()`
+plot!(plt, args...; kw...) # modifies Plot `plt`
+```
+
+More precisely, the signature of `plot`, to plot a solution, is as follows.
+
+```@docs; canonical=false
+plot(::CTModels.Solution, ::Symbol...)
+plot!(::CTModels.Solution, ::Symbol...)
+plot!(::Plots.Plot, ::CTModels.Solution, ::Symbol...)
+```
+
+## Argument Overview
+
+The table below summarizes the main plotting arguments and links to the corresponding documentation sections for detailed explanations:
+
+| Section | Relevant Arguments |
+| :---------------------------------------------------| :-------------------------------------------------------------------------------------------- |
+| [Basic concepts](@ref manual-plot-basic) | `size`, `state_style`, `costate_style`, `control_style`, `time_style`, `kwargs...` |
+| [Split vs. group layout](@ref manual-plot-layout) | `layout` |
+| [Plotting control norm](@ref manual-plot-control) | `control` |
+| [Normalised time](@ref manual-plot-time) | `time` |
+| [Constraints](@ref manual-plot-constraints) | `state_bounds_style`, `control_bounds_style`, `path_style`, `path_bounds_style`, `dual_style` |
+| [What to plot](@ref manual-plot-select) | `description...` |
+
+You can plot solutions obtained from the `solve` function or from a flow computed using an optimal control problem and a control law. See the [Basic Concepts](@ref manual-plot-basic) and [From Flow function](@ref manual-plot-flow) sections for details.
+
+To overlay a new plot on an existing one, use the `plot!` function (see [Add a plot](@ref manual-plot-add)).
+
+If you prefer full control over the visualisation, you can extract the state, costate, and control to create your own plots. Refer to the [Custom plot](@ref manual-plot-custom) section for guidance. You can also access the subplots.
+
+## The problem and the solution
+
+Let us start by importing the packages needed to define and solve the problem.
+
+```@example main
+using OptimalControl
+using NLPModelsIpopt
+```
+
+We consider the simple optimal control problem from the [basic example page](@ref example-double-integrator-energy).
+
+```@example main
+t0 = 0 # initial time
+tf = 1 # final time
+x0 = [-1, 0] # initial condition
+xf = [ 0, 0] # final condition
+
+ocp = @def begin
+ t ∈ [t0, tf], time
+ x ∈ R², state
+ u ∈ R, control
+ x(t0) == x0
+ x(tf) == xf
+ ẋ(t) == [x₂(t), u(t)]
+ ∫( 0.5u(t)^2 ) → min
+end
+
+sol = solve(ocp, display=false)
+nothing # hide
+```
+
+## [Basic concepts](@id manual-plot-basic)
+
+The simplest way to plot the solution is to use the `plot` function with the solution as the only argument.
+
+!!! caveat
+
+ The `plot` function for a solution of an optimal control problem extends the `plot` function from Plots.jl. Therefore, you need to import this package in order to plot a solution.
+
+```@example main
+using Plots
+plot(sol)
+```
+
+In the figure above, we have a grid of subplots: the left column displays the state component trajectories, the right column shows the costate component trajectories, and the bottom row contains the control component trajectory.
+
+As in Plots.jl, input data is passed positionally (for example, `sol` in `plot(sol)`), and attributes are passed as keyword arguments (for example, `plot(sol; color = :blue)`). After executing `using Plots` in the REPL, you can use the `plotattr()` function to print a list of all available attributes for series, plots, subplots, or axes.
+
+```julia
+# Valid Operations
+plotattr(:Plot)
+plotattr(:Series)
+plotattr(:Subplot)
+plotattr(:Axis)
+```
+
+Once you have the list of attributes, you can either use the aliases of a specific attribute or inspect a specific attribute to display its aliases and description.
+
+```@repl main
+plotattr("color") # Specific Attribute Example
+```
+
+!!! warning
+
+ Some attributes have different default values in OptimalControl.jl compared to Plots.jl. For instance, the default figure size is 600x400 in Plots.jl, while in OptimalControl.jl, it depends on the number of states and controls.
+
+You can also visit the Plot documentation online to get the descriptions of the attributes:
+
+- To pass attributes to the plot, see the [attributes plot](https://docs.juliaplots.org/latest/generated/attributes_plot/) documentation. For instance, you can specify the size of the figure.
+```@raw html
+