diff --git a/Project.toml b/Project.toml
index 3261df2c..7883d909 100644
--- a/Project.toml
+++ b/Project.toml
@@ -38,6 +38,7 @@ Documenter = "1.17.0"
 ExaModels = "0.9"
 ForwardDiff = "0.10, 1.0"
 LinearAlgebra = "1"
+Literate = "2"
 MadNCL = "0.2"
 MadNLP = "0.9"
 MadNLPGPU = "0.8"
@@ -59,6 +60,7 @@ BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
 DifferentiationInterface = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 MadNCL = "434a0bcb-5a7c-42b2-a9d3-9e3f760e7af0"
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
diff --git a/binder/Project.toml b/binder/Project.toml
new file mode 100644
index 00000000..534c1862
--- /dev/null
+++ b/binder/Project.toml
@@ -0,0 +1,9 @@
+[deps]
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
+MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
+MadNLPGPU = "d72a61cc-809d-412f-99be-fd81f4b8a598"
+NLPModelsIpopt = "f4238b75-b362-5c4c-b852-0801c9a21d71"
+NonlinearSolve = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
+OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/binder/apt.txt b/binder/apt.txt
new file mode 100644
index 00000000..dd4ad866
--- /dev/null
+++ b/binder/apt.txt
@@ -0,0 +1,4 @@
+build-essential
+libgmp-dev
+libblas-dev
+liblapack-dev
diff --git a/binder/postBuild b/binder/postBuild
new file mode 100755
index 00000000..8cd8361e
--- /dev/null
+++ b/binder/postBuild
@@ -0,0 +1,5 @@
+#!/bin/bash
+set -e
+
+# Install OptimalControl.jl from the repo root into the Binder environment
+julia --project=binder -e 'using Pkg; Pkg.develop(path="."); Pkg.precompile()'
diff --git a/docs/Project.toml b/docs/Project.toml
index 0310baf6..e534f88e 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
 CTDirect = "790bbbee-bee9-49ee-8912-a9de031322d5"
 CTFlows = "1c39547c-7794-42f7-af83-d98194f657c2"
@@ -15,6 +16,7 @@ DocumenterMermaid = "a078cd44-4d9c-4618-b545-3ab9d77f9177"
 ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MadNCL = "434a0bcb-5a7c-42b2-a9d3-9e3f760e7af0"
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
 MadNLPGPU = "d72a61cc-809d-412f-99be-fd81f4b8a598"
@@ -28,10 +30,11 @@ UnoSolver = "1baa60ac-02f7-4b39-a7a8-2f4f58486b05"
 
 [compat]
 ADNLPModels = "0.8"
-CTBase = "0.18"
+BenchmarkTools = "1"
+CTBase = "=0.18.8"
 CTDirect = "1"
 CTFlows = "0.8"
-CTModels = "0.10"
+CTModels = "=0.10.1"
 CTParser = "0.8"
 CTSolvers = "0.4"
 CUDA = "5"
@@ -43,6 +46,7 @@ DocumenterMermaid = "0.2"
 ExaModels = "0.9"
 JLD2 = "0.6"
 JSON3 = "1"
+Literate = "2"
 MadNCL = "0.2"
 MadNLP = "0.9"
 MadNLPGPU = "0.8"
diff --git a/docs/make.jl b/docs/make.jl
index 5ad8885f..183c49bb 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -153,7 +153,7 @@ cp(
 Draft = false
 ```
 =#
-draft = false  # Draft mode: if true, @example blocks in markdown are not executed
+draft = false  # Draft mode: skip @example execution globally; tutorial overrides below
 
 # ═══════════════════════════════════════════════════════════════════════════════
 # Load extensions
@@ -201,6 +201,29 @@ function Logging.handle_message(
 end
 global_logger(ExampleSizeThresholdFilter(global_logger()))
 
+# ═══════════════════════════════════════════════════════════════════════════════
+# Literate: generate tutorial.md, tutorial.ipynb, tutorial.jl
+# ═══════════════════════════════════════════════════════════════════════════════
+using Literate
+
+LITERATE_DIR = joinpath(@__DIR__, "src-literate")
+MD_OUTPUT    = joinpath(@__DIR__, "src")
+NB_OUTPUT    = joinpath(@__DIR__, "src", "notebooks")
+JL_OUTPUT    = joinpath(@__DIR__, "src", "scripts")
+mkpath(NB_OUTPUT)
+mkpath(JL_OUTPUT)
+
+for file in ["tutorial.jl"]
+    INPUT = joinpath(LITERATE_DIR, file)
+    # Inject @meta Draft=false so the tutorial executes even with global draft=true
+    function tutorial_postprocess(content)
+        return "```@meta\nDraft = false\n```\n\n" * content
+    end
+    Literate.markdown(INPUT, MD_OUTPUT; documenter=true, postprocess=tutorial_postprocess)
+    Literate.notebook(INPUT, NB_OUTPUT; execute=false)
+    Literate.script(INPUT, JL_OUTPUT)
+end
+
 # ═══════════════════════════════════════════════════════════════════════════════
 # Build documentation
 # ═══════════════════════════════════════════════════════════════════════════════
@@ -229,11 +252,13 @@ with_api_reference(src_dir, ext_dir) do api_pages
                 joinpath("api", "private.md"),
                 joinpath("api", "public.md"),
                 "manual-macro-free.md",
+                "tutorial.md",
             ],
         ),
         pages=[
             "Introduction" => "index.md",
-            "Basic Examples" => [
+            "Guided tour" => "tutorial.md",
+            "Examples" => [
                 "Energy minimisation" => "example-double-integrator-energy.md",
                 "Time mininimisation" => "example-double-integrator-time.md",
                 "Control-free problems" => "example-control-free.md",
diff --git a/docs/src-literate/tutorial.jl b/docs/src-literate/tutorial.jl
new file mode 100644
index 00000000..7f9d9cf6
--- /dev/null
+++ b/docs/src-literate/tutorial.jl
@@ -0,0 +1,474 @@
+# # OptimalControl.jl — a guided tour
+#
+#src ============================================================================
+#src SKELETON — v1 "socle" (~45 min). See .reports/tutorial-brainstorming.md.
+#src Sections form the self-contained socle. Extensions are out of v1.
+#src
+#src Literate conventions used here:
+#src   `# text`   → markdown line (rendered in .md and .ipynb)
+#src   plain code → executed Julia
+#src   `## text`  → stays a code comment inside a code block
+#src   `#src ...`  → source-only note, stripped from BOTH .md and .ipynb (TODOs live here)
+#src   `#md ...`   → markdown output only     `#nb ...` → notebook output only
+#src
+#src LINKS: doc-internal links are split — `#md` emits `@ref`/`@extref` (cross-refs for the
+#src   Documenter build), `#nb` emits plain https URLs (for the Binder notebook). External
+#src   homepages (jso.dev, juliaplots, control-toolbox.org root) stay as single plain URLs.
+#src ============================================================================
+#
+# This tutorial is a guided tour of [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl), part of the [control-toolbox](https://control-toolbox.org) ecosystem. We follow two problems end to end: a simple **double integrator** for modelling, initialisation and the **indirect** (Pontryagin) method, and the **Goddard rocket** for the **direct** method in depth, grid continuation and GPU solving. Advanced topics are linked at the end.
+#
+# It is written for readers with a background in optimal control, ODEs or optimisation. By the end you will be able to define an optimal control problem, solve it by both the direct and indirect methods, and visualise the result — all in a few lines of code.
+#
+#md # !!! note "Run online"
+#md #     You can run this tutorial interactively in your browser — no installation required — by clicking the Binder badge below:
+#md #
+#md #     [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/control-toolbox/OptimalControl.jl/paderborn?urlpath=%2Fdoc%2Ftree%2Fdocs%2Fsrc%2Fnotebooks%2Ftutorial.ipynb)
+
+#src ============================================================================
+# ## The problem, and installing the tools
+#src ============================================================================
+#
+# An **optimal control problem** (OCP) in Bolza form reads
+#
+# ```math
+# J(x, u) = g(x(t_0), x(t_f)) + \int_{t_0}^{t_f} f^{0}(t, x(t), u(t))\,\mathrm{d}t \;\to\; \min,
+# ```
+#
+# subject to the controlled dynamics $\dot{x}(t) = f(t, x(t), u(t))$ and, possibly, box / path / boundary constraints. When $g = 0$ the cost is of **Lagrange** form; when $f^0 = 0$, of **Mayer** form.
+#
+# More generally, the times $t_0$ and $t_f$ may be free (optimisation variables), and a vector $v$ of additional parameters can enter the cost, dynamics and constraints. The full problem then reads
+#
+# ```math
+# \min_{x,u,v}\; g(x(t_0), x(t_f), v) + \int_{t_0}^{t_f} f^{0}(t, x(t), u(t), v)\,\mathrm{d}t,
+# ```
+#
+# subject to $\dot{x}(t) = f(t, x(t), u(t), v)$, box / path / boundary constraints.
+#
+# OptimalControl.jl is the core of the [control-toolbox](https://control-toolbox.org) ecosystem, a modular suite of Julia packages — CTBase (base types & exceptions), CTParser (DSL parsing), CTModels (problem data structures), CTDirect (discretisation & NLP transcription), CTFlows (Hamiltonian flows for indirect methods), and CTSolvers (solver orchestration) — that can also be used individually.
+#
+# Installation is a single package:
+#
+# ```julia
+# import Pkg
+# Pkg.add("OptimalControl")
+# ```
+#
+# We load OptimalControl.jl to model the problem, a solver backend ([NLPModelsIpopt.jl](https://jso.dev/NLPModelsIpopt.jl)), and [Plots.jl](https://docs.juliaplots.org).
+
+using OptimalControl
+using NLPModelsIpopt
+using Plots
+
+#src ============================================================================
+# ## Defining a problem: `@def` vs macro-free
+#src ============================================================================
+#
+# Our running example: a wagon of unit mass on a frictionless rail, state $x = (q, v)$ (position, velocity), acceleration controlled by a force $u$. We start at $(-1, 0)$, must reach $(0, 0)$ at $t_f = 1$, and minimise the transfer energy
+#
+# ```math
+# \frac{1}{2}\int_0^1 u^2(t)\,\mathrm{d}t,
+# ```
+#
+# subject to the dynamics
+#
+# ```math
+# \dot{q}(t) = v(t), \qquad \dot{v}(t) = u(t).
+# ```
+
+t0 = 0; tf = 1; x0 = [-1, 0]; xf = [0, 0];
+
+# ### The `@def` macro
+#
+#md # The [`@def`](@ref manual-abstract-syntax) macro lets us write the problem almost exactly as the mathematics:
+#nb # The [`@def`](https://control-toolbox.org/OptimalControl.jl/stable/manual-abstract.html) macro lets us write the problem almost exactly as the mathematics:
+#
+# Each line of the `@def` block mirrors a piece of the mathematical formulation — time, state, control, dynamics, boundary conditions, then cost — in the same order one would write them on paper. Unicode symbols (`∈`, `R²`, `ẋ`, `∫`, `→`) make the code read like the maths; plain ASCII alternatives (`R^2`, `derivative`, `integral`, `=>`) are available for keyboards or workflows that prefer them.
+
+ocp = @def begin
+    t ∈ [t0, tf], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+
+    x(t0) == x0
+    x(tf) == xf
+
+    ẋ(t) == [v(t), u(t)]
+
+    0.5∫( u(t)^2 ) → min
+end
+#md nothing # hide
+
+# ### The same problem with the macro-free (functional) API
+#
+#md # The [functional API](@ref manual-macro-free) builds the *same* model step by step with plain functions — useful for programmatic problem generation or macro-free library code.
+#nb # The [functional API](https://control-toolbox.org/OptimalControl.jl/stable/manual-macro-free.html) builds the *same* model step by step with plain functions — useful for programmatic problem generation or macro-free library code.
+
+pre = OptimalControl.PreModel()
+
+time!(pre; t0=t0, tf=tf)
+state!(pre, 2, "x", ["q", "v"])
+control!(pre, 1)
+
+function f_energy!(dx, t, x, u, v)
+    dx[1] = x[2]
+    dx[2] = u[1]
+    return nothing
+end
+dynamics!(pre, f_energy!)
+
+function boundary_energy!(b, x0_, xf_, v)
+    b[1] = x0_[1] - x0[1]
+    b[2] = x0_[2] - x0[2]
+    b[3] = xf_[1] - xf[1]
+    b[4] = xf_[2] - xf[2]
+    return nothing
+end
+constraint!(pre, :boundary; f=boundary_energy!, lb=zeros(4), ub=zeros(4), label=:endpoint)
+
+lagrange_energy(t, x, u, v) = 0.5 * u[1]^2
+objective!(pre, :min; lagrange=lagrange_energy)
+
+time_dependence!(pre; autonomous=true)
+
+ocp_func = build(pre)
+#md nothing # hide
+
+# ### What the macro actually does
+#
+# **Key message:** `@def` *translates the expression* into the very same functional calls, and **additionally records the symbolic definition**. We can see the difference directly: the macro keeps the DSL expression, whereas the functional API stores an empty definition.
+
+definition(ocp)          # the macro records the full DSL expression
+
+#-
+
+has_abstract_definition(ocp_func)   # false: functional API stores no abstract definition
+
+#md # !!! warning "Two things to keep in mind"
+#md #     - In the functional API, callbacks are **always vector-valued**: even when the control is scalar, one writes `u[1]` — not `u` — inside `f_energy!` or `lagrange_energy`.
+#md #     - The functional API currently works only with the `:adnlp` modeler; it does **not** support the `:exa` modeler needed for GPU solving — one more reason to prefer `@def` when GPU execution is contemplated (more in the GPU section).
+#nb # **Two things to keep in mind:** (1) in the functional API, callbacks are always vector-valued — even when the control is scalar, one writes `u[1]` — not `u` — inside `f_energy!` or `lagrange_energy`; (2) the functional API currently works only with the `:adnlp` modeler, **not** `:exa` (needed for GPU solving) — one more reason to prefer `@def` when GPU execution is contemplated (more in the GPU section).
+
+#src ============================================================================
+# ## First solve, initial guess, and the costate
+#src ============================================================================
+#
+# Solving is one call, plotting another.
+
+direct_sol = solve(ocp)
+#md nothing # hide
+
+#md #-
+#md direct_sol # hide
+
+#-
+
+plot(direct_sol; size=(800, 600))
+
+# ### The default initial guess
+#
+# With no initial guess, every variable is initialised to `0.1`. We can *see* the initial guess without optimising, by stopping the solver immediately with `max_iter=0`:
+
+sol_init = solve(ocp; init=nothing, max_iter=0, display=false)
+plot(sol_init; size=(800, 600))
+
+#md # !!! note "Notice the right-hand column: the costate is already there"
+#md #     Even though we only ever provide the state, control and (optional) variable, the solver initialises the **adjoint** internally. After optimisation, this right-column costate is exactly the **adjoint $p$ of Pontryagin's Maximum Principle** — the same $p$ we will reuse to start the indirect method in the indirect section. This closes the loop between the direct and indirect methods.
+#nb # **Notice the right-hand column: the costate is already there.** Even though we only ever provide the state, control and (optional) variable, the solver initialises the **adjoint** internally. After optimisation, this right-column costate is exactly the **adjoint $p$ of Pontryagin's Maximum Principle** — the same $p$ we will reuse to start the indirect method in the indirect section. This closes the loop between the direct and indirect methods.
+
+# ### Providing our own initial guess
+#
+# The recommended way to provide an initial guess is the `@init` macro, using the labels from the `@def` block (`q`, `v`, `u` here):
+
+ig = @init ocp begin
+    q(t) := -1 + t
+    v(t) := 0
+    u(t) := 0
+end
+
+sol = solve(ocp; init=ig, display=false)
+println("iterations, default guess: ", iterations(direct_sol))
+println("iterations, @init guess:   ", iterations(sol))
+
+# In this case both guesses give **1 iteration**: the double integrator is a *linear-quadratic* problem, so the NLP is quadratic and Ipopt solves it in a single step regardless of the starting point. Warm-starting only pays off on genuinely nonlinear problems — we will see this with the **Goddard rocket** in the next section.
+
+#md # For all the ways to specify an initial guess, see [Set an initial guess](@ref manual-initial-guess).
+#nb # For all the ways to specify an initial guess, see [Set an initial guess](https://control-toolbox.org/OptimalControl.jl/stable/manual-initial-guess.html).
+#md # !!! note
+#md #     There is currently no way to initialise the costate directly — only state, control and variable can be provided through `@init`. The solver initialises the adjoint internally (as we saw above). Costate initialisation is a planned feature.
+#nb # **Note:** there is currently no way to initialise the costate directly — only state, control and variable can be provided through `@init`. The solver initialises the adjoint internally (as we saw above). Costate initialisation is a planned feature.
+
+#src ============================================================================
+# ## Direct method in depth: Goddard
+#src ============================================================================
+#
+# ### Discretise optimal control problems
+#
+# The **direct** method turns the infinite-dimensional OCP into a finite-dimensional nonlinear program (NLP) by discretising time (Runge–Kutta / collocation) on a grid, then hands the NLP to a solver. It is robust and easy to use.
+#
+# Concretely, time is discretised on a uniform grid $t_0 < t_1 < \dots < t_N = t_f$ with step $h = (t_f - t_0)/N$. The (explicit) Euler scheme, for instance, replaces the dynamics by
+#
+# ```math
+# x_{i} = x_{i-1} + h\,f(t_{i-1}, x_{i-1}, u_{i-1}), \quad i = 1, \dots, N,
+# ```
+#
+# and the integral cost by the corresponding rectangle sum
+#
+# ```math
+# h\sum_{i=0}^{N-1} f^{0}(t_i, x_i, u_i).
+# ```
+#
+# The continuous OCP thus becomes a finite-dimensional NLP in the variables $X = (x_0, \dots, x_N, u_0, \dots, u_N)$, which is passed to an NLP solver such as [Ipopt](https://coin-or.github.io/Ipopt). Higher-order schemes (midpoint, Gauss–Legendre collocation) follow the same principle with different quadrature and interpolation formulas — `solve` defaults to the second-order `:midpoint` scheme, not Euler.
+#
+# ### The Goddard rocket problem
+#
+# To demonstrate convergence behaviour and warm-starting, we need a genuinely nonlinear problem. The **Goddard rocket** — maximise the final altitude, with free final time and a singular arc — is a classic test case.
+
+## Goddard data and dynamics (F0: drift, F1: thrust)
+const r0 = 1
+const v0 = 0
+const m0 = 1
+const mf = 0.6
+const Cd = 310
+const Tmax = 3.5
+const β = 500
+const b = 2
+
+F0(x) = begin
+    r, v, m = x
+    D = Cd * v^2 * exp(-β * (r - 1))
+    [v, -D/m - 1/r^2, 0]
+end
+F1(x) = begin
+    r, v, m = x
+    [0, Tmax/m, -b*Tmax]
+end
+
+goddard = @def begin
+    tf ∈ R, variable
+    t ∈ [t0, tf], time
+    x = (r, v, m) ∈ R³, state
+    u ∈ R, control
+
+    x(t0) == [r0, v0, m0]
+    m(tf) == mf
+    0 ≤ u(t) ≤ 1
+    r(t) ≥ r0
+
+    ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
+
+    r(tf) → max
+end
+#md nothing # hide
+
+# ### Choosing a solver is trivial
+#
+#md # `solve` uses the defaults (collocation, ADNLP modeler, Ipopt, CPU). Switching solver is just loading a package and passing a token (see [Solve a problem](@ref manual-solve)):
+#nb # `solve` uses the defaults (collocation, ADNLP modeler, Ipopt, CPU). Switching solver is just loading a package and passing a token (see [Solve a problem](https://control-toolbox.org/OptimalControl.jl/stable/manual-solve.html)):
+
+using MadNLP
+
+sol_ipopt  = solve(goddard;          grid_size=250, display=false)
+sol_madnlp = solve(goddard, :madnlp; grid_size=250, display=false)
+
+println("Ipopt  : r(tf) = ", objective(sol_ipopt),  ", ", iterations(sol_ipopt),  " iters")
+println("MadNLP : r(tf) = ", objective(sol_madnlp), ", ", iterations(sol_madnlp), " iters")
+
+# The available methods and their options can be inspected with `methods()` and `describe(:collocation)`; we will not dwell on them here.
+
+# ### Grid continuation by warm-starting
+#
+# A solution can be passed **directly** as the initial guess of another solve — it is interpolated onto the new grid. This makes discrete continuation trivial and ties back to the initialisation above. On this nonlinear problem it genuinely **pays**: we compare reaching a fine grid of 1000 two ways:
+#
+# 1. **cold start** — solve `grid_size=1000` directly;
+# 2. **cascade** — solve `grid_size=50` first, then `grid_size=1000` warm-started with that solution.
+
+## solutions computed once, reused for iteration counts and the overlay plot
+sol_cold = solve(goddard; grid_size=1000, display=false)
+
+## warm cascade: grid 50 first, then grid 1000 initialised from it
+s50   = solve(goddard; grid_size=50, display=false)
+s1000 = solve(goddard; grid_size=1000, init=s50, display=false)
+
+println("cold    grid 1000        : ", iterations(sol_cold), " iters")
+println("cascade grid 50 (warm-up): ", iterations(s50),      " iters")
+println("cascade grid 1000 (warm) : ", iterations(s1000),    " iters")
+
+# **Message:** what matters is the iteration count *at the expensive grid* — the warm-started `iterations(s1000)` is well below the cold `iterations(sol_cold)`, even though the cheap `grid_size=50` warm-up adds iterations of its own to the running total; since a grid-50 iteration is far cheaper than a grid-1000 iteration, the cascade still wins on wall-clock time. Overlay the successive solutions to watch convergence:
+
+plt = plot(s50;  label="50", size=(800, 800))
+plot!(plt, s1000; label="1000")
+
+#md # This is grid-refinement warm-starting. The very same mechanism drives **parametric** continuation (homotopy on a physical parameter, e.g. maximum thrust): [Discrete continuation](@extref Tutorials tutorial-continuation).
+#nb # This is grid-refinement warm-starting. The very same mechanism drives **parametric** continuation (homotopy on a physical parameter, e.g. maximum thrust): <https://control-toolbox.org/Tutorials.jl/stable/tutorial-continuation.html>.
+
+# ### Comparison with a bang-bang strategy
+
+# How much better is the optimal solution compared to a naive strategy? We simulate **full thrust until fuel depletion, then coast to apogee** — a bang-bang profile with no optimisation, just two ODE integrations with callbacks.
+
+using OrdinaryDiffEq   # ODE solver (callbacks for the bang-bang simulation)
+
+## Phase 1: u = 1, stop when m = mf (fuel depleted)
+bang1!(dx, x, p, t) = (dx[:] = F0(x) + F1(x))
+cb_fuel = ContinuousCallback((u, t, int) -> u[3] - mf, terminate!)
+sol_bang1 = solve(ODEProblem(bang1!, [r0, v0, m0], (t0, 100.0)), Tsit5(); callback=cb_fuel, reltol=1e-8, abstol=1e-8)
+t1_bang, x1_bang = sol_bang1.t[end], sol_bang1[:, end]
+
+## Phase 2: u = 0, stop when v = 0 (apogee)
+bang2!(dx, x, p, t) = (dx[:] = F0(x))
+cb_apogee = ContinuousCallback((u, t, int) -> u[2], terminate!)
+sol_bang2 = solve(ODEProblem(bang2!, x1_bang, (t1_bang, 1000.0)), Tsit5(); callback=cb_apogee, reltol=1e-8, abstol=1e-8)
+tf_bang, rf_bang = sol_bang2.t[end], sol_bang2[1, end]
+
+println("Bang-bang: r(tf) = ", round(rf_bang, digits=6), "  (t1=", round(t1_bang, digits=4), ", tf=", round(tf_bang, digits=4), ")")
+println("Optimal:   r(tf) = ", round(objective(sol_cold), digits=6), "  (           tf=", round(variable(sol_cold), digits=4), ")")
+
+# The optimal thrust profile uses a **singular arc** — it does not simply push at the maximum. Overlaying the two trajectories on the altitude–velocity plane makes the difference visible:
+
+## assemble the bang-bang trajectory as (t, r, v, m) for plotting
+t_bang = [sol_bang1.t; sol_bang2.t]
+r_bang = [sol_bang1[1, :]; sol_bang2[1, :]]
+
+plt_bang = plot(sol_cold; label="optimal", linewidth=2, color=1)
+plot!(plt_bang[1], t_bang, r_bang; label="bang-bang", linestyle=:dash, linewidth=2, color=2)
+plot(plt_bang[1]; legend=:bottomright, xlabel="time", ylabel="altitude")
+
+#src ============================================================================
+# ## Solving on a GPU
+#src ============================================================================
+#
+# Moving to the GPU is a single token, `:gpu`, which auto-completes to `(:collocation, :exa, :madnlp, :gpu)`. It requires the `:exa` modeler (hence `@def`, not the macro-free API — cf. the definition section) plus a CUDA-capable GPU.
+#
+# In a seminar or on Binder there is usually **no functional GPU**, so the call is *expected to fail* — that is the pedagogical point: the `:gpu` token needs a specific setup. We wrap it in a `try/catch` so the tour keeps running and shows the raised exception.
+
+using MadNLPGPU
+using CUDA
+
+try
+    global sol_gpu = solve(goddard, :gpu; grid_size=1000, display=false)
+    println("GPU solve succeeded — a functional GPU is available.")
+catch e
+    println("GPU solve failed, as expected without a functional GPU.")
+    println("CUDA.functional() = ", CUDA.functional())
+    println("Exception: ", first(sprint(showerror, e), 400))
+end
+
+#md # For the full GPU setup, see [Solve on GPU](@ref manual-solve-gpu).
+#nb # For the full GPU setup, see [Solve on GPU](https://control-toolbox.org/OptimalControl.jl/stable/manual-solve-gpu.html).
+
+#src ============================================================================
+# ## The indirect method
+#src ============================================================================
+#
+# We now return to the **double integrator** `ocp` from the earlier sections. Its shooting has just two unknowns and is initialised by the direct costate above, which makes it ideal to *see* the indirect method. (The Goddard shooting is a *structured multi-arc* problem — see the links in the last section.)
+#
+# In control-toolbox we systematically pair the direct method with the **indirect** one, based on Pontryagin's Maximum Principle (PMP), with pseudo-Hamiltonian
+#
+# ```math
+# H(x,p,u) = p\,f(x,u) + p^0 f^0(x,u) \qquad (\text{normal case } p^0 = -1).
+# ```
+#
+# The PMP gives the maximising control in feedback form
+#
+# ```math
+# u(x,p) = \arg\max_u H,
+# ```
+#
+# and the optimal trajectory solves a boundary value problem that we recast as a **shooting equation**
+#
+# ```math
+# S(p_0) = 0.
+# ```
+#
+# The indirect method proceeds in three steps:
+#
+# 1. **Maximising control.** The PMP yields the control in feedback form $u(x, p) = \arg\max_u H(x, p, u)$. Substituting back gives the maximised Hamiltonian
+#
+#    ```math
+#    \mathbf{H}(x, p) = H(x, p, u(x, p)).
+#    ```
+#
+# 2. **Boundary value problem.** The optimal trajectory satisfies the Hamiltonian system
+#
+#    ```math
+#    \dot{x} = \nabla_p \mathbf{H}, \qquad \dot{p} = -\nabla_x \mathbf{H},
+#    ```
+#
+#    with boundary conditions $x(t_0) = x_0$, $x(t_f) = x_f$.
+#
+# 3. **Shooting function.** Let $\varphi_{t_0, x_0, p_0}(\cdot)$ denote the flow of the Hamiltonian vector field from $(x_0, p_0)$. The shooting function
+#
+#    ```math
+#    S(p_0) = \pi(\varphi_{t_0, x_0, p_0}(t_f)) - x_f, \qquad \pi(x, p) = x,
+#    ```
+#
+#    measures the miss at $t_f$: solving the BVP reduces to finding $p_0$ such that $S(p_0) = 0$.
+#
+# For the energy problem, $H = p_1 v + p_2 u - u^2/2$, so the maximiser is $u = p_2$.
+
+using OrdinaryDiffEq   # ODE solver (Hamiltonian flow)
+using NonlinearSolve   # nonlinear equations (shooting)
+
+## maximising control in feedback form
+u_max(x, p) = p[2]
+
+## Hamiltonian flow of the OCP
+φ = Flow(ocp, u_max);
+
+## state projection π(x, p) = x
+proj((x, p)) = x
+
+## shooting function
+S(p0) = proj(φ(t0, x0, p0, tf)) - xf
+#md nothing # hide
+
+# **The shooting is initialised with the costate of the direct solution** — the very adjoint we highlighted above:
+
+nle!(s, p0, _) = (s[:] = S(p0))
+
+p_of_t   = costate(direct_sol)     # costate as a function of time
+p0_guess = p_of_t(t0)              # initial costate from the direct method
+
+prob = NonlinearProblem(nle!, p0_guess)
+shooting_sol = solve(prob; show_trace=Val(true))
+p0_sol = shooting_sol.u
+
+println("costate p0 = ", p0_sol)
+println("shoot S(p0) = ", S(p0_sol))
+
+# Reconstruct the indirect solution from the flow and overlay it with the direct solution:
+
+indirect_sol = φ((t0, tf), x0, p0_sol; saveat=range(t0, tf, 100))
+
+plt_compare = plot(direct_sol; label="direct", size=(800, 600))
+plot!(plt_compare, indirect_sol; label="indirect")
+
+#md # See [Compute flows from optimal control problems](@ref manual-flow-ocp) for the flow construction, and the [indirect simple shooting tutorial](@extref tutorial-indirect-simple-shooting).
+#nb # See [Compute flows from optimal control problems](https://control-toolbox.org/OptimalControl.jl/stable/manual-flow-ocp.html) for the flow construction, and the [indirect simple shooting tutorial](https://control-toolbox.org/Tutorials.jl/stable/).
+
+#src ============================================================================
+# ## Going further
+#src ============================================================================
+#
+# **Variables & parameters.** Beyond the control, one can optimise **parameters** naturally, both in an OCP (the `variable` keyword of the DSL) and in a differential-constraint optimisation problem **without any control** (a *control-free* problem).
+#md # See [control-free problems](@ref example-control-free).
+#nb # See [control-free problems](https://control-toolbox.org/OptimalControl.jl/stable/example-control-free.html).
+#src NOTE(v1): parameter estimation is only *mentioned* here. The worked example is an
+#src   extension (out of v1) — see .reports/tutorial-brainstorming.md "Extensions futures".
+#
+# **Advanced examples** (each does both direct and indirect):
+#
+#md # - Singular control (control-affine systems) — [singular control](@ref example-singular-control)
+#md # - State constraint — [state constraint](@ref example-state-constraint)
+#md # - Goddard problem — free final time, a singular arc, a state constraint and a structured shooting all at once — [Goddard tutorial](@extref Tutorials tutorial-goddard)
+#nb # - Singular control (control-affine systems) — <https://control-toolbox.org/OptimalControl.jl/stable/example-singular-control.html>
+#nb # - State constraint — <https://control-toolbox.org/OptimalControl.jl/stable/example-state-constraint.html>
+#nb # - Goddard problem — free final time, a singular arc, a state constraint and a structured shooting all at once — <https://control-toolbox.org/Tutorials.jl/stable/tutorial-goddard.html>
+#
+#md # **Discrete continuation** — warm-starting across a family of problems (homotopy on a physical parameter), the grown-up version of the grid continuation above: [Discrete continuation](@extref Tutorials tutorial-continuation).
+#nb # **Discrete continuation** — warm-starting across a family of problems (homotopy on a physical parameter), the grown-up version of the grid continuation above: <https://control-toolbox.org/Tutorials.jl/stable/tutorial-continuation.html>.
+#
+#src ============================================================================
+#src END OF SOCLE (v1). Extensions (parameter estimation worked out, singular, state
+#src constraint, Goddard) go below in later versions — see brainstorming doc.
+#src ============================================================================
diff --git a/docs/src-literate/tutorial_pre.jl b/docs/src-literate/tutorial_pre.jl
new file mode 100644
index 00000000..9bcc0400
--- /dev/null
+++ b/docs/src-literate/tutorial_pre.jl
@@ -0,0 +1,2 @@
+pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))
+pushfirst!(LOAD_PATH, joinpath(@__DIR__, "..", ".."))
diff --git a/docs/src/assets/Manifest.toml b/docs/src/assets/Manifest.toml
index b74944f0..6cd86e40 100644
--- a/docs/src/assets/Manifest.toml
+++ b/docs/src/assets/Manifest.toml
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.12.5"
+julia_version = "1.12.6"
 manifest_format = "2.0"
-project_hash = "505b0cb4420cde5e120adb378ac9d10f3e75d872"
+project_hash = "6e7706d496f1f0c71e587a20bcffec2e751daac3"
 
 [[deps.ADNLPModels]]
 deps = ["ADTypes", "ForwardDiff", "LinearAlgebra", "NLPModels", "Requires", "ReverseDiff", "SparseArrays", "SparseConnectivityTracer", "SparseMatrixColorings"]
@@ -11,9 +11,9 @@ uuid = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a"
 version = "0.8.13"
 
 [[deps.ADTypes]]
-git-tree-sha1 = "f7304359109c768cf32dc5fa2d371565bb63b68a"
+git-tree-sha1 = "d9aaef7c63466eee4de23b4d9dad03629df54bea"
 uuid = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
-version = "1.21.0"
+version = "1.22.1"
 weakdeps = ["ChainRulesCore", "ConstructionBase", "EnzymeCore"]
 
     [deps.ADTypes.extensions]
@@ -56,9 +56,9 @@ version = "0.4.5"
 
 [[deps.Accessors]]
 deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "MacroTools"]
-git-tree-sha1 = "2eeb2c9bef11013efc6f8f97f32ee59b146b09fb"
+git-tree-sha1 = "7063ad1083578215c7c4bf410368150abe8d5524"
 uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
-version = "0.1.44"
+version = "0.1.45"
 
     [deps.Accessors.extensions]
     AxisKeysExt = "AxisKeys"
@@ -79,10 +79,10 @@ version = "0.1.44"
     Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [[deps.Adapt]]
-deps = ["LinearAlgebra", "Requires"]
-git-tree-sha1 = "0761717147821d696c9470a7a86364b2fbd22fd8"
+deps = ["LinearAlgebra"]
+git-tree-sha1 = "daa72978cd7a624246e894a4f4f067706d4e17e2"
 uuid = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
-version = "4.5.2"
+version = "4.7.0"
 weakdeps = ["SparseArrays", "StaticArrays"]
 
     [deps.Adapt.extensions]
@@ -101,9 +101,9 @@ version = "1.1.2"
 
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra"]
-git-tree-sha1 = "54f895554d05c83e3dd59f6a396671dae8999573"
+git-tree-sha1 = "75757da5d9f771ef5909fc84f81d2f9d24127315"
 uuid = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
-version = "7.24.0"
+version = "7.27.0"
 
     [deps.ArrayInterface.extensions]
     ArrayInterfaceAMDGPUExt = "AMDGPU"
@@ -113,6 +113,7 @@ version = "7.24.0"
     ArrayInterfaceCUDSSExt = ["CUDSS", "CUDA"]
     ArrayInterfaceChainRulesCoreExt = "ChainRulesCore"
     ArrayInterfaceChainRulesExt = "ChainRules"
+    ArrayInterfaceFillArraysExt = "FillArrays"
     ArrayInterfaceGPUArraysCoreExt = "GPUArraysCore"
     ArrayInterfaceMetalExt = "Metal"
     ArrayInterfaceReverseDiffExt = "ReverseDiff"
@@ -128,6 +129,7 @@ version = "7.24.0"
     CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
     ChainRules = "082447d4-558c-5d27-93f4-14fc19e9eca2"
     ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
     Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
@@ -167,10 +169,16 @@ version = "0.6.1"
 uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 version = "1.11.0"
 
+[[deps.BenchmarkTools]]
+deps = ["Compat", "JSON", "Logging", "PrecompileTools", "Printf", "Profile", "Statistics", "UUIDs"]
+git-tree-sha1 = "9670d3febc2b6da60a0ae57846ba74670290653f"
+uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
+version = "1.8.0"
+
 [[deps.BitFlags]]
-git-tree-sha1 = "0691e34b3bb8be9307330f88d1a3c3f25466c24d"
+git-tree-sha1 = "bbe1079eecf9c9fbb52765193ad2bae27ae09bc8"
 uuid = "d1d4a3ce-64b1-5f1a-9ba4-7e7e69966f35"
-version = "0.1.9"
+version = "0.1.10"
 
 [[deps.BitTwiddlingConvenienceFunctions]]
 deps = ["Static"]
@@ -208,9 +216,9 @@ version = "0.2.7"
 
 [[deps.CTBase]]
 deps = ["DocStringExtensions"]
-git-tree-sha1 = "2f689cdb1c28d452f143f2c76897cd5273e89cf7"
+git-tree-sha1 = "d0b5663ca6343c475bb9fc54520c9a4a34425c30"
 uuid = "54762871-cc72-4466-b8e8-f6c8b58076cd"
-version = "0.18.7"
+version = "0.18.8"
 
     [deps.CTBase.extensions]
     CoveragePostprocessing = ["Coverage"]
@@ -226,9 +234,9 @@ version = "0.18.7"
 
 [[deps.CTDirect]]
 deps = ["ADNLPModels", "CTModels", "CTSolvers", "DocStringExtensions", "ExaModels", "SolverCore", "SparseArrays", "SparseConnectivityTracer"]
-git-tree-sha1 = "21b9ea456769176849cdcbd20f04a017526daa76"
+git-tree-sha1 = "be35471ebed785f222afaaecad16ba78ba626a6d"
 uuid = "790bbbee-bee9-49ee-8912-a9de031322d5"
-version = "1.0.10"
+version = "1.0.12"
 
 [[deps.CTFlows]]
 deps = ["CTBase", "CTModels", "DocStringExtensions", "ForwardDiff", "LinearAlgebra", "MLStyle", "MacroTools"]
@@ -242,9 +250,9 @@ weakdeps = ["OrdinaryDiffEq"]
 
 [[deps.CTModels]]
 deps = ["CTBase", "DocStringExtensions", "LinearAlgebra", "MLStyle", "MacroTools", "OrderedCollections", "Parameters", "RecipesBase"]
-git-tree-sha1 = "a183ce43a9f912a4b62a8af88ba1415c6e2b1521"
+git-tree-sha1 = "77f34892988d4a5f7e0ca1fe6e0c21255a9e4478"
 uuid = "34c4fa32-2049-4079-8329-de33c2a22e2d"
-version = "0.10.0"
+version = "0.10.1"
 weakdeps = ["JLD2", "JSON3", "Plots"]
 
     [deps.CTModels.extensions]
@@ -254,15 +262,15 @@ weakdeps = ["JLD2", "JSON3", "Plots"]
 
 [[deps.CTParser]]
 deps = ["CTBase", "DocStringExtensions", "MLStyle", "OrderedCollections", "Parameters", "Unicode"]
-git-tree-sha1 = "67c193f4f1e5167123b917bbbdb26947757ed4b0"
+git-tree-sha1 = "690b4ebe0d321779c2f4f867631a5ddb3ba42e09"
 uuid = "32681960-a1b1-40db-9bff-a1ca817385d1"
-version = "0.8.14"
+version = "0.8.15"
 
 [[deps.CTSolvers]]
 deps = ["ADNLPModels", "CTBase", "CTModels", "CommonSolve", "DocStringExtensions", "ExaModels", "KernelAbstractions", "NLPModels", "SolverCore"]
-git-tree-sha1 = "cc69a7ad7e18ca089e0168c6b361ee1e50a3cfc2"
+git-tree-sha1 = "72b7ad3cd1430fc3ec1a6290b607d6f8bd3904b7"
 uuid = "d3e8d392-8e4b-4d9b-8e92-d7d4e3650ef6"
-version = "0.4.16"
+version = "0.4.17"
 
     [deps.CTSolvers.extensions]
     CTSolversCUDA = "CUDA"
@@ -288,9 +296,9 @@ version = "0.4.16"
 
 [[deps.CUDA]]
 deps = ["AbstractFFTs", "Adapt", "BFloat16s", "CEnum", "CUDA_Compiler_jll", "CUDA_Driver_jll", "CUDA_Runtime_Discovery", "CUDA_Runtime_jll", "Crayons", "ExprTools", "GPUArrays", "GPUCompiler", "GPUToolbox", "KernelAbstractions", "LLVM", "LLVMLoopInfo", "LazyArtifacts", "Libdl", "LinearAlgebra", "Logging", "NVTX", "Preferences", "PrettyTables", "Printf", "Random", "Random123", "RandomNumbers", "Reexport", "SparseArrays", "StaticArrays", "Statistics", "demumble_jll"]
-git-tree-sha1 = "ea6a2ab8307059b6c9ea186ff7dfcd032a13b731"
+git-tree-sha1 = "54751d09f9acf05ea7b7ee6baa6a99677c788880"
 uuid = "052768ef-5323-5732-b1bb-66c8b64840ba"
-version = "5.11.0"
+version = "5.11.3"
 
     [deps.CUDA.extensions]
     ChainRulesCoreExt = "ChainRulesCore"
@@ -306,21 +314,21 @@ version = "5.11.0"
 
 [[deps.CUDA_Compiler_jll]]
 deps = ["Artifacts", "CUDA_Driver_jll", "CUDA_Runtime_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
-git-tree-sha1 = "b977706846cb0a75d3842a1fed810ab2e6ab2f94"
+git-tree-sha1 = "c32d22f2f563ce192c88a44b09c2b569f1e7a980"
 uuid = "d1e2174e-dfdc-576e-b43e-73b79eb1aca8"
-version = "0.4.3+0"
+version = "0.4.4+1"
 
 [[deps.CUDA_Driver_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "TOML"]
-git-tree-sha1 = "3b759ec65ac87ad192c2925114fa5c126657a5bd"
+git-tree-sha1 = "2d6222474d868469a72de5bd47c5c25c0e1fe518"
 uuid = "4ee394cb-3365-5eb0-8335-949819d2adfc"
-version = "13.2.1+0"
+version = "13.3.0+0"
 
 [[deps.CUDA_Runtime_Discovery]]
 deps = ["Libdl"]
-git-tree-sha1 = "f9a521f52d236fe49f1028d69e549e7f2644bb72"
+git-tree-sha1 = "79312abe5261a660f94e746e449d2cb2fe3284d9"
 uuid = "1af6417a-86b4-443c-805f-a4643ffb695f"
-version = "1.0.0"
+version = "2.1.0"
 
 [[deps.CUDA_Runtime_jll]]
 deps = ["Artifacts", "CUDA_Driver_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
@@ -329,10 +337,10 @@ uuid = "76a88914-d11a-5bdc-97e0-2f5a05c973a2"
 version = "0.21.0+1"
 
 [[deps.CUDSS]]
-deps = ["CEnum", "CUDA", "CUDA_Runtime_Discovery", "CUDSS_jll", "GPUToolbox", "LinearAlgebra", "SparseArrays"]
-git-tree-sha1 = "715c81db8d627a18d3d1806b068e370dbf724bfc"
+deps = ["CEnum", "CUDA", "CUDSS_jll", "GPUToolbox", "LinearAlgebra", "SparseArrays"]
+git-tree-sha1 = "2d193125a15d15db83d16ab9a02ae4bf7027a18e"
 uuid = "45b445bb-4962-46a0-9369-b4df9d0f772e"
-version = "0.6.7"
+version = "0.6.5"
 
 [[deps.CUDSS_jll]]
 deps = ["Artifacts", "CUDA_Runtime_jll", "CompilerSupportLibraries_jll", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
@@ -342,9 +350,9 @@ version = "0.7.1+0"
 
 [[deps.Cairo_jll]]
 deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "d0efe2c6fdcdaa1c161d206aa8b933788397ec71"
+git-tree-sha1 = "1fa950ebc3e37eccd51c6a8fe1f92f7d86263522"
 uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
-version = "1.18.6+0"
+version = "1.18.7+0"
 
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra"]
@@ -363,9 +371,9 @@ version = "1.0.1"
 
 [[deps.ChunkCodecLibZlib]]
 deps = ["ChunkCodecCore", "Zlib_jll"]
-git-tree-sha1 = "cee8104904c53d39eb94fd06cbe60cb5acde7177"
+git-tree-sha1 = "d4101e848e8d3f585d61d244c2fe0c80a70e6b3b"
 uuid = "4c0bbee4-addc-4d73-81a0-b6caacae83c8"
-version = "1.0.0"
+version = "1.1.0"
 
 [[deps.ChunkCodecLibZstd]]
 deps = ["ChunkCodecCore", "Zstd_jll"]
@@ -418,9 +426,9 @@ uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
 version = "0.13.1"
 
 [[deps.CommonSolve]]
-git-tree-sha1 = "78ea4ddbcf9c241827e7035c3a03e2e456711470"
+git-tree-sha1 = "99ee296f88c12485402e37c2fd025f95ae097637"
 uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
-version = "0.2.6"
+version = "0.2.9"
 
 [[deps.CommonSubexpressions]]
 deps = ["MacroTools"]
@@ -429,9 +437,9 @@ uuid = "bbf7d656-a473-5ed7-a52c-81e309532950"
 version = "0.3.1"
 
 [[deps.CommonWorldInvalidations]]
-git-tree-sha1 = "ae52d1c52048455e85a387fbee9be553ec2b68d0"
+git-tree-sha1 = "f1697a56da59e8a2cefcbbfe71c13354a6f18c61"
 uuid = "f70d9fcc-98c5-4d4a-abd7-e4cdeebd8ca8"
-version = "1.0.0"
+version = "1.1.0"
 
 [[deps.Compat]]
 deps = ["TOML", "UUIDs"]
@@ -458,9 +466,9 @@ weakdeps = ["InverseFunctions"]
     CompositionsBaseInverseFunctionsExt = "InverseFunctions"
 
 [[deps.ConcreteStructs]]
-git-tree-sha1 = "f749037478283d372048690eb3b5f92a79432b34"
+git-tree-sha1 = "1988532cb3b4525bb718b99ba07e433bbf0e600b"
 uuid = "2569d6c7-a4a2-43d3-a901-331e8e4be471"
-version = "0.2.3"
+version = "0.2.5"
 
 [[deps.ConcurrentUtilities]]
 deps = ["Serialization", "Sockets"]
@@ -512,9 +520,9 @@ version = "1.8.2"
 
 [[deps.DataStructures]]
 deps = ["OrderedCollections"]
-git-tree-sha1 = "e86f4a2805f7f19bec5129bc9150c38208e5dc23"
+git-tree-sha1 = "6fb53a69613a0b2b68a0d12671717d307ab8b24e"
 uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.19.4"
+version = "0.19.5"
 
 [[deps.DataValueInterfaces]]
 git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
@@ -589,15 +597,15 @@ version = "1.1.0"
 
 [[deps.DiffRules]]
 deps = ["IrrationalConstants", "LogExpFunctions", "NaNMath", "Random", "SpecialFunctions"]
-git-tree-sha1 = "23163d55f885173722d1e4cf0f6110cdbaf7e272"
+git-tree-sha1 = "79a2aca180a85c690c58a020d47b426954b590f8"
 uuid = "b552c78f-8df3-52c6-915a-8e097449b14b"
-version = "1.15.1"
+version = "1.16.0"
 
 [[deps.DifferentiationInterface]]
 deps = ["ADTypes", "LinearAlgebra"]
-git-tree-sha1 = "7ae99144ea44715402c6c882bfef2adbeadbc4ce"
+git-tree-sha1 = "2348cf722e74ef177e165986c765b61217f7433a"
 uuid = "a0c0ee7d-e4b9-4e03-894e-1c5f64a51d63"
-version = "0.7.16"
+version = "0.7.19"
 
     [deps.DifferentiationInterface.extensions]
     DifferentiationInterfaceChainRulesCoreExt = "ChainRulesCore"
@@ -607,8 +615,9 @@ version = "0.7.16"
     DifferentiationInterfaceFiniteDiffExt = "FiniteDiff"
     DifferentiationInterfaceFiniteDifferencesExt = "FiniteDifferences"
     DifferentiationInterfaceForwardDiffExt = ["ForwardDiff", "DiffResults"]
-    DifferentiationInterfaceGPUArraysCoreExt = "GPUArraysCore"
+    DifferentiationInterfaceGPUArraysCoreExt = ["GPUArraysCore", "Adapt"]
     DifferentiationInterfaceGTPSAExt = "GTPSA"
+    DifferentiationInterfaceHyperHessiansExt = "HyperHessians"
     DifferentiationInterfaceMooncakeExt = "Mooncake"
     DifferentiationInterfacePolyesterForwardDiffExt = ["PolyesterForwardDiff", "ForwardDiff", "DiffResults"]
     DifferentiationInterfaceReverseDiffExt = ["ReverseDiff", "DiffResults"]
@@ -621,6 +630,7 @@ version = "0.7.16"
     DifferentiationInterfaceZygoteExt = ["Zygote", "ForwardDiff"]
 
     [deps.DifferentiationInterface.weakdeps]
+    Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
     ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
     DiffResults = "163ba53b-c6d8-5494-b064-1a9d43ac40c5"
     Diffractor = "9f5e2b26-1114-432f-b630-d3fe2085c51c"
@@ -632,6 +642,7 @@ version = "0.7.16"
     ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
     GPUArraysCore = "46192b85-c4d5-4398-a991-12ede77f4527"
     GTPSA = "b27dd330-f138-47c5-815b-40db9dd9b6e8"
+    HyperHessians = "06b494a0-c8e0-40cc-ad32-d99506a00a6c"
     Mooncake = "da2b9cff-9c12-43a0-ae48-6db2b0edb7d6"
     PolyesterForwardDiff = "98d1487c-24ca-40b6-b7ab-df2af84e126b"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
@@ -688,9 +699,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.7"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "24bbb6fc8fb87eb71c1f8d00184a60fc22c63903"
+git-tree-sha1 = "971d7831cc85f43bc9f51d615a3f7f21270c2f1d"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.8.19"
+version = "0.8.21"
 weakdeps = ["Adapt", "ChainRulesCore"]
 
     [deps.EnzymeCore.extensions]
@@ -744,15 +755,15 @@ version = "0.1.11"
 
 [[deps.Expat_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "27af30de8b5445644e8ffe3bcb0d72049c089cf1"
+git-tree-sha1 = "c307cd83373868391f3ac30b41530bc5d5d05d08"
 uuid = "2e619515-83b5-522b-bb60-26c02a35a201"
-version = "2.7.3+0"
+version = "2.8.1+0"
 
 [[deps.ExponentialUtilities]]
 deps = ["Adapt", "ArrayInterface", "GPUArraysCore", "GenericSchur", "LinearAlgebra", "PrecompileTools", "Printf", "SparseArrays", "libblastrampoline_jll"]
-git-tree-sha1 = "cc294ead6a85e975a8519dd4a0a6cb294eeb18d1"
+git-tree-sha1 = "a3c8b66f8813b721a45bfe3cf32acf793aed78bb"
 uuid = "d4d017d3-3776-5f7e-afef-a10c40355c18"
-version = "1.30.0"
+version = "1.31.0"
 weakdeps = ["StaticArrays"]
 
     [deps.ExponentialUtilities.extensions]
@@ -776,15 +787,15 @@ version = "0.4.5"
 
 [[deps.FFMPEG_jll]]
 deps = ["Artifacts", "Bzip2_jll", "FreeType2_jll", "FriBidi_jll", "JLLWrappers", "LAME_jll", "Libdl", "Ogg_jll", "OpenSSL_jll", "Opus_jll", "PCRE2_jll", "Zlib_jll", "libaom_jll", "libass_jll", "libfdk_aac_jll", "libva_jll", "libvorbis_jll", "x264_jll", "x265_jll"]
-git-tree-sha1 = "66381d7059b5f3f6162f28831854008040a4e905"
+git-tree-sha1 = "7a58e45171b63ed4782f2d36fdee8713a469e6e0"
 uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
-version = "8.0.1+1"
+version = "8.1.2+0"
 
 [[deps.FastBroadcast]]
 deps = ["ArrayInterface", "LinearAlgebra"]
-git-tree-sha1 = "e3e64918b1604ba8b1734c4a27febdfe5d09e235"
+git-tree-sha1 = "8b77b1a6d20b051c223c4907d92dbf9af7b23fa5"
 uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
-version = "1.3.1"
+version = "1.3.3"
 weakdeps = ["Polyester", "Static"]
 
     [deps.FastBroadcast.extensions]
@@ -798,14 +809,14 @@ version = "0.3.2"
 
 [[deps.FastGaussQuadrature]]
 deps = ["LinearAlgebra", "SpecialFunctions", "StaticArrays"]
-git-tree-sha1 = "0044e9f5e49a57e88205e8f30ab73928b05fe5b6"
+git-tree-sha1 = "4916117dd032ec5959b7633aedbbac408ca5ddeb"
 uuid = "442a2c76-b920-505d-bb47-c5924d526838"
-version = "1.1.0"
+version = "1.3.0"
 
 [[deps.FastPower]]
-git-tree-sha1 = "862831f78c7a48681a074ecc9aac09f2de563f71"
+git-tree-sha1 = "1b078ec113773ef9a5b07ea71ff90630eff88f58"
 uuid = "a4df4552-cc26-4903-aec0-212e50a0e84b"
-version = "1.3.1"
+version = "1.3.3"
 
     [deps.FastPower.extensions]
     FastPowerEnzymeExt = "Enzyme"
@@ -827,9 +838,9 @@ version = "1.3.1"
 
 [[deps.FileIO]]
 deps = ["Pkg", "Requires", "UUIDs"]
-git-tree-sha1 = "6522cfb3b8fe97bec632252263057996cbd3de20"
+git-tree-sha1 = "8e9c059d6857607253e837730dbf780b6b151acd"
 uuid = "5789e2e9-d7fb-5bc7-8068-2c6fae9b9549"
-version = "1.18.0"
+version = "1.19.0"
 weakdeps = ["HTTP"]
 
     [deps.FileIO.extensions]
@@ -859,9 +870,9 @@ version = "1.16.0"
 
 [[deps.FiniteDiff]]
 deps = ["ArrayInterface", "LinearAlgebra", "Setfield"]
-git-tree-sha1 = "73e879af0e767bd6dfade7c5b09d7b05657a8284"
+git-tree-sha1 = "2e5742cda07276e1834281ada4cc71b6bfc87586"
 uuid = "6a86dc24-6348-571c-b903-95158fe2bd41"
-version = "2.30.0"
+version = "2.31.1"
 
     [deps.FiniteDiff.extensions]
     FiniteDiffBandedMatricesExt = "BandedMatrices"
@@ -876,10 +887,10 @@ version = "2.30.0"
     StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 
 [[deps.FixedPointNumbers]]
-deps = ["Statistics"]
-git-tree-sha1 = "05882d6995ae5c12bb5f36dd2ed3f61c98cbb172"
+deps = ["Random", "Statistics"]
+git-tree-sha1 = "59af96b98217c6ef4ae0dfe065ac7c20831d1a84"
 uuid = "53c48c17-4a7d-5ca2-90c5-79b7896eea93"
-version = "0.8.5"
+version = "0.8.6"
 
 [[deps.Fontconfig_jll]]
 deps = ["Artifacts", "Bzip2_jll", "Expat_jll", "FreeType2_jll", "JLLWrappers", "Libdl", "Libuuid_jll", "Zlib_jll"]
@@ -894,9 +905,9 @@ version = "1.3.7"
 
 [[deps.ForwardDiff]]
 deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"]
-git-tree-sha1 = "cddeab6487248a39dae1a960fff0ac17b2a28888"
+git-tree-sha1 = "2c5d0b0e12088cde2cf84afb2784415b1ea3dfee"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "1.3.3"
+version = "1.4.1"
 weakdeps = ["StaticArrays"]
 
     [deps.ForwardDiff.extensions]
@@ -921,9 +932,9 @@ version = "1.1.3"
 
 [[deps.FunctionWrappersWrappers]]
 deps = ["FunctionWrappers", "PrecompileTools", "TruncatedStacktraces"]
-git-tree-sha1 = "ce6762f8f0e7542534f01523ae051e625cbf0468"
+git-tree-sha1 = "ffa0d00178e2004db4da70d96e17ee02e85e8966"
 uuid = "77dc65aa-8811-40c2-897b-53d922fa7daf"
-version = "1.5.0"
+version = "1.10.0"
 
     [deps.FunctionWrappersWrappers.extensions]
     FunctionWrappersWrappersEnzymeExt = ["Enzyme", "EnzymeCore"]
@@ -947,9 +958,9 @@ version = "3.4.1+1"
 
 [[deps.GPUArrays]]
 deps = ["Adapt", "GPUArraysCore", "KernelAbstractions", "LLVM", "LinearAlgebra", "Printf", "Random", "Reexport", "ScopedValues", "Serialization", "SparseArrays", "Statistics"]
-git-tree-sha1 = "99ebe50e3b6537de14d43f1c88963a2e13eff3b7"
+git-tree-sha1 = "4ea5e2aecfd52e595ff6c343410e32f83f5b9bbf"
 uuid = "0c68f7d7-f131-5f86-a1c3-88cf8149b2d7"
-version = "11.5.1"
+version = "11.5.8"
 weakdeps = ["JLD2"]
 
     [deps.GPUArrays.extensions]
@@ -963,9 +974,12 @@ version = "0.2.0"
 
 [[deps.GPUCompiler]]
 deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "PrecompileTools", "Preferences", "Scratch", "Serialization", "TOML", "Tracy", "UUIDs"]
-git-tree-sha1 = "fedfe5e7db7035271c3f58359007f971da1dde87"
+git-tree-sha1 = "b863a2e71f89328e2af69069490aec2448106e59"
 uuid = "61eb1bfa-7361-4325-ad38-22787b887f55"
-version = "1.9.1"
+version = "1.17.1"
+
+    [deps.GPUCompiler.weakdeps]
+    LLVMDowngrader_jll = "f52de702-fb25-5922-94ba-81dd59b07444"
 
 [[deps.GPUToolbox]]
 deps = ["LLVM"]
@@ -975,9 +989,9 @@ version = "1.1.1"
 
 [[deps.GR]]
 deps = ["Artifacts", "Base64", "DelimitedFiles", "Downloads", "GR_jll", "HTTP", "JSON", "Libdl", "LinearAlgebra", "Preferences", "Printf", "Qt6Wayland_jll", "Random", "Serialization", "Sockets", "TOML", "Tar", "Test", "p7zip_jll"]
-git-tree-sha1 = "44716a1a667cb867ee0e9ec8edc31c3e4aa5afdc"
+git-tree-sha1 = "f954322d5de03ec630d177cda203dcd92b6be399"
 uuid = "28b8d3ca-fb5f-59d9-8090-bfdbd6d07a71"
-version = "0.73.24"
+version = "0.73.26"
 
     [deps.GR.extensions]
     IJuliaExt = "IJulia"
@@ -987,9 +1001,9 @@ version = "0.73.24"
 
 [[deps.GR_jll]]
 deps = ["Artifacts", "Bzip2_jll", "Cairo_jll", "FFMPEG_jll", "Fontconfig_jll", "FreeType2_jll", "GLFW_jll", "JLLWrappers", "JpegTurbo_jll", "Libdl", "Libtiff_jll", "Pixman_jll", "Qt6Base_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "be8a1b8065959e24fdc1b51402f39f3b6f0f6653"
+git-tree-sha1 = "6fada551286ab6ea4ca1628cb2de9f166a2ec966"
 uuid = "d2c73de3-f751-5644-a686-071e5b155ba9"
-version = "0.73.24+0"
+version = "0.73.26+0"
 
 [[deps.GenericSchur]]
 deps = ["LinearAlgebra", "Printf"]
@@ -1017,15 +1031,15 @@ version = "1.5.0"
 
 [[deps.Git_LFS_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "bb8471f313ed941f299aa53d32a94ab3bee08844"
+git-tree-sha1 = "8c66e385d631bb934ff05e76d4a566c640c8df69"
 uuid = "020c3dae-16b3-5ae5-87b3-4cb189e250b2"
-version = "3.7.0+0"
+version = "3.7.1+0"
 
 [[deps.Git_jll]]
 deps = ["Artifacts", "Expat_jll", "JLLWrappers", "LibCURL_jll", "Libdl", "Libiconv_jll", "OpenSSL_jll", "PCRE2_jll", "Zlib_jll"]
-git-tree-sha1 = "dc34a3e3d96b4ed305b641e626dc14c12b7824b8"
+git-tree-sha1 = "0dd4cfb426924210c8f42742751cbde74b27bfa3"
 uuid = "f8c6e375-362e-5223-8a59-34ff63f689eb"
-version = "2.53.0+0"
+version = "2.54.0+0"
 
 [[deps.Glib_jll]]
 deps = ["Artifacts", "GettextRuntime_jll", "JLLWrappers", "Libdl", "Libffi_jll", "Libiconv_jll", "Libmount_jll", "PCRE2_jll", "Zlib_jll"]
@@ -1035,9 +1049,9 @@ version = "2.86.3+0"
 
 [[deps.Graphite2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "8a6dbda1fd736d60cc477d99f2e7a042acfa46e8"
+git-tree-sha1 = "69ffb934a5c5b7e086a0b4fee3427db2556fba6e"
 uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
-version = "1.3.15+0"
+version = "1.3.16+0"
 
 [[deps.Grisu]]
 git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
@@ -1046,9 +1060,9 @@ version = "1.0.2"
 
 [[deps.HSL_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "4b34f6e368aa509b244847e6b0c9b370791bab09"
+git-tree-sha1 = "4c7f8674babdebebb99ba9b20686711c5f48a395"
 uuid = "017b0a0e-03f4-516a-9b91-836bbd1904dd"
-version = "4.0.4+0"
+version = "4.0.6+0"
 
 [[deps.HTTP]]
 deps = ["Base64", "CodecZlib", "ConcurrentUtilities", "Dates", "ExceptionUnwrapping", "Logging", "LoggingExtras", "MbedTLS", "NetworkOptions", "OpenSSL", "PrecompileTools", "Random", "SimpleBufferStream", "Sockets", "URIs", "UUIDs"]
@@ -1069,15 +1083,15 @@ version = "0.2.0"
 
 [[deps.HiGHS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "Zlib_jll", "libblastrampoline_jll"]
-git-tree-sha1 = "50ed12dc8c37ebb8d2b759f21755259d8512f2bd"
+git-tree-sha1 = "5b9e398d6e4ed7e1f9a801a6f50ecca92dab0e10"
 uuid = "8fd58aa0-07eb-5a78-9b36-339c94fd15ea"
-version = "1.14.0+0"
+version = "1.15.0+0"
 
 [[deps.Hwloc_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "XML2_jll", "Xorg_libpciaccess_jll"]
-git-tree-sha1 = "baaaebd42ed9ee1bd9173cfd56910e55a8622ee1"
+git-tree-sha1 = "c35847ca5b4997fc8418836354a56c459bcf48d8"
 uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
-version = "2.13.0+1"
+version = "2.14.0+0"
 
 [[deps.IOCapture]]
 deps = ["Logging", "Random"]
@@ -1131,9 +1145,9 @@ version = "1.3.1"
 
 [[deps.Ipopt]]
 deps = ["Ipopt_jll", "LinearAlgebra", "OpenBLAS32_jll", "PrecompileTools"]
-git-tree-sha1 = "f1b9bf4b24fa1844f25fe570836f75cdb9f5245e"
+git-tree-sha1 = "f8443766032a81e1f2cddfd4f624a5650067f0d0"
 uuid = "b6b21f68-93f8-5de0-b562-5493be1d77c9"
-version = "1.14.1"
+version = "1.15.0"
 
     [deps.Ipopt.extensions]
     IpoptMathOptInterfaceExt = "MathOptInterface"
@@ -1159,9 +1173,9 @@ version = "1.0.0"
 
 [[deps.JLD2]]
 deps = ["ChunkCodecLibZlib", "ChunkCodecLibZstd", "FileIO", "MacroTools", "Mmap", "OrderedCollections", "PrecompileTools", "ScopedValues"]
-git-tree-sha1 = "941f87a0ae1b14d1ac2fa57245425b23a9d7a516"
+git-tree-sha1 = "9ebadf3f8f0de07031359917549bbdadc23f5dc3"
 uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-version = "0.6.4"
+version = "0.6.5"
 weakdeps = ["UnPack"]
 
     [deps.JLD2.extensions]
@@ -1175,15 +1189,15 @@ version = "0.1.11"
 
 [[deps.JLLWrappers]]
 deps = ["Artifacts", "Preferences"]
-git-tree-sha1 = "0533e564aae234aff59ab625543145446d8b6ec2"
+git-tree-sha1 = "7204148362dafe5fe6a273f855b8ccbe4df8173e"
 uuid = "692b3bcd-3c85-4b1f-b108-f13ce0eb3210"
-version = "1.7.1"
+version = "1.8.0"
 
 [[deps.JSON]]
 deps = ["Dates", "Logging", "Parsers", "PrecompileTools", "StructUtils", "UUIDs", "Unicode"]
-git-tree-sha1 = "67c6f1f085cb2671c93fe34244c9cccde30f7a26"
+git-tree-sha1 = "c89d196f5ffb64bfbf80985b699ea913b0d2c211"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "1.5.0"
+version = "1.6.1"
 
     [deps.JSON.extensions]
     JSONArrowExt = ["ArrowTypes"]
@@ -1211,9 +1225,9 @@ version = "0.2.1"
 
 [[deps.JpegTurbo_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "c0c9b76f3520863909825cbecdef58cd63de705a"
+git-tree-sha1 = "1dae3057da6f2b9c857afef03177bbdc7c4afe92"
 uuid = "aacddb02-875f-59d6-b918-886e6ef4fbf8"
-version = "3.1.5+0"
+version = "3.2.0+0"
 
 [[deps.JuliaNVTXCallbacks_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Pkg"]
@@ -1246,9 +1260,9 @@ version = "15.1.0"
 
 [[deps.KernelAbstractions]]
 deps = ["Adapt", "Atomix", "InteractiveUtils", "MacroTools", "PrecompileTools", "Requires", "StaticArrays", "UUIDs"]
-git-tree-sha1 = "f2e76d3ced51a2a9e185abc0b97494c7273f649f"
+git-tree-sha1 = "a5b87110fa95d711355af44832497745aa93fb52"
 uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
-version = "0.9.41"
+version = "0.9.42"
 weakdeps = ["EnzymeCore", "LinearAlgebra", "SparseArrays"]
 
     [deps.KernelAbstractions.extensions]
@@ -1258,9 +1272,9 @@ weakdeps = ["EnzymeCore", "LinearAlgebra", "SparseArrays"]
 
 [[deps.Krylov]]
 deps = ["LinearAlgebra", "Printf", "SparseArrays"]
-git-tree-sha1 = "c4d19f51afc7ba2afbe32031b8f2d21b11c9e26e"
+git-tree-sha1 = "fc2e5bc665dfa1be33fac60b5762d462bccfae7b"
 uuid = "ba0b0d4f-ebba-5204-a429-3ac8c609bfb7"
-version = "0.10.6"
+version = "0.10.8"
 
 [[deps.LAME_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1281,10 +1295,10 @@ uuid = "88015f11-f218-50d7-93a8-a6af411a945d"
 version = "4.1.0+0"
 
 [[deps.LLVM]]
-deps = ["CEnum", "LLVMExtra_jll", "Libdl", "Preferences", "Printf", "Unicode"]
-git-tree-sha1 = "69e4739502b7ab5176117e97e1664ed181c35036"
+deps = ["CEnum", "LLVMExtra_jll", "Libdl", "PrecompileTools", "Preferences", "Printf", "Unicode"]
+git-tree-sha1 = "f74a9668f02e33399baa5ed3a092b3f7a93f192e"
 uuid = "929cbde3-209d-540e-8aea-75f648917ca0"
-version = "9.4.6"
+version = "9.10.0"
 weakdeps = ["BFloat16s"]
 
     [deps.LLVM.extensions]
@@ -1292,9 +1306,9 @@ weakdeps = ["BFloat16s"]
 
 [[deps.LLVMExtra_jll]]
 deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl", "TOML"]
-git-tree-sha1 = "8e76807afb59ebb833e9b131ebf1a8c006510f33"
+git-tree-sha1 = "70c96f133c78c3cdc06234157144fab3744c6b38"
 uuid = "dad2f222-ce93-54a1-a47d-0025e8a3acab"
-version = "0.0.38+0"
+version = "0.0.43+1"
 
 [[deps.LLVMLoopInfo]]
 git-tree-sha1 = "2e5c102cfc41f48ae4740c7eca7743cc7e7b75ea"
@@ -1303,9 +1317,9 @@ version = "1.0.0"
 
 [[deps.LLVMOpenMP_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "eb62a3deb62fc6d8822c0c4bef73e4412419c5d8"
+git-tree-sha1 = "3ac157462e1e800777cc97d0eafd1bdb5356a470"
 uuid = "1d63c593-3942-5779-bab2-d838dc0a180e"
-version = "18.1.8+0"
+version = "21.1.8+0"
 
 [[deps.LaTeXStrings]]
 git-tree-sha1 = "dda21b8cbd6a6c40d9d02a73230f9d70fed6918c"
@@ -1407,9 +1421,9 @@ version = "2.42.0+0"
 
 [[deps.Libtiff_jll]]
 deps = ["Artifacts", "JLLWrappers", "JpegTurbo_jll", "LERC_jll", "Libdl", "XZ_jll", "Zlib_jll", "Zstd_jll"]
-git-tree-sha1 = "f04133fe05eff1667d2054c53d59f9122383fe05"
+git-tree-sha1 = "aebd334d06cee9f24cea70bd19a39749daf73881"
 uuid = "89763e89-9b03-5906-acba-b20f662cd828"
-version = "4.7.2+0"
+version = "4.7.3+0"
 
 [[deps.Libuuid_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1419,9 +1433,9 @@ version = "2.42.0+0"
 
 [[deps.LineSearch]]
 deps = ["ADTypes", "CommonSolve", "ConcreteStructs", "FastClosures", "LinearAlgebra", "MaybeInplace", "PrecompileTools", "SciMLBase", "SciMLJacobianOperators", "StaticArraysCore"]
-git-tree-sha1 = "25454bc65c3eec4656cbe201c3d9336af49178c7"
+git-tree-sha1 = "a7553de02835e996f80bf895a291ab981ba1f796"
 uuid = "87fe0de2-c867-4266-b59a-2f0a94fc965b"
-version = "0.1.8"
+version = "0.1.10"
 weakdeps = ["LineSearches"]
 
     [deps.LineSearch.extensions]
@@ -1429,9 +1443,9 @@ weakdeps = ["LineSearches"]
 
 [[deps.LineSearches]]
 deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Printf"]
-git-tree-sha1 = "a666999510c794fe1d9dfa629ef33366f11103aa"
+git-tree-sha1 = "cef1ba655e8c1f65af9d96c4fffe18bf1a3a3291"
 uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
-version = "7.6.1"
+version = "7.7.1"
 
 [[deps.LinearAlgebra]]
 deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
@@ -1440,9 +1454,9 @@ version = "1.12.0"
 
 [[deps.LinearOperators]]
 deps = ["FastClosures", "LinearAlgebra", "Printf", "SparseArrays", "TimerOutputs"]
-git-tree-sha1 = "ddd5a43cff2692c26f09952d33c9746cfc740d60"
+git-tree-sha1 = "4170853dfdb5ac1374ffb5fcf79c24ba5f0bb8e3"
 uuid = "5c8ed15e-5a4c-59e4-a42b-c7e8811fb125"
-version = "2.13.0"
+version = "2.14.1"
 
     [deps.LinearOperators.extensions]
     LinearOperatorsAMDGPUExt = "AMDGPU"
@@ -1451,20 +1465,24 @@ version = "2.13.0"
     LinearOperatorsJLArraysExt = "JLArrays"
     LinearOperatorsLDLFactorizationsExt = "LDLFactorizations"
     LinearOperatorsMetalExt = "Metal"
+    LinearOperatorsOpNormExt = ["Arpack", "TSVD", "GenericLinearAlgebra"]
 
     [deps.LinearOperators.weakdeps]
     AMDGPU = "21141c5a-9bdb-4563-92ae-f87d6854732e"
+    Arpack = "7d9fca2a-8960-54d3-9f78-7d1dccf2cb97"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
     ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    GenericLinearAlgebra = "14197337-ba66-59df-a3e3-ca00e7dcff7a"
     JLArrays = "27aeb0d3-9eb9-45fb-866b-73c2ecf80fcb"
     LDLFactorizations = "40e66cde-538c-5869-a4ad-c39174c6795b"
     Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
+    TSVD = "9449cd9e-2762-5aa3-a617-5413e99d722e"
 
 [[deps.LinearSolve]]
 deps = ["ArrayInterface", "ConcreteStructs", "DocStringExtensions", "EnumX", "GPUArraysCore", "InteractiveUtils", "Krylov", "Libdl", "LinearAlgebra", "MKL_jll", "Markdown", "OpenBLAS_jll", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLLogging", "SciMLOperators", "Setfield", "StaticArraysCore"]
-git-tree-sha1 = "42b5cb44317e89ef75dd841c9c8eba9045bf9ff0"
+git-tree-sha1 = "97a6bf19ef32518268ba15f166fc030df43b619a"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-version = "3.75.0"
+version = "3.79.0"
 
     [deps.LinearSolve.extensions]
     LinearSolveAMDGPUExt = "AMDGPU"
@@ -1489,12 +1507,12 @@ version = "3.75.0"
     LinearSolveKrylovKitExt = "KrylovKit"
     LinearSolveMetalExt = "Metal"
     LinearSolveMooncakeExt = "Mooncake"
-    LinearSolvePETScCSRExt = ["PETSc", "SparseArrays", "SparseMatricesCSR"]
-    LinearSolvePETScExt = ["PETSc", "SparseArrays"]
+    LinearSolvePETScExt = ["PETSc", "SparseArrays", "SparseMatricesCSR"]
     LinearSolvePETScMPIExt = ["PETSc", "PartitionedArrays", "SparseArrays", "SparseMatricesCSR"]
     LinearSolveParUExt = ["ParU_jll", "SparseArrays"]
     LinearSolvePardisoExt = ["Pardiso", "SparseArrays"]
     LinearSolveRecursiveFactorizationExt = "RecursiveFactorization"
+    LinearSolveSTRUMPACKExt = ["SparseArrays", "STRUMPACK_jll"]
     LinearSolveSparseArraysExt = "SparseArrays"
     LinearSolveSparspakExt = ["SparseArrays", "Sparspak"]
 
@@ -1526,11 +1544,18 @@ version = "3.75.0"
     Pardiso = "46dd5b70-b6fb-5a00-ae2d-e8fea33afaf2"
     PartitionedArrays = "5a9dfac6-5c52-46f7-8278-5e2210713be9"
     RecursiveFactorization = "f2c3362d-daeb-58d1-803e-2bc74f2840b4"
+    STRUMPACK_jll = "86fbd0b9-476f-557c-b766-62c724b42d8c"
     SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
     SparseMatricesCSR = "a0a7dd2c-ebf4-11e9-1f05-cf50bc540ca1"
     Sparspak = "e56a9233-b9d6-4f03-8d0f-1825330902ac"
     blis_jll = "6136c539-28a5-5bf0-87cc-b183200dce32"
 
+[[deps.Literate]]
+deps = ["Base64", "IOCapture", "JSON", "REPL"]
+git-tree-sha1 = "bb26d8b8ed0fa451ce3511e99c950653a2f31fe1"
+uuid = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
+version = "2.21.0"
+
 [[deps.LogExpFunctions]]
 deps = ["DocStringExtensions", "IrrationalConstants", "LinearAlgebra"]
 git-tree-sha1 = "13ca9e2586b89836fd20cccf56e57e2b9ae7f38f"
@@ -1597,9 +1622,9 @@ weakdeps = ["MadNLPGPU"]
 
 [[deps.MadNLP]]
 deps = ["LDLFactorizations", "LinearAlgebra", "Logging", "MUMPS_seq_jll", "NLPModels", "OpenBLAS32_jll", "Pkg", "PrecompileTools", "Printf", "SolverCore", "SparseArrays", "SuiteSparse"]
-git-tree-sha1 = "287850e28c5460a6a6ee3b83c3b2f90dbab27a79"
+git-tree-sha1 = "badbb19b4a53f15f44bbbd082676db999c388927"
 uuid = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
-version = "0.9.1"
+version = "0.9.2"
 
     [deps.MadNLP.extensions]
     MadNLPMOI = "MathOptInterface"
@@ -1637,9 +1662,9 @@ version = "0.1.3"
 
 [[deps.MaybeInplace]]
 deps = ["ArrayInterface", "LinearAlgebra", "MacroTools"]
-git-tree-sha1 = "54e2fdc38130c05b42be423e90da3bade29b74bd"
+git-tree-sha1 = "bf287c2abdd365d73c9ee86b262c6d8af0d6e2a1"
 uuid = "bb5d69b7-63fc-4a16-80bd-7e42200c7bdb"
-version = "0.1.4"
+version = "0.1.5"
 weakdeps = ["SparseArrays"]
 
     [deps.MaybeInplace.extensions]
@@ -1690,18 +1715,19 @@ version = "1.11.0"
 
 [[deps.Moshi]]
 deps = ["ExproniconLite", "Jieko"]
-git-tree-sha1 = "53f817d3e84537d84545e0ad749e483412dd6b2a"
+git-tree-sha1 = "999dfa2b4f8334c1c23bd307c2b0cb6f97c54591"
 uuid = "2e0e35c7-a2e4-4343-998d-7ef72827ed2d"
-version = "0.3.7"
+version = "0.3.8"
 
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
 version = "2025.11.4"
 
 [[deps.MuladdMacro]]
-git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab"
+deps = ["PrecompileTools"]
+git-tree-sha1 = "e8dcbeef032ba2f9051a44ac22b4e54e3a1a0099"
 uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
-version = "0.2.4"
+version = "0.2.6"
 
 [[deps.NLPModels]]
 deps = ["FastClosures", "LinearAlgebra", "LinearOperators", "Printf", "SparseArrays"]
@@ -1717,9 +1743,9 @@ version = "0.11.3"
 
 [[deps.NLPModelsKnitro]]
 deps = ["KNITRO", "NLPModels", "NLPModelsModifiers", "SolverCore"]
-git-tree-sha1 = "763bda401fa32460dae1caa48f6ab3c195d3bf3e"
+git-tree-sha1 = "0bc8ffe6b7caba17f0273a003ef38a73ffa519dd"
 uuid = "bec4dd0d-7755-52d5-9a02-22f0ffc7efcb"
-version = "0.10.3"
+version = "0.10.6"
 
 [[deps.NLPModelsModifiers]]
 deps = ["FastClosures", "LinearAlgebra", "LinearOperators", "NLPModels", "Printf", "SparseArrays"]
@@ -1751,19 +1777,19 @@ version = "3.2.2+0"
 
 [[deps.NaNMath]]
 deps = ["OpenLibm_jll"]
-git-tree-sha1 = "9b8215b1ee9e78a293f99797cd31375471b2bcae"
+git-tree-sha1 = "dbd2e8cd2c1c27f0b584f6661b4309609c5a685e"
 uuid = "77ba4419-2d1f-58cd-9bb1-8ffee604a2e3"
-version = "1.1.3"
+version = "1.1.4"
 
 [[deps.NetworkOptions]]
 uuid = "ca575930-c2e3-43a9-ace4-1e988b2c1908"
 version = "1.3.0"
 
 [[deps.NonlinearSolve]]
-deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolveQuasiNewton", "NonlinearSolveSpectralMethods", "PrecompileTools", "Preferences", "Reexport", "SciMLBase", "SciMLLogging", "Setfield", "SimpleNonlinearSolve", "StaticArraysCore", "SymbolicIndexingInterface"]
-git-tree-sha1 = "a2db21951cd1cd46a3ef8ba4bbfc84ddc9a5b3fb"
+deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "NonlinearSolveBase", "NonlinearSolveFirstOrder", "NonlinearSolveQuasiNewton", "NonlinearSolveSpectralMethods", "PrecompileTools", "Preferences", "Reexport", "SciMLBase", "Setfield", "SimpleNonlinearSolve", "StaticArraysCore", "SymbolicIndexingInterface"]
+git-tree-sha1 = "a6c5719bbb42985c72f4cacbfa49e86bab850d66"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-version = "4.18.0"
+version = "4.19.1"
 
     [deps.NonlinearSolve.extensions]
     NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
@@ -1794,9 +1820,9 @@ version = "4.18.0"
 
 [[deps.NonlinearSolveBase]]
 deps = ["ADTypes", "Adapt", "ArrayInterface", "CommonSolve", "Compat", "ConcreteStructs", "DifferentiationInterface", "EnzymeCore", "FastClosures", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "LogExpFunctions", "Markdown", "MaybeInplace", "PreallocationTools", "PrecompileTools", "Preferences", "Printf", "RecursiveArrayTools", "SciMLBase", "SciMLJacobianOperators", "SciMLLogging", "SciMLOperators", "SciMLStructures", "Setfield", "StaticArraysCore", "SymbolicIndexingInterface", "TimerOutputs"]
-git-tree-sha1 = "a19a5df29ef2b197499fc631fa1a59385ae15262"
+git-tree-sha1 = "cb824030a0d5769555704db95ee3cbe0f76ed53d"
 uuid = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
-version = "2.25.0"
+version = "2.30.3"
 
     [deps.NonlinearSolveBase.extensions]
     NonlinearSolveBaseBandedMatricesExt = "BandedMatrices"
@@ -1857,10 +1883,10 @@ uuid = "e7412a2a-1a6e-54c0-be00-318e2571c051"
 version = "1.3.6+0"
 
 [[deps.OpenBLAS32_jll]]
-deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "46cce8b42186882811da4ce1f4c7208b02deb716"
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl", "libblastrampoline_jll"]
+git-tree-sha1 = "8b492aefdd20fb9dc1ebc377ff7e5fa1591c9acc"
 uuid = "656ef2d0-ae68-5445-9ca0-591084a874a2"
-version = "0.3.30+0"
+version = "0.3.33+2"
 
 [[deps.OpenBLAS_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
@@ -1874,9 +1900,9 @@ version = "0.8.7+0"
 
 [[deps.OpenSSH_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "OpenSSL_jll", "Zlib_jll"]
-git-tree-sha1 = "57baa4b81a24c2910afbb6d853aa0685e4312bf7"
+git-tree-sha1 = "b862f484cf659efabffd601c5c920e981c942b63"
 uuid = "9bd350c2-7e96-507f-8002-3f2e150b4e1b"
-version = "10.3.1+0"
+version = "10.4.1+0"
 
 [[deps.OpenSSL]]
 deps = ["BitFlags", "Dates", "MozillaCACerts_jll", "NetworkOptions", "OpenSSL_jll", "Sockets"]
@@ -1902,9 +1928,9 @@ uuid = "91d4177d-7536-5919-b921-800302f37372"
 version = "1.6.1+0"
 
 [[deps.OrderedCollections]]
-git-tree-sha1 = "05868e21324cede2207c6f0f466b4bfef6d5e7ee"
+git-tree-sha1 = "94ba93778373a53bfd5a0caaf7d809c445292ff4"
 uuid = "bac558e1-5e72-5ebc-8fee-abe8a469f55d"
-version = "1.8.1"
+version = "1.8.2"
 
 [[deps.OrdinaryDiffEq]]
 deps = ["ADTypes", "Adapt", "ArrayInterface", "CommonSolve", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "ExponentialUtilities", "FastBroadcast", "FastClosures", "FillArrays", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "InteractiveUtils", "LineSearches", "LinearAlgebra", "LinearSolve", "Logging", "MacroTools", "MuladdMacro", "NonlinearSolve", "OrdinaryDiffEqAdamsBashforthMoulton", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqDefault", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqExplicitRK", "OrdinaryDiffEqExponentialRK", "OrdinaryDiffEqExtrapolation", "OrdinaryDiffEqFIRK", "OrdinaryDiffEqFeagin", "OrdinaryDiffEqFunctionMap", "OrdinaryDiffEqHighOrderRK", "OrdinaryDiffEqIMEXMultistep", "OrdinaryDiffEqLinear", "OrdinaryDiffEqLowOrderRK", "OrdinaryDiffEqLowStorageRK", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqNordsieck", "OrdinaryDiffEqPDIRK", "OrdinaryDiffEqPRK", "OrdinaryDiffEqQPRK", "OrdinaryDiffEqRKN", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqSDIRK", "OrdinaryDiffEqSSPRK", "OrdinaryDiffEqStabilizedIRK", "OrdinaryDiffEqStabilizedRK", "OrdinaryDiffEqSymplecticRK", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqVerner", "Polyester", "PreallocationTools", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleNonlinearSolve", "SparseArrays", "Static", "StaticArrayInterface", "StaticArrays", "TruncatedStacktraces"]
@@ -1926,9 +1952,9 @@ version = "1.26.0"
 
 [[deps.OrdinaryDiffEqCore]]
 deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "ConcreteStructs", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "FastPower", "FunctionWrappersWrappers", "InteractiveUtils", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "Polyester", "PrecompileTools", "Preferences", "Random", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLLogging", "SciMLOperators", "SciMLStructures", "Static", "SymbolicIndexingInterface", "TruncatedStacktraces"]
-git-tree-sha1 = "2b55a12e68d2e9eb9f06c00cde116356c2d3823d"
+git-tree-sha1 = "c0a22f012d3e98fd78afafe8e5152ace7a2276bc"
 uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
-version = "3.33.0"
+version = "3.33.1"
 
     [deps.OrdinaryDiffEqCore.extensions]
     OrdinaryDiffEqCoreMooncakeExt = "Mooncake"
@@ -2058,9 +2084,9 @@ version = "1.12.0"
 
 [[deps.OrdinaryDiffEqRosenbrock]]
 deps = ["ADTypes", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"]
-git-tree-sha1 = "dc4c25c30e63b9ac50e332b030c6fa99f04034ab"
+git-tree-sha1 = "3579d9272c286ffcceee794c5752e8cc7c3ae514"
 uuid = "43230ef6-c299-4910-a778-202eb28ce4ce"
-version = "1.31.0"
+version = "1.31.1"
 
 [[deps.OrdinaryDiffEqSDIRK]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"]
@@ -2123,15 +2149,15 @@ version = "0.12.3"
 
 [[deps.Parsers]]
 deps = ["Dates", "PrecompileTools", "UUIDs"]
-git-tree-sha1 = "7d2f8f21da5db6a806faf7b9b292296da42b2810"
+git-tree-sha1 = "32a4e09c5f29402573d673901778a0e03b0807b9"
 uuid = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
-version = "2.8.3"
+version = "2.8.6"
 
 [[deps.Pixman_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl"]
-git-tree-sha1 = "db76b1ecd5e9715f3d043cec13b2ec93ce015d53"
+git-tree-sha1 = "e4a6721aa89e62e5d4217c0b21bd714263779dda"
 uuid = "30392449-352a-5448-841d-b1acce4e97dc"
-version = "0.44.2+0"
+version = "0.46.4+0"
 
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"]
@@ -2194,9 +2220,9 @@ version = "1.4.3"
 
 [[deps.PreallocationTools]]
 deps = ["Adapt", "ArrayInterface", "PrecompileTools"]
-git-tree-sha1 = "e16b73bf892c55d16d53c9c0dbd0fb31cb7e25da"
+git-tree-sha1 = "0c319df479a33799d3b7566fbf14498e4e38a6f3"
 uuid = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
-version = "1.2.0"
+version = "1.2.1"
 weakdeps = ["ForwardDiff", "ReverseDiff", "SparseConnectivityTracer"]
 
     [deps.PreallocationTools.extensions]
@@ -2206,9 +2232,9 @@ weakdeps = ["ForwardDiff", "ReverseDiff", "SparseConnectivityTracer"]
 
 [[deps.PrecompileTools]]
 deps = ["Preferences"]
-git-tree-sha1 = "07a921781cab75691315adc645096ed5e370cb77"
+git-tree-sha1 = "edbeefc7a4889f528644251bdb5fc9ab5348bc2c"
 uuid = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
-version = "1.3.3"
+version = "1.3.4"
 
 [[deps.Preferences]]
 deps = ["TOML"]
@@ -2218,21 +2244,28 @@ version = "1.5.2"
 
 [[deps.PrettyTables]]
 deps = ["Crayons", "LaTeXStrings", "Markdown", "PrecompileTools", "Printf", "REPL", "Reexport", "StringManipulation", "Tables"]
-git-tree-sha1 = "624de6279ab7d94fc9f672f0068107eb6619732c"
+git-tree-sha1 = "ebf455bb866ee6737030e3d3816bb6a0683c4325"
 uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
-version = "3.3.2"
+version = "3.4.0"
 
     [deps.PrettyTables.extensions]
+    PrettyTablesExcelExt = "XLSX"
     PrettyTablesTypstryExt = "Typstry"
 
     [deps.PrettyTables.weakdeps]
     Typstry = "f0ed7684-a786-439e-b1e3-3b82803b501e"
+    XLSX = "fdbf4ff8-1666-58a4-91e7-1b58723a45e0"
 
 [[deps.Printf]]
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 version = "1.11.0"
 
+[[deps.Profile]]
+deps = ["StyledStrings"]
+uuid = "9abbd945-dff8-562f-b5e8-e1ebf5ef1b79"
+version = "1.11.0"
+
 [[deps.PtrArrays]]
 git-tree-sha1 = "4fbbafbc6251b883f4d2705356f3641f3652a7fe"
 uuid = "43287f4e-b6f4-7ad1-bb20-aadabca52c3d"
@@ -2240,15 +2273,15 @@ version = "1.4.0"
 
 [[deps.Qt6Base_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Fontconfig_jll", "Glib_jll", "JLLWrappers", "Libdl", "Libglvnd_jll", "OpenSSL_jll", "Vulkan_Loader_jll", "Xorg_libSM_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Xorg_libxcb_jll", "Xorg_xcb_util_cursor_jll", "Xorg_xcb_util_image_jll", "Xorg_xcb_util_keysyms_jll", "Xorg_xcb_util_renderutil_jll", "Xorg_xcb_util_wm_jll", "Zlib_jll", "libinput_jll", "xkbcommon_jll"]
-git-tree-sha1 = "d7a4bff94f42208ce3cf6bc8e4e7d1d663e7ee8b"
+git-tree-sha1 = "144895f6166994730ee7ff8113b981fc360638f1"
 uuid = "c0090381-4147-56d7-9ebc-da0b1113ec56"
-version = "6.10.2+1"
+version = "6.10.2+2"
 
 [[deps.Qt6Declarative_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Qt6Base_jll", "Qt6ShaderTools_jll", "Qt6Svg_jll"]
-git-tree-sha1 = "d5b7dd0e226774cbd87e2790e34def09245c7eab"
+git-tree-sha1 = "159d253ab126d5b29230cf53521899bea4ef4648"
 uuid = "629bc702-f1f5-5709-abd5-49b8460ea067"
-version = "6.10.2+1"
+version = "6.10.2+2"
 
 [[deps.Qt6ShaderTools_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Qt6Base_jll"]
@@ -2363,15 +2396,15 @@ version = "1.3.1"
 
 [[deps.ReverseDiff]]
 deps = ["ChainRulesCore", "DiffResults", "DiffRules", "ForwardDiff", "FunctionWrappers", "LinearAlgebra", "LogExpFunctions", "MacroTools", "NaNMath", "Random", "SpecialFunctions", "StaticArrays", "Statistics"]
-git-tree-sha1 = "f1b07322a8cdc0d46812473b37fb72f69ec07b22"
+git-tree-sha1 = "f5f3c09103951beae390666a18e084655996370b"
 uuid = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
-version = "1.16.2"
+version = "1.17.0"
 
 [[deps.RuntimeGeneratedFunctions]]
 deps = ["ExprTools", "SHA", "Serialization"]
-git-tree-sha1 = "cfcdc949c4660544ab0fdeed169561cb22f835f4"
+git-tree-sha1 = "bd9d6c153d0c8a120b504bfb2f3be42308cc857a"
 uuid = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47"
-version = "0.5.18"
+version = "0.5.21"
 
 [[deps.SHA]]
 uuid = "ea8e919c-243c-51af-8825-aaa63cd721ce"
@@ -2390,9 +2423,9 @@ version = "2025.9.18+0"
 
 [[deps.SciMLBase]]
 deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "CommonSolve", "ConstructionBase", "Distributed", "DocStringExtensions", "EnumX", "FunctionWrappersWrappers", "IteratorInterfaceExtensions", "LinearAlgebra", "Logging", "Markdown", "Moshi", "PreallocationTools", "PrecompileTools", "Preferences", "Printf", "RecipesBase", "RecursiveArrayTools", "Reexport", "RuntimeGeneratedFunctions", "SciMLLogging", "SciMLOperators", "SciMLPublic", "SciMLStructures", "StaticArraysCore", "Statistics", "SymbolicIndexingInterface"]
-git-tree-sha1 = "4ab0641a10918b1a2d8f59c22ff641fc478332e3"
+git-tree-sha1 = "a017ed325ac5e11438c888864fe83b124bb171b7"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "2.155.0"
+version = "2.155.1"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -2435,47 +2468,52 @@ version = "2.155.0"
 
 [[deps.SciMLJacobianOperators]]
 deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DifferentiationInterface", "FastClosures", "LinearAlgebra", "SciMLBase", "SciMLOperators"]
-git-tree-sha1 = "7156a5b51cba1bea33a82a036939ead4131f92bc"
+git-tree-sha1 = "fcdbc0214fd43dcb9201bddda77bc4cb2ddf5da7"
 uuid = "19f34311-ddf3-4b8b-af20-060888a46c0e"
-version = "0.1.13"
+version = "0.1.14"
 
 [[deps.SciMLLogging]]
 deps = ["Logging", "LoggingExtras", "Preferences"]
-git-tree-sha1 = "0161be062570af4042cf6f69e3d5d0b0555b6927"
+git-tree-sha1 = "4e1e21f14a284f892eb62923a356c70a2a0c68e1"
 uuid = "a6db7da4-7206-11f0-1eab-35f2a5dbe1d1"
-version = "1.9.1"
+version = "1.10.1"
 weakdeps = ["Tracy"]
 
     [deps.SciMLLogging.extensions]
     SciMLLoggingTracyExt = "Tracy"
 
 [[deps.SciMLOperators]]
-deps = ["Accessors", "ArrayInterface", "DocStringExtensions", "LinearAlgebra"]
-git-tree-sha1 = "234869cf9fee9258a95464b7a7065cc7be84db00"
+deps = ["Accessors", "Adapt", "ArrayInterface", "DocStringExtensions", "LinearAlgebra"]
+git-tree-sha1 = "da2b6aeedd65761a5a8bb6667acfb9c22069bee6"
 uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
-version = "1.16.0"
-weakdeps = ["SparseArrays", "StaticArraysCore"]
+version = "1.23.0"
 
     [deps.SciMLOperators.extensions]
+    SciMLOperatorsLoopVectorizationExt = "LoopVectorization"
     SciMLOperatorsSparseArraysExt = "SparseArrays"
     SciMLOperatorsStaticArraysCoreExt = "StaticArraysCore"
 
+    [deps.SciMLOperators.weakdeps]
+    LoopVectorization = "bdcacae8-1622-11e9-2a5c-532679323890"
+    SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+    StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
+
 [[deps.SciMLPublic]]
-git-tree-sha1 = "0ba076dbdce87ba230fff48ca9bca62e1f345c9b"
+git-tree-sha1 = "2b1b64add566435a768abdb3b053cac17d19ff3c"
 uuid = "431bcebd-1456-4ced-9d72-93c2757fff0b"
-version = "1.0.1"
+version = "1.2.1"
 
 [[deps.SciMLStructures]]
 deps = ["ArrayInterface", "PrecompileTools"]
-git-tree-sha1 = "607f6867d0b0553e98fc7f725c9f9f13b4d01a32"
+git-tree-sha1 = "1419128e9816e2876f7c4e93a519683b6d1d211e"
 uuid = "53ae85a6-f571-4167-b2af-e1d143709226"
-version = "1.10.0"
+version = "1.10.1"
 
 [[deps.ScopedValues]]
 deps = ["HashArrayMappedTries", "Logging"]
-git-tree-sha1 = "ac4b837d89a58c848e85e698e2a2514e9d59d8f6"
+git-tree-sha1 = "67a144433c4ce877ee6d1ada69a124d6b1ecf7be"
 uuid = "7e506255-f358-4e82-b7e4-beb19740aa63"
-version = "1.6.0"
+version = "1.6.2"
 
 [[deps.Scratch]]
 deps = ["Dates"]
@@ -2485,9 +2523,9 @@ version = "1.3.0"
 
 [[deps.SentinelArrays]]
 deps = ["Dates", "Random"]
-git-tree-sha1 = "ebe7e59b37c400f694f52b58c93d26201387da70"
+git-tree-sha1 = "084c47c7c5ce5cfecefa0a98dff69eb3646b5a80"
 uuid = "91c51154-3ec4-41a3-a24f-3f23e20d615c"
-version = "1.4.9"
+version = "1.4.10"
 
 [[deps.Serialization]]
 uuid = "9e88b42a-f829-5b0c-bbe9-9e923198166b"
@@ -2512,9 +2550,9 @@ version = "1.2.0"
 
 [[deps.SimpleNonlinearSolve]]
 deps = ["ADTypes", "ArrayInterface", "BracketingNonlinearSolve", "CommonSolve", "ConcreteStructs", "DifferentiationInterface", "FastClosures", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "MaybeInplace", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase", "Setfield", "StaticArraysCore"]
-git-tree-sha1 = "d688de789b7e643326caf9a1051376dadbcd8873"
+git-tree-sha1 = "08d8bafd57b7ea16dc1f69a8bd57db23a2ceb686"
 uuid = "727e6d20-b764-4bd8-a329-72de5adea6c7"
-version = "2.11.1"
+version = "2.12.0"
 
     [deps.SimpleNonlinearSolve.extensions]
     SimpleNonlinearSolveChainRulesCoreExt = "ChainRulesCore"
@@ -2542,9 +2580,9 @@ weakdeps = ["NLPModels"]
 
 [[deps.SortingAlgorithms]]
 deps = ["DataStructures"]
-git-tree-sha1 = "64d974c2e6fdf07f8155b5b2ca2ffa9069b608d9"
+git-tree-sha1 = "13cd91cc9be159e3f4d95b857fa2aa383b53772a"
 uuid = "a2af1166-a08f-5f64-846c-94a0d3cef48c"
-version = "1.2.2"
+version = "1.2.3"
 
 [[deps.SparseArrays]]
 deps = ["Libdl", "LinearAlgebra", "Random", "Serialization", "SuiteSparse_jll"]
@@ -2553,9 +2591,9 @@ version = "1.12.0"
 
 [[deps.SparseConnectivityTracer]]
 deps = ["ADTypes", "DocStringExtensions", "FillArrays", "LinearAlgebra", "Random", "SparseArrays"]
-git-tree-sha1 = "590b72143436e443888124aaf4026a636049e3f5"
+git-tree-sha1 = "ad4d1275eeb223cbd4d362563954a661fe12d2f7"
 uuid = "9f842d2f-2579-4b1d-911e-f412cf18a3f5"
-version = "1.2.1"
+version = "1.2.2"
 
     [deps.SparseConnectivityTracer.extensions]
     SparseConnectivityTracerChainRulesCoreExt = "ChainRulesCore"
@@ -2592,9 +2630,9 @@ version = "0.4.26"
 
 [[deps.SpecialFunctions]]
 deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "2700b235561b0335d5bef7097a111dc513b8655e"
+git-tree-sha1 = "6547cbdd8ce32efba0d21c5a40fa96d1a3548f9f"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.7.2"
+version = "2.8.0"
 weakdeps = ["ChainRulesCore"]
 
     [deps.SpecialFunctions.extensions]
@@ -2608,15 +2646,15 @@ version = "1.0.4"
 
 [[deps.Static]]
 deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools", "SciMLPublic"]
-git-tree-sha1 = "49440414711eddc7227724ae6e570c7d5559a086"
+git-tree-sha1 = "b151f033556272891e184d7d36c62518b56bbaac"
 uuid = "aedffcd0-7271-4cad-89d0-dc628f76c6d3"
-version = "1.3.1"
+version = "1.4.2"
 
 [[deps.StaticArrayInterface]]
 deps = ["ArrayInterface", "Compat", "IfElse", "LinearAlgebra", "PrecompileTools", "SciMLPublic", "Static"]
-git-tree-sha1 = "aa1ea41b3d45ac449d10477f65e2b40e3197a0d2"
+git-tree-sha1 = "2a635e15d5035c53b345077c947f31ff91744078"
 uuid = "0d7ed370-da01-4f52-bd93-41d350b8b718"
-version = "1.9.0"
+version = "1.10.0"
 
     [deps.StaticArrayInterface.extensions]
     StaticArrayInterfaceOffsetArraysExt = "OffsetArrays"
@@ -2660,15 +2698,15 @@ version = "1.8.0"
 
 [[deps.StatsBase]]
 deps = ["AliasTables", "DataAPI", "DataStructures", "IrrationalConstants", "LinearAlgebra", "LogExpFunctions", "Missings", "Printf", "Random", "SortingAlgorithms", "SparseArrays", "Statistics", "StatsAPI"]
-git-tree-sha1 = "aceda6f4e598d331548e04cc6b2124a6148138e3"
+git-tree-sha1 = "e4d7a1a0edc20af42689ea6f4f3587a2175d50ee"
 uuid = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91"
-version = "0.34.10"
+version = "0.34.12"
 
 [[deps.StrideArraysCore]]
 deps = ["ArrayInterface", "CloseOpenIntervals", "IfElse", "LayoutPointers", "LinearAlgebra", "ManualMemory", "SIMDTypes", "Static", "StaticArrayInterface", "ThreadingUtilities"]
-git-tree-sha1 = "83151ba8065a73f53ca2ae98bc7274d817aa30f2"
+git-tree-sha1 = "5316097111523c9a970596a5b33cfea5f92e8581"
 uuid = "7792a7ef-975c-4747-a70f-980b88e8d1da"
-version = "0.5.8"
+version = "0.5.9"
 
 [[deps.StringManipulation]]
 deps = ["PrecompileTools"]
@@ -2684,9 +2722,9 @@ version = "1.11.0"
 
 [[deps.StructUtils]]
 deps = ["Dates", "UUIDs"]
-git-tree-sha1 = "aab80fbf866600f3299dd7f6656d80e7be177cfe"
+git-tree-sha1 = "82bee338d650aa515f31866c460cb7e3bcef90b8"
 uuid = "ec057cc2-7a8d-4b58-b3b3-92acb9f63b42"
-version = "2.7.2"
+version = "2.8.2"
 
     [deps.StructUtils.extensions]
     StructUtilsMeasurementsExt = ["Measurements"]
@@ -2713,9 +2751,9 @@ version = "7.8.3+2"
 
 [[deps.SymbolicIndexingInterface]]
 deps = ["Accessors", "ArrayInterface", "RuntimeGeneratedFunctions", "StaticArraysCore"]
-git-tree-sha1 = "94c58884e013efff548002e8dc2fdd1cb74dfce5"
+git-tree-sha1 = "d4751bc16b120dc617719f7901a3b4e69c85b7bf"
 uuid = "2efcf032-c050-4f8e-a9bb-153293bab1f5"
-version = "0.3.46"
+version = "0.3.49"
 weakdeps = ["PrettyTables"]
 
     [deps.SymbolicIndexingInterface.extensions]
@@ -2734,9 +2772,9 @@ version = "1.0.1"
 
 [[deps.Tables]]
 deps = ["DataAPI", "DataValueInterfaces", "IteratorInterfaceExtensions", "OrderedCollections", "TableTraits"]
-git-tree-sha1 = "f2c1efbc8f3a609aadf318094f8fc5204bdaf344"
+git-tree-sha1 = "0f38a06c83f0007bbab3cf911262841c9a0f07e0"
 uuid = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
-version = "1.12.1"
+version = "1.13.0"
 
 [[deps.Tar]]
 deps = ["ArgTools", "SHA"]
@@ -2756,9 +2794,9 @@ version = "1.11.0"
 
 [[deps.ThreadingUtilities]]
 deps = ["ManualMemory"]
-git-tree-sha1 = "d969183d3d244b6c33796b5ed01ab97328f2db85"
+git-tree-sha1 = "7c73336785b21f723f5b143f6e99cf6c43b37dc1"
 uuid = "8290d209-cae3-49c0-8002-c8c24d57dab5"
-version = "0.5.5"
+version = "0.5.6"
 
 [[deps.TimerOutputs]]
 deps = ["ExprTools", "Printf"]
@@ -2952,9 +2990,9 @@ version = "0.9.12+0"
 
 [[deps.Xorg_libpciaccess_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "4909eb8f1cbf6bd4b1c30dd18b2ead9019ef2fad"
+git-tree-sha1 = "58972370b81423fc546c56a60ed1a009450177c3"
 uuid = "a65dc6b1-eb27-53a1-bb3e-dea574b5389e"
-version = "0.18.1+0"
+version = "0.19.0+0"
 
 [[deps.Xorg_libxcb_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libXau_jll", "Xorg_libXdmcp_jll"]
@@ -3012,9 +3050,9 @@ version = "1.4.7+0"
 
 [[deps.Xorg_xkeyboard_config_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_xkbcomp_jll"]
-git-tree-sha1 = "00af7ebdc563c9217ecc67776d1bbf037dbcebf4"
+git-tree-sha1 = "ed349d26affcacafbc7fc2941ace1fb98f71e715"
 uuid = "33bec58e-1273-512f-9401-5d533626f822"
-version = "2.44.0+0"
+version = "2.47.0+1"
 
 [[deps.Xorg_xtrans_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -3129,9 +3167,9 @@ version = "1.64.0+1"
 
 [[deps.oneTBB_jll]]
 deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"]
-git-tree-sha1 = "1350188a69a6e46f799d3945beef36435ed7262f"
+git-tree-sha1 = "da8c1f6eee04831f14edcfa5dae611d309807e57"
 uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e"
-version = "2022.0.0+1"
+version = "2022.3.0+0"
 
 [[deps.p7zip_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml
index 0310baf6..e534f88e 100644
--- a/docs/src/assets/Project.toml
+++ b/docs/src/assets/Project.toml
@@ -1,5 +1,6 @@
 [deps]
 ADNLPModels = "54578032-b7ea-4c30-94aa-7cbd1cce6c9a"
+BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
 CTDirect = "790bbbee-bee9-49ee-8912-a9de031322d5"
 CTFlows = "1c39547c-7794-42f7-af83-d98194f657c2"
@@ -15,6 +16,7 @@ DocumenterMermaid = "a078cd44-4d9c-4618-b545-3ab9d77f9177"
 ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2"
 JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
 JSON3 = "0f8b85d8-7281-11e9-16c2-39a750bddbf1"
+Literate = "98b081ad-f1c9-55d3-8b20-4c87d4299306"
 MadNCL = "434a0bcb-5a7c-42b2-a9d3-9e3f760e7af0"
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
 MadNLPGPU = "d72a61cc-809d-412f-99be-fd81f4b8a598"
@@ -28,10 +30,11 @@ UnoSolver = "1baa60ac-02f7-4b39-a7a8-2f4f58486b05"
 
 [compat]
 ADNLPModels = "0.8"
-CTBase = "0.18"
+BenchmarkTools = "1"
+CTBase = "=0.18.8"
 CTDirect = "1"
 CTFlows = "0.8"
-CTModels = "0.10"
+CTModels = "=0.10.1"
 CTParser = "0.8"
 CTSolvers = "0.4"
 CUDA = "5"
@@ -43,6 +46,7 @@ DocumenterMermaid = "0.2"
 ExaModels = "0.9"
 JLD2 = "0.6"
 JSON3 = "1"
+Literate = "2"
 MadNCL = "0.2"
 MadNLP = "0.9"
 MadNLPGPU = "0.8"
diff --git a/docs/src/notebooks/tutorial.ipynb b/docs/src/notebooks/tutorial.ipynb
new file mode 100644
index 00000000..652eba19
--- /dev/null
+++ b/docs/src/notebooks/tutorial.ipynb
@@ -0,0 +1,716 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# OptimalControl.jl — a guided tour\n",
+    "\n",
+    "\n",
+    "This tutorial is a guided tour of [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl), part of the [control-toolbox](https://control-toolbox.org) ecosystem. We follow two problems end to end: a simple **double integrator** for modelling, initialisation and the **indirect** (Pontryagin) method, and the **Goddard rocket** for the **direct** method in depth, grid continuation and GPU solving. Advanced topics are linked at the end.\n",
+    "\n",
+    "It is written for readers with a background in optimal control, ODEs or optimisation. By the end you will be able to define an optimal control problem, solve it by both the direct and indirect methods, and visualise the result — all in a few lines of code."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The problem, and installing the tools\n",
+    "\n",
+    "An **optimal control problem** (OCP) in Bolza form reads\n",
+    "\n",
+    "$$\n",
+    "J(x, u) = g(x(t_0), x(t_f)) + \\int_{t_0}^{t_f} f^{0}(t, x(t), u(t))\\,\\mathrm{d}t \\;\\to\\; \\min,\n",
+    "$$\n",
+    "\n",
+    "subject to the controlled dynamics $\\dot{x}(t) = f(t, x(t), u(t))$ and, possibly, box / path / boundary constraints. When $g = 0$ the cost is of **Lagrange** form; when $f^0 = 0$, of **Mayer** form.\n",
+    "\n",
+    "More generally, the times $t_0$ and $t_f$ may be free (optimisation variables), and a vector $v$ of additional parameters can enter the cost, dynamics and constraints. The full problem then reads\n",
+    "\n",
+    "$$\n",
+    "\\min_{x,u,v}\\; g(x(t_0), x(t_f), v) + \\int_{t_0}^{t_f} f^{0}(t, x(t), u(t), v)\\,\\mathrm{d}t,\n",
+    "$$\n",
+    "\n",
+    "subject to $\\dot{x}(t) = f(t, x(t), u(t), v)$, box / path / boundary constraints.\n",
+    "\n",
+    "OptimalControl.jl is the core of the [control-toolbox](https://control-toolbox.org) ecosystem, a modular suite of Julia packages — CTBase (base types & exceptions), CTParser (DSL parsing), CTModels (problem data structures), CTDirect (discretisation & NLP transcription), CTFlows (Hamiltonian flows for indirect methods), and CTSolvers (solver orchestration) — that can also be used individually.\n",
+    "\n",
+    "Installation is a single package:\n",
+    "\n",
+    "```julia\n",
+    "import Pkg\n",
+    "Pkg.add(\"OptimalControl\")\n",
+    "```\n",
+    "\n",
+    "We load OptimalControl.jl to model the problem, a solver backend ([NLPModelsIpopt.jl](https://jso.dev/NLPModelsIpopt.jl)), and [Plots.jl](https://docs.juliaplots.org)."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using OptimalControl\n",
+    "using NLPModelsIpopt\n",
+    "using Plots"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Defining a problem: `@def` vs macro-free\n",
+    "\n",
+    "Our running example: a wagon of unit mass on a frictionless rail, state $x = (q, v)$ (position, velocity), acceleration controlled by a force $u$. We start at $(-1, 0)$, must reach $(0, 0)$ at $t_f = 1$, and minimise the transfer energy\n",
+    "\n",
+    "$$\n",
+    "\\frac{1}{2}\\int_0^1 u^2(t)\\,\\mathrm{d}t,\n",
+    "$$\n",
+    "\n",
+    "subject to the dynamics\n",
+    "\n",
+    "$$\n",
+    "\\dot{q}(t) = v(t), \\qquad \\dot{v}(t) = u(t).\n",
+    "$$"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "t0 = 0; tf = 1; x0 = [-1, 0]; xf = [0, 0];"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### The `@def` macro\n",
+    "\n",
+    "The [`@def`](https://control-toolbox.org/OptimalControl.jl/stable/manual-abstract.html) macro lets us write the problem almost exactly as the mathematics:\n",
+    "\n",
+    "Each line of the `@def` block mirrors a piece of the mathematical formulation — time, state, control, dynamics, boundary conditions, then cost — in the same order one would write them on paper. Unicode symbols (`∈`, `R²`, `ẋ`, `∫`, `→`) make the code read like the maths; plain ASCII alternatives (`R^2`, `derivative`, `integral`, `=>`) are available for keyboards or workflows that prefer them."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ocp = @def begin\n",
+    "    t ∈ [t0, tf], time\n",
+    "    x = (q, v) ∈ R², state\n",
+    "    u ∈ R, control\n",
+    "\n",
+    "    x(t0) == x0\n",
+    "    x(tf) == xf\n",
+    "\n",
+    "    ẋ(t) == [v(t), u(t)]\n",
+    "\n",
+    "    0.5∫( u(t)^2 ) → min\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### The same problem with the macro-free (functional) API\n",
+    "\n",
+    "The [functional API](https://control-toolbox.org/OptimalControl.jl/stable/manual-macro-free.html) builds the *same* model step by step with plain functions — useful for programmatic problem generation or macro-free library code."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pre = OptimalControl.PreModel()\n",
+    "\n",
+    "time!(pre; t0=t0, tf=tf)\n",
+    "state!(pre, 2, \"x\", [\"q\", \"v\"])\n",
+    "control!(pre, 1)\n",
+    "\n",
+    "function f_energy!(dx, t, x, u, v)\n",
+    "    dx[1] = x[2]\n",
+    "    dx[2] = u[1]\n",
+    "    return nothing\n",
+    "end\n",
+    "dynamics!(pre, f_energy!)\n",
+    "\n",
+    "function boundary_energy!(b, x0_, xf_, v)\n",
+    "    b[1] = x0_[1] - x0[1]\n",
+    "    b[2] = x0_[2] - x0[2]\n",
+    "    b[3] = xf_[1] - xf[1]\n",
+    "    b[4] = xf_[2] - xf[2]\n",
+    "    return nothing\n",
+    "end\n",
+    "constraint!(pre, :boundary; f=boundary_energy!, lb=zeros(4), ub=zeros(4), label=:endpoint)\n",
+    "\n",
+    "lagrange_energy(t, x, u, v) = 0.5 * u[1]^2\n",
+    "objective!(pre, :min; lagrange=lagrange_energy)\n",
+    "\n",
+    "time_dependence!(pre; autonomous=true)\n",
+    "\n",
+    "ocp_func = build(pre)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### What the macro actually does\n",
+    "\n",
+    "**Key message:** `@def` *translates the expression* into the very same functional calls, and **additionally records the symbolic definition**. We can see the difference directly: the macro keeps the DSL expression, whereas the functional API stores an empty definition."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "definition(ocp)          # the macro records the full DSL expression"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "has_abstract_definition(ocp_func)   # false: functional API stores no abstract definition"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Two things to keep in mind:** (1) in the functional API, callbacks are always vector-valued — even when the control is scalar, one writes `u[1]` — not `u` — inside `f_energy!` or `lagrange_energy`; (2) the functional API currently works only with the `:adnlp` modeler, **not** `:exa` (needed for GPU solving) — one more reason to prefer `@def` when GPU execution is contemplated (more in the GPU section)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## First solve, initial guess, and the costate\n",
+    "\n",
+    "Solving is one call, plotting another."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "direct_sol = solve(ocp)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plot(direct_sol; size=(800, 600))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### The default initial guess\n",
+    "\n",
+    "With no initial guess, every variable is initialised to `0.1`. We can *see* the initial guess without optimising, by stopping the solver immediately with `max_iter=0`:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "sol_init = solve(ocp; init=nothing, max_iter=0, display=false)\n",
+    "plot(sol_init; size=(800, 600))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Notice the right-hand column: the costate is already there.** Even though we only ever provide the state, control and (optional) variable, the solver initialises the **adjoint** internally. After optimisation, this right-column costate is exactly the **adjoint $p$ of Pontryagin's Maximum Principle** — the same $p$ we will reuse to start the indirect method in the indirect section. This closes the loop between the direct and indirect methods."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Providing our own initial guess\n",
+    "\n",
+    "The recommended way to provide an initial guess is the `@init` macro, using the labels from the `@def` block (`q`, `v`, `u` here):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "ig = @init ocp begin\n",
+    "    q(t) := -1 + t\n",
+    "    v(t) := 0\n",
+    "    u(t) := 0\n",
+    "end\n",
+    "\n",
+    "sol = solve(ocp; init=ig, display=false)\n",
+    "println(\"iterations, default guess: \", iterations(direct_sol))\n",
+    "println(\"iterations, @init guess:   \", iterations(sol))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "In this case both guesses give **1 iteration**: the double integrator is a *linear-quadratic* problem, so the NLP is quadratic and Ipopt solves it in a single step regardless of the starting point. Warm-starting only pays off on genuinely nonlinear problems — we will see this with the **Goddard rocket** in the next section."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For all the ways to specify an initial guess, see [Set an initial guess](https://control-toolbox.org/OptimalControl.jl/stable/manual-initial-guess.html).\n",
+    "**Note:** there is currently no way to initialise the costate directly — only state, control and variable can be provided through `@init`. The solver initialises the adjoint internally (as we saw above). Costate initialisation is a planned feature."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Direct method in depth: Goddard\n",
+    "\n",
+    "### Discretise optimal control problems\n",
+    "\n",
+    "The **direct** method turns the infinite-dimensional OCP into a finite-dimensional nonlinear program (NLP) by discretising time (Runge–Kutta / collocation) on a grid, then hands the NLP to a solver. It is robust and easy to use.\n",
+    "\n",
+    "Concretely, time is discretised on a uniform grid $t_0 < t_1 < \\dots < t_N = t_f$ with step $h = (t_f - t_0)/N$. The (explicit) Euler scheme, for instance, replaces the dynamics by\n",
+    "\n",
+    "$$\n",
+    "x_{i} = x_{i-1} + h\\,f(t_{i-1}, x_{i-1}, u_{i-1}), \\quad i = 1, \\dots, N,\n",
+    "$$\n",
+    "\n",
+    "and the integral cost by the corresponding rectangle sum\n",
+    "\n",
+    "$$\n",
+    "h\\sum_{i=0}^{N-1} f^{0}(t_i, x_i, u_i).\n",
+    "$$\n",
+    "\n",
+    "The continuous OCP thus becomes a finite-dimensional NLP in the variables $X = (x_0, \\dots, x_N, u_0, \\dots, u_N)$, which is passed to an NLP solver such as [Ipopt](https://coin-or.github.io/Ipopt). Higher-order schemes (midpoint, Gauss–Legendre collocation) follow the same principle with different quadrature and interpolation formulas — `solve` defaults to the second-order `:midpoint` scheme, not Euler.\n",
+    "\n",
+    "### The Goddard rocket problem\n",
+    "\n",
+    "To demonstrate convergence behaviour and warm-starting, we need a genuinely nonlinear problem. The **Goddard rocket** — maximise the final altitude, with free final time and a singular arc — is a classic test case."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Goddard data and dynamics (F0: drift, F1: thrust)\n",
+    "const r0 = 1\n",
+    "const v0 = 0\n",
+    "const m0 = 1\n",
+    "const mf = 0.6\n",
+    "const Cd = 310\n",
+    "const Tmax = 3.5\n",
+    "const β = 500\n",
+    "const b = 2\n",
+    "\n",
+    "F0(x) = begin\n",
+    "    r, v, m = x\n",
+    "    D = Cd * v^2 * exp(-β * (r - 1))\n",
+    "    [v, -D/m - 1/r^2, 0]\n",
+    "end\n",
+    "F1(x) = begin\n",
+    "    r, v, m = x\n",
+    "    [0, Tmax/m, -b*Tmax]\n",
+    "end\n",
+    "\n",
+    "goddard = @def begin\n",
+    "    tf ∈ R, variable\n",
+    "    t ∈ [t0, tf], time\n",
+    "    x = (r, v, m) ∈ R³, state\n",
+    "    u ∈ R, control\n",
+    "\n",
+    "    x(t0) == [r0, v0, m0]\n",
+    "    m(tf) == mf\n",
+    "    0 ≤ u(t) ≤ 1\n",
+    "    r(t) ≥ r0\n",
+    "\n",
+    "    ẋ(t) == F0(x(t)) + u(t) * F1(x(t))\n",
+    "\n",
+    "    r(tf) → max\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Choosing a solver is trivial\n",
+    "\n",
+    "`solve` uses the defaults (collocation, ADNLP modeler, Ipopt, CPU). Switching solver is just loading a package and passing a token (see [Solve a problem](https://control-toolbox.org/OptimalControl.jl/stable/manual-solve.html)):"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using MadNLP\n",
+    "\n",
+    "sol_ipopt  = solve(goddard;          grid_size=250, display=false)\n",
+    "sol_madnlp = solve(goddard, :madnlp; grid_size=250, display=false)\n",
+    "\n",
+    "println(\"Ipopt  : r(tf) = \", objective(sol_ipopt),  \", \", iterations(sol_ipopt),  \" iters\")\n",
+    "println(\"MadNLP : r(tf) = \", objective(sol_madnlp), \", \", iterations(sol_madnlp), \" iters\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The available methods and their options can be inspected with `methods()` and `describe(:collocation)`; we will not dwell on them here."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Grid continuation by warm-starting\n",
+    "\n",
+    "A solution can be passed **directly** as the initial guess of another solve — it is interpolated onto the new grid. This makes discrete continuation trivial and ties back to the initialisation above. On this nonlinear problem it genuinely **pays**: we compare reaching a fine grid of 1000 two ways:\n",
+    "\n",
+    "1. **cold start** — solve `grid_size=1000` directly;\n",
+    "2. **cascade** — solve `grid_size=50` first, then `grid_size=1000` warm-started with that solution."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# solutions computed once, reused for iteration counts and the overlay plot\n",
+    "sol_cold = solve(goddard; grid_size=1000, display=false)\n",
+    "\n",
+    "# warm cascade: grid 50 first, then grid 1000 initialised from it\n",
+    "s50   = solve(goddard; grid_size=50, display=false)\n",
+    "s1000 = solve(goddard; grid_size=1000, init=s50, display=false)\n",
+    "\n",
+    "println(\"cold    grid 1000        : \", iterations(sol_cold), \" iters\")\n",
+    "println(\"cascade grid 50 (warm-up): \", iterations(s50),      \" iters\")\n",
+    "println(\"cascade grid 1000 (warm) : \", iterations(s1000),    \" iters\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**Message:** what matters is the iteration count *at the expensive grid* — the warm-started `iterations(s1000)` is well below the cold `iterations(sol_cold)`, even though the cheap `grid_size=50` warm-up adds iterations of its own to the running total; since a grid-50 iteration is far cheaper than a grid-1000 iteration, the cascade still wins on wall-clock time. Overlay the successive solutions to watch convergence:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "plt = plot(s50;  label=\"50\", size=(800, 800))\n",
+    "plot!(plt, s1000; label=\"1000\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This is grid-refinement warm-starting. The very same mechanism drives **parametric** continuation (homotopy on a physical parameter, e.g. maximum thrust): <https://control-toolbox.org/Tutorials.jl/stable/tutorial-continuation.html>."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Comparison with a bang-bang strategy"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "How much better is the optimal solution compared to a naive strategy? We simulate **full thrust until fuel depletion, then coast to apogee** — a bang-bang profile with no optimisation, just two ODE integrations with callbacks."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using OrdinaryDiffEq   # ODE solver (callbacks for the bang-bang simulation)\n",
+    "\n",
+    "# Phase 1: u = 1, stop when m = mf (fuel depleted)\n",
+    "bang1!(dx, x, p, t) = (dx[:] = F0(x) + F1(x))\n",
+    "cb_fuel = ContinuousCallback((u, t, int) -> u[3] - mf, terminate!)\n",
+    "sol_bang1 = solve(ODEProblem(bang1!, [r0, v0, m0], (t0, 100.0)), Tsit5(); callback=cb_fuel, reltol=1e-8, abstol=1e-8)\n",
+    "t1_bang, x1_bang = sol_bang1.t[end], sol_bang1[:, end]\n",
+    "\n",
+    "# Phase 2: u = 0, stop when v = 0 (apogee)\n",
+    "bang2!(dx, x, p, t) = (dx[:] = F0(x))\n",
+    "cb_apogee = ContinuousCallback((u, t, int) -> u[2], terminate!)\n",
+    "sol_bang2 = solve(ODEProblem(bang2!, x1_bang, (t1_bang, 1000.0)), Tsit5(); callback=cb_apogee, reltol=1e-8, abstol=1e-8)\n",
+    "tf_bang, rf_bang = sol_bang2.t[end], sol_bang2[1, end]\n",
+    "\n",
+    "println(\"Bang-bang: r(tf) = \", round(rf_bang, digits=6), \"  (t1=\", round(t1_bang, digits=4), \", tf=\", round(tf_bang, digits=4), \")\")\n",
+    "println(\"Optimal:   r(tf) = \", round(objective(sol_cold), digits=6), \"  (           tf=\", round(variable(sol_cold), digits=4), \")\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The optimal thrust profile uses a **singular arc** — it does not simply push at the maximum. Overlaying the two trajectories on the altitude–velocity plane makes the difference visible:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# assemble the bang-bang trajectory as (t, r, v, m) for plotting\n",
+    "t_bang = [sol_bang1.t; sol_bang2.t]\n",
+    "r_bang = [sol_bang1[1, :]; sol_bang2[1, :]]\n",
+    "\n",
+    "plt_bang = plot(sol_cold; label=\"optimal\", linewidth=2, color=1)\n",
+    "plot!(plt_bang[1], t_bang, r_bang; label=\"bang-bang\", linestyle=:dash, linewidth=2, color=2)\n",
+    "plot(plt_bang[1]; legend=:bottomright, xlabel=\"time\", ylabel=\"altitude\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Solving on a GPU\n",
+    "\n",
+    "Moving to the GPU is a single token, `:gpu`, which auto-completes to `(:collocation, :exa, :madnlp, :gpu)`. It requires the `:exa` modeler (hence `@def`, not the macro-free API — cf. the definition section) plus a CUDA-capable GPU.\n",
+    "\n",
+    "In a seminar or on Binder there is usually **no functional GPU**, so the call is *expected to fail* — that is the pedagogical point: the `:gpu` token needs a specific setup. We wrap it in a `try/catch` so the tour keeps running and shows the raised exception."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using MadNLPGPU\n",
+    "using CUDA\n",
+    "\n",
+    "try\n",
+    "    global sol_gpu = solve(goddard, :gpu; grid_size=1000, display=false)\n",
+    "    println(\"GPU solve succeeded — a functional GPU is available.\")\n",
+    "catch e\n",
+    "    println(\"GPU solve failed, as expected without a functional GPU.\")\n",
+    "    println(\"CUDA.functional() = \", CUDA.functional())\n",
+    "    println(\"Exception: \", first(sprint(showerror, e), 400))\n",
+    "end"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "For the full GPU setup, see [Solve on GPU](https://control-toolbox.org/OptimalControl.jl/stable/manual-solve-gpu.html)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The indirect method\n",
+    "\n",
+    "We now return to the **double integrator** `ocp` from the earlier sections. Its shooting has just two unknowns and is initialised by the direct costate above, which makes it ideal to *see* the indirect method. (The Goddard shooting is a *structured multi-arc* problem — see the links in the last section.)\n",
+    "\n",
+    "In control-toolbox we systematically pair the direct method with the **indirect** one, based on Pontryagin's Maximum Principle (PMP), with pseudo-Hamiltonian\n",
+    "\n",
+    "$$\n",
+    "H(x,p,u) = p\\,f(x,u) + p^0 f^0(x,u) \\qquad (\\text{normal case } p^0 = -1).\n",
+    "$$\n",
+    "\n",
+    "The PMP gives the maximising control in feedback form\n",
+    "\n",
+    "$$\n",
+    "u(x,p) = \\arg\\max_u H,\n",
+    "$$\n",
+    "\n",
+    "and the optimal trajectory solves a boundary value problem that we recast as a **shooting equation**\n",
+    "\n",
+    "$$\n",
+    "S(p_0) = 0.\n",
+    "$$\n",
+    "\n",
+    "The indirect method proceeds in three steps:\n",
+    "\n",
+    "1. **Maximising control.** The PMP yields the control in feedback form $u(x, p) = \\arg\\max_u H(x, p, u)$. Substituting back gives the maximised Hamiltonian\n",
+    "\n",
+    "   $$\n",
+    "   \\mathbf{H}(x, p) = H(x, p, u(x, p)).\n",
+    "   $$\n",
+    "\n",
+    "2. **Boundary value problem.** The optimal trajectory satisfies the Hamiltonian system\n",
+    "\n",
+    "   $$\n",
+    "   \\dot{x} = \\nabla_p \\mathbf{H}, \\qquad \\dot{p} = -\\nabla_x \\mathbf{H},\n",
+    "   $$\n",
+    "\n",
+    "   with boundary conditions $x(t_0) = x_0$, $x(t_f) = x_f$.\n",
+    "\n",
+    "3. **Shooting function.** Let $\\varphi_{t_0, x_0, p_0}(\\cdot)$ denote the flow of the Hamiltonian vector field from $(x_0, p_0)$. The shooting function\n",
+    "\n",
+    "   $$\n",
+    "   S(p_0) = \\pi(\\varphi_{t_0, x_0, p_0}(t_f)) - x_f, \\qquad \\pi(x, p) = x,\n",
+    "   $$\n",
+    "\n",
+    "   measures the miss at $t_f$: solving the BVP reduces to finding $p_0$ such that $S(p_0) = 0$.\n",
+    "\n",
+    "For the energy problem, $H = p_1 v + p_2 u - u^2/2$, so the maximiser is $u = p_2$."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "using OrdinaryDiffEq   # ODE solver (Hamiltonian flow)\n",
+    "using NonlinearSolve   # nonlinear equations (shooting)\n",
+    "\n",
+    "# maximising control in feedback form\n",
+    "u_max(x, p) = p[2]\n",
+    "\n",
+    "# Hamiltonian flow of the OCP\n",
+    "φ = Flow(ocp, u_max);\n",
+    "\n",
+    "# state projection π(x, p) = x\n",
+    "proj((x, p)) = x\n",
+    "\n",
+    "# shooting function\n",
+    "S(p0) = proj(φ(t0, x0, p0, tf)) - xf"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "**The shooting is initialised with the costate of the direct solution** — the very adjoint we highlighted above:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "nle!(s, p0, _) = (s[:] = S(p0))\n",
+    "\n",
+    "p_of_t   = costate(direct_sol)     # costate as a function of time\n",
+    "p0_guess = p_of_t(t0)              # initial costate from the direct method\n",
+    "\n",
+    "prob = NonlinearProblem(nle!, p0_guess)\n",
+    "shooting_sol = solve(prob; show_trace=Val(true))\n",
+    "p0_sol = shooting_sol.u\n",
+    "\n",
+    "println(\"costate p0 = \", p0_sol)\n",
+    "println(\"shoot S(p0) = \", S(p0_sol))"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Reconstruct the indirect solution from the flow and overlay it with the direct solution:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "indirect_sol = φ((t0, tf), x0, p0_sol; saveat=range(t0, tf, 100))\n",
+    "\n",
+    "plt_compare = plot(direct_sol; label=\"direct\", size=(800, 600))\n",
+    "plot!(plt_compare, indirect_sol; label=\"indirect\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "See [Compute flows from optimal control problems](https://control-toolbox.org/OptimalControl.jl/stable/manual-flow-ocp.html) for the flow construction, and the [indirect simple shooting tutorial](https://control-toolbox.org/Tutorials.jl/stable/)."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Going further\n",
+    "\n",
+    "**Variables & parameters.** Beyond the control, one can optimise **parameters** naturally, both in an OCP (the `variable` keyword of the DSL) and in a differential-constraint optimisation problem **without any control** (a *control-free* problem).\n",
+    "See [control-free problems](https://control-toolbox.org/OptimalControl.jl/stable/example-control-free.html).\n",
+    "\n",
+    "**Advanced examples** (each does both direct and indirect):\n",
+    "\n",
+    "- Singular control (control-affine systems) — <https://control-toolbox.org/OptimalControl.jl/stable/example-singular-control.html>\n",
+    "- State constraint — <https://control-toolbox.org/OptimalControl.jl/stable/example-state-constraint.html>\n",
+    "- Goddard problem — free final time, a singular arc, a state constraint and a structured shooting all at once — <https://control-toolbox.org/Tutorials.jl/stable/tutorial-goddard.html>\n",
+    "\n",
+    "**Discrete continuation** — warm-starting across a family of problems (homotopy on a physical parameter), the grown-up version of the grid continuation above: <https://control-toolbox.org/Tutorials.jl/stable/tutorial-continuation.html>."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "---\n",
+    "\n",
+    "*This notebook was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Julia 1.12.6",
+   "language": "julia",
+   "name": "julia-1.12"
+  },
+  "language_info": {
+   "file_extension": ".jl",
+   "mimetype": "application/julia",
+   "name": "julia",
+   "version": "1.12.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 3
+}
\ No newline at end of file
diff --git a/docs/src/scripts/tutorial.jl b/docs/src/scripts/tutorial.jl
new file mode 100644
index 00000000..15878d3d
--- /dev/null
+++ b/docs/src/scripts/tutorial.jl
@@ -0,0 +1,197 @@
+using OptimalControl
+using NLPModelsIpopt
+using Plots
+
+t0 = 0; tf = 1; x0 = [-1, 0]; xf = [0, 0];
+
+ocp = @def begin
+    t ∈ [t0, tf], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+
+    x(t0) == x0
+    x(tf) == xf
+
+    ẋ(t) == [v(t), u(t)]
+
+    0.5∫( u(t)^2 ) → min
+end
+
+pre = OptimalControl.PreModel()
+
+time!(pre; t0=t0, tf=tf)
+state!(pre, 2, "x", ["q", "v"])
+control!(pre, 1)
+
+function f_energy!(dx, t, x, u, v)
+    dx[1] = x[2]
+    dx[2] = u[1]
+    return nothing
+end
+dynamics!(pre, f_energy!)
+
+function boundary_energy!(b, x0_, xf_, v)
+    b[1] = x0_[1] - x0[1]
+    b[2] = x0_[2] - x0[2]
+    b[3] = xf_[1] - xf[1]
+    b[4] = xf_[2] - xf[2]
+    return nothing
+end
+constraint!(pre, :boundary; f=boundary_energy!, lb=zeros(4), ub=zeros(4), label=:endpoint)
+
+lagrange_energy(t, x, u, v) = 0.5 * u[1]^2
+objective!(pre, :min; lagrange=lagrange_energy)
+
+time_dependence!(pre; autonomous=true)
+
+ocp_func = build(pre)
+
+definition(ocp)          # the macro records the full DSL expression
+
+has_abstract_definition(ocp_func)   # false: functional API stores no abstract definition
+
+direct_sol = solve(ocp)
+
+plot(direct_sol; size=(800, 600))
+
+sol_init = solve(ocp; init=nothing, max_iter=0, display=false)
+plot(sol_init; size=(800, 600))
+
+ig = @init ocp begin
+    q(t) := -1 + t
+    v(t) := 0
+    u(t) := 0
+end
+
+sol = solve(ocp; init=ig, display=false)
+println("iterations, default guess: ", iterations(direct_sol))
+println("iterations, @init guess:   ", iterations(sol))
+
+# Goddard data and dynamics (F0: drift, F1: thrust)
+const r0 = 1
+const v0 = 0
+const m0 = 1
+const mf = 0.6
+const Cd = 310
+const Tmax = 3.5
+const β = 500
+const b = 2
+
+F0(x) = begin
+    r, v, m = x
+    D = Cd * v^2 * exp(-β * (r - 1))
+    [v, -D/m - 1/r^2, 0]
+end
+F1(x) = begin
+    r, v, m = x
+    [0, Tmax/m, -b*Tmax]
+end
+
+goddard = @def begin
+    tf ∈ R, variable
+    t ∈ [t0, tf], time
+    x = (r, v, m) ∈ R³, state
+    u ∈ R, control
+
+    x(t0) == [r0, v0, m0]
+    m(tf) == mf
+    0 ≤ u(t) ≤ 1
+    r(t) ≥ r0
+
+    ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
+
+    r(tf) → max
+end
+
+using MadNLP
+
+sol_ipopt  = solve(goddard;          grid_size=250, display=false)
+sol_madnlp = solve(goddard, :madnlp; grid_size=250, display=false)
+
+println("Ipopt  : r(tf) = ", objective(sol_ipopt),  ", ", iterations(sol_ipopt),  " iters")
+println("MadNLP : r(tf) = ", objective(sol_madnlp), ", ", iterations(sol_madnlp), " iters")
+
+# solutions computed once, reused for iteration counts and the overlay plot
+sol_cold = solve(goddard; grid_size=1000, display=false)
+
+# warm cascade: grid 50 first, then grid 1000 initialised from it
+s50   = solve(goddard; grid_size=50, display=false)
+s1000 = solve(goddard; grid_size=1000, init=s50, display=false)
+
+println("cold    grid 1000        : ", iterations(sol_cold), " iters")
+println("cascade grid 50 (warm-up): ", iterations(s50),      " iters")
+println("cascade grid 1000 (warm) : ", iterations(s1000),    " iters")
+
+plt = plot(s50;  label="50", size=(800, 800))
+plot!(plt, s1000; label="1000")
+
+using OrdinaryDiffEq   # ODE solver (callbacks for the bang-bang simulation)
+
+# Phase 1: u = 1, stop when m = mf (fuel depleted)
+bang1!(dx, x, p, t) = (dx[:] = F0(x) + F1(x))
+cb_fuel = ContinuousCallback((u, t, int) -> u[3] - mf, terminate!)
+sol_bang1 = solve(ODEProblem(bang1!, [r0, v0, m0], (t0, 100.0)), Tsit5(); callback=cb_fuel, reltol=1e-8, abstol=1e-8)
+t1_bang, x1_bang = sol_bang1.t[end], sol_bang1[:, end]
+
+# Phase 2: u = 0, stop when v = 0 (apogee)
+bang2!(dx, x, p, t) = (dx[:] = F0(x))
+cb_apogee = ContinuousCallback((u, t, int) -> u[2], terminate!)
+sol_bang2 = solve(ODEProblem(bang2!, x1_bang, (t1_bang, 1000.0)), Tsit5(); callback=cb_apogee, reltol=1e-8, abstol=1e-8)
+tf_bang, rf_bang = sol_bang2.t[end], sol_bang2[1, end]
+
+println("Bang-bang: r(tf) = ", round(rf_bang, digits=6), "  (t1=", round(t1_bang, digits=4), ", tf=", round(tf_bang, digits=4), ")")
+println("Optimal:   r(tf) = ", round(objective(sol_cold), digits=6), "  (           tf=", round(variable(sol_cold), digits=4), ")")
+
+# assemble the bang-bang trajectory as (t, r, v, m) for plotting
+t_bang = [sol_bang1.t; sol_bang2.t]
+r_bang = [sol_bang1[1, :]; sol_bang2[1, :]]
+
+plt_bang = plot(sol_cold; label="optimal", linewidth=2, color=1)
+plot!(plt_bang[1], t_bang, r_bang; label="bang-bang", linestyle=:dash, linewidth=2, color=2)
+plot(plt_bang[1]; legend=:bottomright, xlabel="time", ylabel="altitude")
+
+using MadNLPGPU
+using CUDA
+
+try
+    global sol_gpu = solve(goddard, :gpu; grid_size=1000, display=false)
+    println("GPU solve succeeded — a functional GPU is available.")
+catch e
+    println("GPU solve failed, as expected without a functional GPU.")
+    println("CUDA.functional() = ", CUDA.functional())
+    println("Exception: ", first(sprint(showerror, e), 400))
+end
+
+using OrdinaryDiffEq   # ODE solver (Hamiltonian flow)
+using NonlinearSolve   # nonlinear equations (shooting)
+
+# maximising control in feedback form
+u_max(x, p) = p[2]
+
+# Hamiltonian flow of the OCP
+φ = Flow(ocp, u_max);
+
+# state projection π(x, p) = x
+proj((x, p)) = x
+
+# shooting function
+S(p0) = proj(φ(t0, x0, p0, tf)) - xf
+
+nle!(s, p0, _) = (s[:] = S(p0))
+
+p_of_t   = costate(direct_sol)     # costate as a function of time
+p0_guess = p_of_t(t0)              # initial costate from the direct method
+
+prob = NonlinearProblem(nle!, p0_guess)
+shooting_sol = solve(prob; show_trace=Val(true))
+p0_sol = shooting_sol.u
+
+println("costate p0 = ", p0_sol)
+println("shoot S(p0) = ", S(p0_sol))
+
+indirect_sol = φ((t0, tf), x0, p0_sol; saveat=range(t0, tf, 100))
+
+plt_compare = plot(direct_sol; label="direct", size=(800, 600))
+plot!(plt_compare, indirect_sol; label="indirect")
+
+# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl
diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md
new file mode 100644
index 00000000..e204eb2f
--- /dev/null
+++ b/docs/src/tutorial.md
@@ -0,0 +1,474 @@
+```@meta
+Draft = false
+```
+
+```@meta
+EditURL = "../src-literate/tutorial.jl"
+```
+
+# OptimalControl.jl — a guided tour
+
+
+This tutorial is a guided tour of [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl), part of the [control-toolbox](https://control-toolbox.org) ecosystem. We follow two problems end to end: a simple **double integrator** for modelling, initialisation and the **indirect** (Pontryagin) method, and the **Goddard rocket** for the **direct** method in depth, grid continuation and GPU solving. Advanced topics are linked at the end.
+
+It is written for readers with a background in optimal control, ODEs or optimisation. By the end you will be able to define an optimal control problem, solve it by both the direct and indirect methods, and visualise the result — all in a few lines of code.
+
+!!! note "Run online"
+    You can run this tutorial interactively in your browser — no installation required — by clicking the Binder badge below:
+
+    [![Binder](https://mybinder.org/badge_logo.svg)](https://mybinder.org/v2/gh/control-toolbox/OptimalControl.jl/paderborn?urlpath=%2Fdoc%2Ftree%2Fdocs%2Fsrc%2Fnotebooks%2Ftutorial.ipynb)
+
+## The problem, and installing the tools
+
+An **optimal control problem** (OCP) in Bolza form reads
+
+```math
+J(x, u) = g(x(t_0), x(t_f)) + \int_{t_0}^{t_f} f^{0}(t, x(t), u(t))\,\mathrm{d}t \;\to\; \min,
+```
+
+subject to the controlled dynamics $\dot{x}(t) = f(t, x(t), u(t))$ and, possibly, box / path / boundary constraints. When $g = 0$ the cost is of **Lagrange** form; when $f^0 = 0$, of **Mayer** form.
+
+More generally, the times $t_0$ and $t_f$ may be free (optimisation variables), and a vector $v$ of additional parameters can enter the cost, dynamics and constraints. The full problem then reads
+
+```math
+\min_{x,u,v}\; g(x(t_0), x(t_f), v) + \int_{t_0}^{t_f} f^{0}(t, x(t), u(t), v)\,\mathrm{d}t,
+```
+
+subject to $\dot{x}(t) = f(t, x(t), u(t), v)$, box / path / boundary constraints.
+
+OptimalControl.jl is the core of the [control-toolbox](https://control-toolbox.org) ecosystem, a modular suite of Julia packages — CTBase (base types & exceptions), CTParser (DSL parsing), CTModels (problem data structures), CTDirect (discretisation & NLP transcription), CTFlows (Hamiltonian flows for indirect methods), and CTSolvers (solver orchestration) — that can also be used individually.
+
+Installation is a single package:
+
+```julia
+import Pkg
+Pkg.add("OptimalControl")
+```
+
+We load OptimalControl.jl to model the problem, a solver backend ([NLPModelsIpopt.jl](https://jso.dev/NLPModelsIpopt.jl)), and [Plots.jl](https://docs.juliaplots.org).
+
+````@example tutorial
+using OptimalControl
+using NLPModelsIpopt
+using Plots
+````
+
+## Defining a problem: `@def` vs macro-free
+
+Our running example: a wagon of unit mass on a frictionless rail, state $x = (q, v)$ (position, velocity), acceleration controlled by a force $u$. We start at $(-1, 0)$, must reach $(0, 0)$ at $t_f = 1$, and minimise the transfer energy
+
+```math
+\frac{1}{2}\int_0^1 u^2(t)\,\mathrm{d}t,
+```
+
+subject to the dynamics
+
+```math
+\dot{q}(t) = v(t), \qquad \dot{v}(t) = u(t).
+```
+
+````@example tutorial
+t0 = 0; tf = 1; x0 = [-1, 0]; xf = [0, 0];
+nothing #hide
+````
+
+### The `@def` macro
+
+The [`@def`](@ref manual-abstract-syntax) macro lets us write the problem almost exactly as the mathematics:
+
+Each line of the `@def` block mirrors a piece of the mathematical formulation — time, state, control, dynamics, boundary conditions, then cost — in the same order one would write them on paper. Unicode symbols (`∈`, `R²`, `ẋ`, `∫`, `→`) make the code read like the maths; plain ASCII alternatives (`R^2`, `derivative`, `integral`, `=>`) are available for keyboards or workflows that prefer them.
+
+````@example tutorial
+ocp = @def begin
+    t ∈ [t0, tf], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+
+    x(t0) == x0
+    x(tf) == xf
+
+    ẋ(t) == [v(t), u(t)]
+
+    0.5∫( u(t)^2 ) → min
+end
+nothing # hide
+````
+
+### The same problem with the macro-free (functional) API
+
+The [functional API](@ref manual-macro-free) builds the *same* model step by step with plain functions — useful for programmatic problem generation or macro-free library code.
+
+````@example tutorial
+pre = OptimalControl.PreModel()
+
+time!(pre; t0=t0, tf=tf)
+state!(pre, 2, "x", ["q", "v"])
+control!(pre, 1)
+
+function f_energy!(dx, t, x, u, v)
+    dx[1] = x[2]
+    dx[2] = u[1]
+    return nothing
+end
+dynamics!(pre, f_energy!)
+
+function boundary_energy!(b, x0_, xf_, v)
+    b[1] = x0_[1] - x0[1]
+    b[2] = x0_[2] - x0[2]
+    b[3] = xf_[1] - xf[1]
+    b[4] = xf_[2] - xf[2]
+    return nothing
+end
+constraint!(pre, :boundary; f=boundary_energy!, lb=zeros(4), ub=zeros(4), label=:endpoint)
+
+lagrange_energy(t, x, u, v) = 0.5 * u[1]^2
+objective!(pre, :min; lagrange=lagrange_energy)
+
+time_dependence!(pre; autonomous=true)
+
+ocp_func = build(pre)
+nothing # hide
+````
+
+### What the macro actually does
+
+**Key message:** `@def` *translates the expression* into the very same functional calls, and **additionally records the symbolic definition**. We can see the difference directly: the macro keeps the DSL expression, whereas the functional API stores an empty definition.
+
+````@example tutorial
+definition(ocp)          # the macro records the full DSL expression
+````
+
+````@example tutorial
+has_abstract_definition(ocp_func)   # false: functional API stores no abstract definition
+````
+
+!!! warning "Two things to keep in mind"
+    - In the functional API, callbacks are **always vector-valued**: even when the control is scalar, one writes `u[1]` — not `u` — inside `f_energy!` or `lagrange_energy`.
+    - The functional API currently works only with the `:adnlp` modeler; it does **not** support the `:exa` modeler needed for GPU solving — one more reason to prefer `@def` when GPU execution is contemplated (more in the GPU section).
+
+## First solve, initial guess, and the costate
+
+Solving is one call, plotting another.
+
+````@example tutorial
+direct_sol = solve(ocp)
+nothing # hide
+````
+
+````@example tutorial
+direct_sol # hide
+````
+
+````@example tutorial
+plot(direct_sol; size=(800, 600))
+````
+
+### The default initial guess
+
+With no initial guess, every variable is initialised to `0.1`. We can *see* the initial guess without optimising, by stopping the solver immediately with `max_iter=0`:
+
+````@example tutorial
+sol_init = solve(ocp; init=nothing, max_iter=0, display=false)
+plot(sol_init; size=(800, 600))
+````
+
+!!! note "Notice the right-hand column: the costate is already there"
+    Even though we only ever provide the state, control and (optional) variable, the solver initialises the **adjoint** internally. After optimisation, this right-column costate is exactly the **adjoint $p$ of Pontryagin's Maximum Principle** — the same $p$ we will reuse to start the indirect method in the indirect section. This closes the loop between the direct and indirect methods.
+
+### Providing our own initial guess
+
+The recommended way to provide an initial guess is the `@init` macro, using the labels from the `@def` block (`q`, `v`, `u` here):
+
+````@example tutorial
+ig = @init ocp begin
+    q(t) := -1 + t
+    v(t) := 0
+    u(t) := 0
+end
+
+sol = solve(ocp; init=ig, display=false)
+println("iterations, default guess: ", iterations(direct_sol))
+println("iterations, @init guess:   ", iterations(sol))
+````
+
+In this case both guesses give **1 iteration**: the double integrator is a *linear-quadratic* problem, so the NLP is quadratic and Ipopt solves it in a single step regardless of the starting point. Warm-starting only pays off on genuinely nonlinear problems — we will see this with the **Goddard rocket** in the next section.
+
+For all the ways to specify an initial guess, see [Set an initial guess](@ref manual-initial-guess).
+!!! note
+    There is currently no way to initialise the costate directly — only state, control and variable can be provided through `@init`. The solver initialises the adjoint internally (as we saw above). Costate initialisation is a planned feature.
+
+## Direct method in depth: Goddard
+
+### Discretise optimal control problems
+
+The **direct** method turns the infinite-dimensional OCP into a finite-dimensional nonlinear program (NLP) by discretising time (Runge–Kutta / collocation) on a grid, then hands the NLP to a solver. It is robust and easy to use.
+
+Concretely, time is discretised on a uniform grid $t_0 < t_1 < \dots < t_N = t_f$ with step $h = (t_f - t_0)/N$. The (explicit) Euler scheme, for instance, replaces the dynamics by
+
+```math
+x_{i} = x_{i-1} + h\,f(t_{i-1}, x_{i-1}, u_{i-1}), \quad i = 1, \dots, N,
+```
+
+and the integral cost by the corresponding rectangle sum
+
+```math
+h\sum_{i=0}^{N-1} f^{0}(t_i, x_i, u_i).
+```
+
+The continuous OCP thus becomes a finite-dimensional NLP in the variables $X = (x_0, \dots, x_N, u_0, \dots, u_N)$, which is passed to an NLP solver such as [Ipopt](https://coin-or.github.io/Ipopt). Higher-order schemes (midpoint, Gauss–Legendre collocation) follow the same principle with different quadrature and interpolation formulas — `solve` defaults to the second-order `:midpoint` scheme, not Euler.
+
+### The Goddard rocket problem
+
+To demonstrate convergence behaviour and warm-starting, we need a genuinely nonlinear problem. The **Goddard rocket** — maximise the final altitude, with free final time and a singular arc — is a classic test case.
+
+````@example tutorial
+# Goddard data and dynamics (F0: drift, F1: thrust)
+const r0 = 1
+const v0 = 0
+const m0 = 1
+const mf = 0.6
+const Cd = 310
+const Tmax = 3.5
+const β = 500
+const b = 2
+
+F0(x) = begin
+    r, v, m = x
+    D = Cd * v^2 * exp(-β * (r - 1))
+    [v, -D/m - 1/r^2, 0]
+end
+F1(x) = begin
+    r, v, m = x
+    [0, Tmax/m, -b*Tmax]
+end
+
+goddard = @def begin
+    tf ∈ R, variable
+    t ∈ [t0, tf], time
+    x = (r, v, m) ∈ R³, state
+    u ∈ R, control
+
+    x(t0) == [r0, v0, m0]
+    m(tf) == mf
+    0 ≤ u(t) ≤ 1
+    r(t) ≥ r0
+
+    ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
+
+    r(tf) → max
+end
+nothing # hide
+````
+
+### Choosing a solver is trivial
+
+`solve` uses the defaults (collocation, ADNLP modeler, Ipopt, CPU). Switching solver is just loading a package and passing a token (see [Solve a problem](@ref manual-solve)):
+
+````@example tutorial
+using MadNLP
+
+sol_ipopt  = solve(goddard;          grid_size=250, display=false)
+sol_madnlp = solve(goddard, :madnlp; grid_size=250, display=false)
+
+println("Ipopt  : r(tf) = ", objective(sol_ipopt),  ", ", iterations(sol_ipopt),  " iters")
+println("MadNLP : r(tf) = ", objective(sol_madnlp), ", ", iterations(sol_madnlp), " iters")
+````
+
+The available methods and their options can be inspected with `methods()` and `describe(:collocation)`; we will not dwell on them here.
+
+### Grid continuation by warm-starting
+
+A solution can be passed **directly** as the initial guess of another solve — it is interpolated onto the new grid. This makes discrete continuation trivial and ties back to the initialisation above. On this nonlinear problem it genuinely **pays**: we compare reaching a fine grid of 1000 two ways:
+
+1. **cold start** — solve `grid_size=1000` directly;
+2. **cascade** — solve `grid_size=50` first, then `grid_size=1000` warm-started with that solution.
+
+````@example tutorial
+# solutions computed once, reused for iteration counts and the overlay plot
+sol_cold = solve(goddard; grid_size=1000, display=false)
+
+# warm cascade: grid 50 first, then grid 1000 initialised from it
+s50   = solve(goddard; grid_size=50, display=false)
+s1000 = solve(goddard; grid_size=1000, init=s50, display=false)
+
+println("cold    grid 1000        : ", iterations(sol_cold), " iters")
+println("cascade grid 50 (warm-up): ", iterations(s50),      " iters")
+println("cascade grid 1000 (warm) : ", iterations(s1000),    " iters")
+````
+
+**Message:** what matters is the iteration count *at the expensive grid* — the warm-started `iterations(s1000)` is well below the cold `iterations(sol_cold)`, even though the cheap `grid_size=50` warm-up adds iterations of its own to the running total; since a grid-50 iteration is far cheaper than a grid-1000 iteration, the cascade still wins on wall-clock time. Overlay the successive solutions to watch convergence:
+
+````@example tutorial
+plt = plot(s50;  label="50", size=(800, 800))
+plot!(plt, s1000; label="1000")
+````
+
+This is grid-refinement warm-starting. The very same mechanism drives **parametric** continuation (homotopy on a physical parameter, e.g. maximum thrust): [Discrete continuation](@extref Tutorials tutorial-continuation).
+
+### Comparison with a bang-bang strategy
+
+How much better is the optimal solution compared to a naive strategy? We simulate **full thrust until fuel depletion, then coast to apogee** — a bang-bang profile with no optimisation, just two ODE integrations with callbacks.
+
+````@example tutorial
+using OrdinaryDiffEq   # ODE solver (callbacks for the bang-bang simulation)
+
+# Phase 1: u = 1, stop when m = mf (fuel depleted)
+bang1!(dx, x, p, t) = (dx[:] = F0(x) + F1(x))
+cb_fuel = ContinuousCallback((u, t, int) -> u[3] - mf, terminate!)
+sol_bang1 = solve(ODEProblem(bang1!, [r0, v0, m0], (t0, 100.0)), Tsit5(); callback=cb_fuel, reltol=1e-8, abstol=1e-8)
+t1_bang, x1_bang = sol_bang1.t[end], sol_bang1[:, end]
+
+# Phase 2: u = 0, stop when v = 0 (apogee)
+bang2!(dx, x, p, t) = (dx[:] = F0(x))
+cb_apogee = ContinuousCallback((u, t, int) -> u[2], terminate!)
+sol_bang2 = solve(ODEProblem(bang2!, x1_bang, (t1_bang, 1000.0)), Tsit5(); callback=cb_apogee, reltol=1e-8, abstol=1e-8)
+tf_bang, rf_bang = sol_bang2.t[end], sol_bang2[1, end]
+
+println("Bang-bang: r(tf) = ", round(rf_bang, digits=6), "  (t1=", round(t1_bang, digits=4), ", tf=", round(tf_bang, digits=4), ")")
+println("Optimal:   r(tf) = ", round(objective(sol_cold), digits=6), "  (           tf=", round(variable(sol_cold), digits=4), ")")
+````
+
+The optimal thrust profile uses a **singular arc** — it does not simply push at the maximum. Overlaying the two trajectories on the altitude–velocity plane makes the difference visible:
+
+````@example tutorial
+# assemble the bang-bang trajectory as (t, r, v, m) for plotting
+t_bang = [sol_bang1.t; sol_bang2.t]
+r_bang = [sol_bang1[1, :]; sol_bang2[1, :]]
+
+plt_bang = plot(sol_cold; label="optimal", linewidth=2, color=1)
+plot!(plt_bang[1], t_bang, r_bang; label="bang-bang", linestyle=:dash, linewidth=2, color=2)
+plot(plt_bang[1]; legend=:bottomright, xlabel="time", ylabel="altitude")
+````
+
+## Solving on a GPU
+
+Moving to the GPU is a single token, `:gpu`, which auto-completes to `(:collocation, :exa, :madnlp, :gpu)`. It requires the `:exa` modeler (hence `@def`, not the macro-free API — cf. the definition section) plus a CUDA-capable GPU.
+
+In a seminar or on Binder there is usually **no functional GPU**, so the call is *expected to fail* — that is the pedagogical point: the `:gpu` token needs a specific setup. We wrap it in a `try/catch` so the tour keeps running and shows the raised exception.
+
+````@example tutorial
+using MadNLPGPU
+using CUDA
+
+try
+    global sol_gpu = solve(goddard, :gpu; grid_size=1000, display=false)
+    println("GPU solve succeeded — a functional GPU is available.")
+catch e
+    println("GPU solve failed, as expected without a functional GPU.")
+    println("CUDA.functional() = ", CUDA.functional())
+    println("Exception: ", first(sprint(showerror, e), 400))
+end
+````
+
+For the full GPU setup, see [Solve on GPU](@ref manual-solve-gpu).
+
+## The indirect method
+
+We now return to the **double integrator** `ocp` from the earlier sections. Its shooting has just two unknowns and is initialised by the direct costate above, which makes it ideal to *see* the indirect method. (The Goddard shooting is a *structured multi-arc* problem — see the links in the last section.)
+
+In control-toolbox we systematically pair the direct method with the **indirect** one, based on Pontryagin's Maximum Principle (PMP), with pseudo-Hamiltonian
+
+```math
+H(x,p,u) = p\,f(x,u) + p^0 f^0(x,u) \qquad (\text{normal case } p^0 = -1).
+```
+
+The PMP gives the maximising control in feedback form
+
+```math
+u(x,p) = \arg\max_u H,
+```
+
+and the optimal trajectory solves a boundary value problem that we recast as a **shooting equation**
+
+```math
+S(p_0) = 0.
+```
+
+The indirect method proceeds in three steps:
+
+1. **Maximising control.** The PMP yields the control in feedback form $u(x, p) = \arg\max_u H(x, p, u)$. Substituting back gives the maximised Hamiltonian
+
+   ```math
+   \mathbf{H}(x, p) = H(x, p, u(x, p)).
+   ```
+
+2. **Boundary value problem.** The optimal trajectory satisfies the Hamiltonian system
+
+   ```math
+   \dot{x} = \nabla_p \mathbf{H}, \qquad \dot{p} = -\nabla_x \mathbf{H},
+   ```
+
+   with boundary conditions $x(t_0) = x_0$, $x(t_f) = x_f$.
+
+3. **Shooting function.** Let $\varphi_{t_0, x_0, p_0}(\cdot)$ denote the flow of the Hamiltonian vector field from $(x_0, p_0)$. The shooting function
+
+   ```math
+   S(p_0) = \pi(\varphi_{t_0, x_0, p_0}(t_f)) - x_f, \qquad \pi(x, p) = x,
+   ```
+
+   measures the miss at $t_f$: solving the BVP reduces to finding $p_0$ such that $S(p_0) = 0$.
+
+For the energy problem, $H = p_1 v + p_2 u - u^2/2$, so the maximiser is $u = p_2$.
+
+````@example tutorial
+using OrdinaryDiffEq   # ODE solver (Hamiltonian flow)
+using NonlinearSolve   # nonlinear equations (shooting)
+
+# maximising control in feedback form
+u_max(x, p) = p[2]
+
+# Hamiltonian flow of the OCP
+φ = Flow(ocp, u_max);
+
+# state projection π(x, p) = x
+proj((x, p)) = x
+
+# shooting function
+S(p0) = proj(φ(t0, x0, p0, tf)) - xf
+nothing # hide
+````
+
+**The shooting is initialised with the costate of the direct solution** — the very adjoint we highlighted above:
+
+````@example tutorial
+nle!(s, p0, _) = (s[:] = S(p0))
+
+p_of_t   = costate(direct_sol)     # costate as a function of time
+p0_guess = p_of_t(t0)              # initial costate from the direct method
+
+prob = NonlinearProblem(nle!, p0_guess)
+shooting_sol = solve(prob; show_trace=Val(true))
+p0_sol = shooting_sol.u
+
+println("costate p0 = ", p0_sol)
+println("shoot S(p0) = ", S(p0_sol))
+````
+
+Reconstruct the indirect solution from the flow and overlay it with the direct solution:
+
+````@example tutorial
+indirect_sol = φ((t0, tf), x0, p0_sol; saveat=range(t0, tf, 100))
+
+plt_compare = plot(direct_sol; label="direct", size=(800, 600))
+plot!(plt_compare, indirect_sol; label="indirect")
+````
+
+See [Compute flows from optimal control problems](@ref manual-flow-ocp) for the flow construction, and the [indirect simple shooting tutorial](@extref tutorial-indirect-simple-shooting).
+
+## Going further
+
+**Variables & parameters.** Beyond the control, one can optimise **parameters** naturally, both in an OCP (the `variable` keyword of the DSL) and in a differential-constraint optimisation problem **without any control** (a *control-free* problem).
+See [control-free problems](@ref example-control-free).
+
+**Advanced examples** (each does both direct and indirect):
+
+- Singular control (control-affine systems) — [singular control](@ref example-singular-control)
+- State constraint — [state constraint](@ref example-state-constraint)
+- Goddard problem — free final time, a singular arc, a state constraint and a structured shooting all at once — [Goddard tutorial](@extref Tutorials tutorial-goddard)
+
+**Discrete continuation** — warm-starting across a family of problems (homotopy on a physical parameter), the grown-up version of the grid continuation above: [Discrete continuation](@extref Tutorials tutorial-continuation).
+
+---
+
+*This page was generated using [Literate.jl](https://github.com/fredrikekre/Literate.jl).*
+