diff --git a/BREAKING.md b/BREAKING.md
index 4f442727b..745577b3a 100644
--- a/BREAKING.md
+++ b/BREAKING.md
@@ -10,7 +10,8 @@ Version 2.0.0 represents a major architectural redesign of OptimalControl.jl, in
 - **GPU/CPU parameter system** for heterogeneous computing
 - **Advanced option routing** with introspection and disambiguation tools
 - **New solver integrations** (Uno, MadNCL)
-- **Control-free problems** support
+- **Control-free problems** support with augmented Hamiltonian approach
+- **CTFlows enhancements** with `augment=true` and direct OCP flow creation
 - **Modernized reexport system** using `@reexport import`
 
 ## Removed Functions
@@ -172,7 +173,7 @@ The old functional approach is no longer supported.
 
 These features are new in v2.0.0 but don't break existing code:
 
-### Control-Free Problems
+### Control-Free Problems Support
 
 Support for optimal control problems without control variables:
 
@@ -220,6 +221,45 @@ using CUDA, MadNLPGPU
 sol = solve(ocp, :collocation, :exa, :madnlp, :gpu)
 ```
 
+## CTFlows Features
+
+### Control-Free Problems
+
+v2.0.0 introduces comprehensive support for control-free problems (optimal control without control variables) with enhanced CTFlows integration:
+
+**Augmented Hamiltonian approach:**
+
+```julia
+# v1.1.6: Manual augmented Hamiltonian construction
+function H_aug(t, x_, p_)
+    x, λ = x_
+    p, _ = p_
+    return H(t, x, p, λ)
+end
+f = Flow(Hamiltonian(H_aug))
+
+# v2.0.0: Direct OCP flow creation
+f = Flow(ocp)
+```
+
+**Automatic costate computation:**
+
+```julia
+# v2.0.0: augment=true automatically computes p_λ(tf)
+function shoot!(s, p0, λ)
+    _, px_tf, pλ_tf = f(t0, x0, p0, tf, λ; augment=true)
+    s[1] = px_tf  # p(tf) = 0
+    s[2] = pλ_tf  # p_λ(tf) = 0
+end
+```
+
+**Mathematical framework:**
+
+- Complete augmented system dynamics with proper transversality conditions
+- Automatic handling of Lagrange costs: $p_\lambda(t_f) = 0$
+- Automatic handling of Mayer costs: $p_\omega(t_f) = -2\omega$
+- Initial conditions: $p_\lambda(t_0) = 0$ by construction
+
 ## Dependency Updates
 
 v2.0.0 requires updated versions of CTX packages:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index bb92cb390..3594a6a1c 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -55,6 +55,10 @@ Versions follow [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   - Optimal control problems without control variables
   - Optimization of constant parameters in dynamical systems
   - Full integration with solve pipeline
+  - **Augmented Hamiltonian approach**: `augment=true` feature in CTFlows for automatic costate computation
+  - **Simplified flow creation**: `Flow(ocp)` directly creates Hamiltonian flow from control-free problems
+  - **Mathematical framework**: Complete transversality conditions for variable parameters
+  - **Documentation**: Comprehensive examples with exponential growth and harmonic oscillator
 
 - **New solvers**:
   - **Uno**: CPU-only nonlinear optimization solver (methods with `:uno`)
@@ -78,6 +82,12 @@ Versions follow [Semantic Versioning](https://semver.org/spec/v2.0.0.html).
   - Organized by source package (ctbase.jl, ctdirect.jl, ctflows.jl, ctmodels.jl, ctparser.jl, ctsolvers.jl)
   - Cleaner separation between imported and exported symbols
 
+- **CTFlows enhancements**:
+  - **Augmented Hamiltonian computation**: `augment=true` automatically computes costates for variable parameters
+  - **Direct OCP flow creation**: `Flow(ocp)` creates Hamiltonian flow without manual Hamiltonian definition
+  - **Transversality conditions**: Automatic handling of $p_\lambda(t_f) = 0$ for Lagrange costs and $p_\omega(t_f) = -2\omega$ for Mayer costs
+  - **Mathematical rigor**: Complete augmented system dynamics with proper initial conditions
+
 - **Strategy registry system**:
   - `StrategyRegistry` with metadata for all strategies
   - `StrategyMetadata` with id, options, and parameter support
diff --git a/docs/src/assets/Manifest.toml b/docs/src/assets/Manifest.toml
index 0f77cd71e..514f44a3e 100644
--- a/docs/src/assets/Manifest.toml
+++ b/docs/src/assets/Manifest.toml
@@ -232,9 +232,9 @@ version = "1.0.6"
 
 [[deps.CTFlows]]
 deps = ["CTBase", "CTModels", "DocStringExtensions", "ForwardDiff", "LinearAlgebra", "MLStyle", "MacroTools"]
-git-tree-sha1 = "536c34374ab27f3d4d6a51d9e92ed89f146ca989"
+git-tree-sha1 = "c4a9adf514aca4e006131a5e39792bcf01b449b4"
 uuid = "1c39547c-7794-42f7-af83-d98194f657c2"
-version = "0.8.20"
+version = "0.8.22-beta"
 weakdeps = ["OrdinaryDiffEq"]
 
     [deps.CTFlows.extensions]
@@ -242,9 +242,9 @@ weakdeps = ["OrdinaryDiffEq"]
 
 [[deps.CTModels]]
 deps = ["CTBase", "DocStringExtensions", "LinearAlgebra", "MLStyle", "MacroTools", "OrderedCollections", "Parameters", "RecipesBase"]
-git-tree-sha1 = "3e9df6b6cb96ccb05051ded9a5d4f76f649f6d0c"
+git-tree-sha1 = "791bcfa3cc9b39d2c1b5bc66581d1739bd8c96c0"
 uuid = "34c4fa32-2049-4079-8329-de33c2a22e2d"
-version = "0.9.12-beta"
+version = "0.9.13"
 weakdeps = ["JLD2", "JSON3", "Plots"]
 
     [deps.CTModels.extensions]
@@ -254,9 +254,9 @@ weakdeps = ["JLD2", "JSON3", "Plots"]
 
 [[deps.CTParser]]
 deps = ["CTBase", "DocStringExtensions", "MLStyle", "OrderedCollections", "Parameters", "Unicode"]
-git-tree-sha1 = "1a7445af1c937c291371e0f7c0a51a80f5b5d0eb"
+git-tree-sha1 = "0a57111b2d95d6272bb76f4f6f7c4fee1e40ce82"
 uuid = "32681960-a1b1-40db-9bff-a1ca817385d1"
-version = "0.8.12-beta"
+version = "0.8.13"
 
 [[deps.CTSolvers]]
 deps = ["ADNLPModels", "CTBase", "CTModels", "CommonSolve", "DocStringExtensions", "ExaModels", "KernelAbstractions", "NLPModels", "SolverCore"]
@@ -921,9 +921,16 @@ version = "1.1.3"
 
 [[deps.FunctionWrappersWrappers]]
 deps = ["FunctionWrappers", "PrecompileTools", "TruncatedStacktraces"]
-git-tree-sha1 = "6874da243fb93e34201d7d4587ffa0e920682f64"
+git-tree-sha1 = "5201523536a43bf8aef3914b7f60b552b098ef8e"
 uuid = "77dc65aa-8811-40c2-897b-53d922fa7daf"
-version = "1.0.0"
+version = "1.1.0"
+
+    [deps.FunctionWrappersWrappers.extensions]
+    FunctionWrappersWrappersEnzymeExt = ["Enzyme", "EnzymeCore"]
+
+    [deps.FunctionWrappersWrappers.weakdeps]
+    Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
+    EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
 
 [[deps.Future]]
 deps = ["Random"]
@@ -1696,21 +1703,21 @@ version = "0.21.12"
 
 [[deps.NLPModelsIpopt]]
 deps = ["Ipopt", "NLPModels", "NLPModelsModifiers", "SolverCore"]
-git-tree-sha1 = "37a2187dc5f1291f6a6cb592acaa59cee685b5fb"
+git-tree-sha1 = "260a33809f49f7ccb011fdc1d9ebcf8879bc9757"
 uuid = "f4238b75-b362-5c4c-b852-0801c9a21d71"
-version = "0.11.2"
+version = "0.11.3"
 
 [[deps.NLPModelsKnitro]]
 deps = ["KNITRO", "NLPModels", "NLPModelsModifiers", "SolverCore"]
-git-tree-sha1 = "90e9e12bf859aa3b9453a852015a552592b8ea95"
+git-tree-sha1 = "763bda401fa32460dae1caa48f6ab3c195d3bf3e"
 uuid = "bec4dd0d-7755-52d5-9a02-22f0ffc7efcb"
-version = "0.10.2"
+version = "0.10.3"
 
 [[deps.NLPModelsModifiers]]
 deps = ["FastClosures", "LinearAlgebra", "LinearOperators", "NLPModels", "Printf", "SparseArrays"]
-git-tree-sha1 = "15bd5325aaa45090019f43b49cb3cb938b195f04"
+git-tree-sha1 = "3a6ec5add6f5a95d598648432ac06126157ba11c"
 uuid = "e01155f1-5c6f-4375-a9d8-616dd036575f"
-version = "0.7.4"
+version = "0.8.0"
 
 [[deps.NLSolversBase]]
 deps = ["ADTypes", "DifferentiationInterface", "FiniteDiff", "LinearAlgebra"]
@@ -1779,9 +1786,9 @@ version = "4.17.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 = "7deb924291e30ef27e8823e9d048ffb98ca6ffce"
+git-tree-sha1 = "e9bad06173df78880b8c0fd51ad9cd940547a195"
 uuid = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
-version = "2.20.0"
+version = "2.21.0"
 
     [deps.NonlinearSolveBase.extensions]
     NonlinearSolveBaseBandedMatricesExt = "BandedMatrices"
@@ -1811,15 +1818,15 @@ version = "2.20.0"
 
 [[deps.NonlinearSolveFirstOrder]]
 deps = ["ADTypes", "ArrayInterface", "CommonSolve", "ConcreteStructs", "FiniteDiff", "ForwardDiff", "LineSearch", "LinearAlgebra", "LinearSolve", "MaybeInplace", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase", "SciMLJacobianOperators", "Setfield", "StaticArraysCore"]
-git-tree-sha1 = "eea7cbe389b168c77df7ff779fb7277019c685c8"
+git-tree-sha1 = "ae36a13005343aed5e731f146637f61a54702444"
 uuid = "5959db7a-ea39-4486-b5fe-2dd0bf03d60d"
-version = "2.0.0"
+version = "2.1.0"
 
 [[deps.NonlinearSolveQuasiNewton]]
 deps = ["ArrayInterface", "CommonSolve", "ConcreteStructs", "LinearAlgebra", "LinearSolve", "MaybeInplace", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase", "SciMLOperators", "StaticArraysCore"]
-git-tree-sha1 = "ade27e8e9566b6cec63ee62f6a6650a11cf9a2eb"
+git-tree-sha1 = "2ada8e9892fd60b3f01a4984704d66fca682bd50"
 uuid = "9a2c21bd-3a47-402d-9113-8faf9a0ee114"
-version = "1.12.0"
+version = "1.13.0"
 weakdeps = ["ForwardDiff"]
 
     [deps.NonlinearSolveQuasiNewton.extensions]
@@ -1827,9 +1834,9 @@ weakdeps = ["ForwardDiff"]
 
 [[deps.NonlinearSolveSpectralMethods]]
 deps = ["CommonSolve", "ConcreteStructs", "LineSearch", "MaybeInplace", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase"]
-git-tree-sha1 = "eafd027b5cd768f19bb5de76c0e908a9065ddd36"
+git-tree-sha1 = "e3e8be34968ad07be8734c3d6e42b3f14218b694"
 uuid = "26075421-4e9a-44e1-8bd1-420ed7ad02b2"
-version = "1.6.0"
+version = "1.7.0"
 weakdeps = ["ForwardDiff"]
 
     [deps.NonlinearSolveSpectralMethods.extensions]
@@ -1931,9 +1938,9 @@ version = "1.14.0"
 
 [[deps.OrdinaryDiffEqDifferentiation]]
 deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqCore", "SciMLBase", "SciMLOperators", "SparseMatrixColorings", "StaticArrays"]
-git-tree-sha1 = "c7d06493d3327cc45f0d8d7e641203d28779e837"
+git-tree-sha1 = "8b18ce8ab49cd04a1443fd63eb569d270a07090c"
 uuid = "4302a76b-040a-498a-8c04-15b101fed76b"
-version = "2.6.0"
+version = "2.7.0"
 weakdeps = ["SparseArrays"]
 
     [deps.OrdinaryDiffEqDifferentiation.extensions]
@@ -2007,9 +2014,9 @@ version = "1.14.0"
 
 [[deps.OrdinaryDiffEqNonlinearSolve]]
 deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "FastClosures", "ForwardDiff", "LinearAlgebra", "LinearSolve", "MuladdMacro", "NonlinearSolve", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "PreallocationTools", "RecursiveArrayTools", "SciMLBase", "SciMLOperators", "SciMLStructures", "SimpleNonlinearSolve", "SparseArrays", "StaticArrays"]
-git-tree-sha1 = "f6722d6a96c27263ac3ea6159a3174914567a807"
+git-tree-sha1 = "7eabb9877144d3cfc894fee899a76439ee2070de"
 uuid = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
-version = "1.25.0"
+version = "1.26.0"
 
 [[deps.OrdinaryDiffEqNordsieck]]
 deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqTsit5", "Polyester", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"]
@@ -2043,9 +2050,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 = "88f722a74df3dadaf5f31c454a280f1201a1b368"
+git-tree-sha1 = "25e5043968c80e0addead2d7533d30a608478842"
 uuid = "43230ef6-c299-4910-a778-202eb28ce4ce"
-version = "1.28.0"
+version = "1.29.0"
 
 [[deps.OrdinaryDiffEqSDIRK]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"]
@@ -2289,12 +2296,13 @@ version = "0.6.12"
 
 [[deps.RecursiveArrayTools]]
 deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "LinearAlgebra", "PrecompileTools", "RecipesBase", "StaticArraysCore", "SymbolicIndexingInterface"]
-git-tree-sha1 = "e4fd3369c78666a65ccec25dba28a0b181434ab2"
+git-tree-sha1 = "d0282d612f22dcad7b81cf487b746e63aa2a6709"
 uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
-version = "3.52.0"
+version = "3.54.0"
 
     [deps.RecursiveArrayTools.extensions]
     RecursiveArrayToolsFastBroadcastExt = "FastBroadcast"
+    RecursiveArrayToolsFastBroadcastPolyesterExt = ["FastBroadcast", "Polyester"]
     RecursiveArrayToolsForwardDiffExt = "ForwardDiff"
     RecursiveArrayToolsKernelAbstractionsExt = "KernelAbstractions"
     RecursiveArrayToolsMeasurementsExt = "Measurements"
@@ -2313,6 +2321,7 @@ version = "3.52.0"
     KernelAbstractions = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
     Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
     MonteCarloMeasurements = "0987c9cc-fe09-11e8-30f0-b96dd679fdca"
+    Polyester = "f517fe37-dbe3-4b94-8317-1923a5111588"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
     SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
     Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
@@ -2434,9 +2443,9 @@ weakdeps = ["Tracy"]
 
 [[deps.SciMLOperators]]
 deps = ["Accessors", "ArrayInterface", "DocStringExtensions", "LinearAlgebra"]
-git-tree-sha1 = "794c760e6aafe9f40dcd7dd30526ea33f0adc8b7"
+git-tree-sha1 = "234869cf9fee9258a95464b7a7065cc7be84db00"
 uuid = "c0aeaf25-5076-4817-a8d5-81caf7dfa961"
-version = "1.15.1"
+version = "1.16.0"
 weakdeps = ["SparseArrays", "StaticArraysCore"]
 
     [deps.SciMLOperators.extensions]
diff --git a/docs/src/example-control-free.md b/docs/src/example-control-free.md
index d4023b3a3..29633349b 100644
--- a/docs/src/example-control-free.md
+++ b/docs/src/example-control-free.md
@@ -9,7 +9,7 @@ This page demonstrates two simple examples with known analytical solutions.
 
 First, we import the necessary packages:
 
-```@example main
+```@example main-growth
 using OptimalControl
 using NLPModelsIpopt
 using Plots
@@ -20,55 +20,166 @@ using Plots
 Consider a system with exponential growth:
 
 ```math
-\dot{x}(t) = p \cdot x(t), \quad x(0) = 2
+\dot{x}(t) = \lambda \cdot x(t), \quad x(0) = 2
 ```
 
-where $p$ is an unknown growth rate parameter. We have observed data following $x_{\text{obs}}(t) = 2 e^{0.5 t}$ and want to estimate $p$ by minimizing the squared error:
+where $\lambda$ is an unknown growth rate parameter. We have observed data with some perturbations and want to estimate $\lambda$ by minimizing the squared error:
 
 ```math
-\min_{p} \int_0^{10} (x(t) - x_{\text{obs}}(t))^2 \, dt
+\min_{\lambda} \int_0^{2} (x(t) - x_{\text{obs}}(t))^2 \, \mathrm{d}t
 ```
 
-The analytical solution is $p = 0.5$.
+The underlying model has $\lambda = 0.5$, but the observed data includes perturbations.
 
 ### [Problem definition](@id example-control-free-problem-1)
 
-```@example main
-# observed data (analytical solution)
-data(t) = 2.0 * exp(0.5 * t)
+```@example main-growth
+# observed data (analytical solution with λ = 0.5)
+λ_true = 0.5
+model(t) = 2 * exp(λ_true * t)
+perturbation(t) = 2e-1*sin(4π*t)
+data(t) = model(t) + perturbation(t)
 
 # optimal control problem (parameter estimation)
-ocp1 = @def begin
-    p ∈ R, variable              # growth rate to estimate
-    t ∈ [0, 10], time
+t0 = 0; tf = 2; x0 = 2
+ocp = @def begin
+    λ ∈ R, variable              # growth rate to estimate
+    t ∈ [t0, tf], time
     x ∈ R, state
     
-    x(0) == 2.0
-    ẋ(t) == p * x(t)
+    x(t0) == x0
+    ẋ(t) == λ * x(t)
     
     ∫((x(t) - data(t))^2) → min  # fit to observed data
 end
 nothing # hide
 ```
 
-### [Solution](@id example-control-free-solution-1)
+### Direct method
 
-```@example main
-sol1 = solve(ocp1; display=false)
-println("Estimated growth rate: p = ", variable(sol1))
-println("Objective value: ", objective(sol1))
-println("Expected: p = 0.5, objective ≈ 0.0")
+```@example main-growth
+direct_sol = solve(ocp; grid_size=20, display=false)
+println("Estimated growth rate: λ = ", variable(direct_sol))
+println("Objective value: ", objective(direct_sol))
 nothing # hide
 ```
 
-```@example main
-plot(sol1; size=(800, 400))
+```@example main-growth
+# plot direct solution
+plt = plot(direct_sol; size=(800, 400), label="Direct")
+
+# Add data on first plot
+t_grid = time_grid(direct_sol)
+plot!(plt, t_grid, data.(t_grid); subplot=1, line=:dot, lw=2, label="Data", color=:black)
+```
+
+The estimated parameter should be close to $\lambda \approx 0.5$.
+
+### Indirect method
+
+We now solve the same problem using an indirect shooting method based on Pontryagin's Maximum Principle. First, we import the necessary packages:
+
+```@example main-growth
+using OrdinaryDiffEq  # ODE solver
+using NonlinearSolve  # Nonlinear solver
 ```
 
-The estimated parameter should be very close to $p = 0.5$, and the objective (squared error) should be nearly zero since we're fitting to the exact analytical solution.
+For control-free problems with a variable parameter, we use an **augmented Hamiltonian** approach. The Hamiltonian for this problem is:
+
+```math
+H(t, x, p, \lambda) = p \lambda x - (x - x_{\text{obs}}(t))^2
+```
+
+To handle the variable parameter $\lambda$, we treat it as an additional state with zero dynamics. This gives us the augmented system with state $(x, \lambda)$ and costate $(p, p_\lambda)$, where:
+
+```math
+\begin{aligned}
+\frac{\mathrm{d}x}{\mathrm{d}t} &= \frac{\partial H}{\partial p} = \lambda x \\
+\frac{\mathrm{d}\lambda}{\mathrm{d}t} &= 0 \quad \text{(constant parameter)} \\
+\frac{\mathrm{d}p}{\mathrm{d}t} &= -\frac{\partial H}{\partial x} = -p\lambda + 2(x - x_{\text{obs}}(t)) \\
+\frac{\mathrm{d}p_\lambda}{\mathrm{d}t} &= -\frac{\partial H}{\partial \lambda} = -px
+\end{aligned}
+```
+
+The transversality condition for the variable parameter requires $p_\lambda(t_f) - p_\lambda(t_0) = 0$. Assuming $p_\lambda(t_0) = 0$, we have to satisfy:
+
+```math
+p_\lambda(t_f) = -\int_{t_0}^{t_f} \frac{\partial H}{\partial \lambda}(t, x(t), p(t), \lambda) \, \mathrm{d}t = 0
+```
+
+We use CTFlows' `augment=true` feature to automatically compute $p_\lambda(t_f)$ without manually constructing the augmented system.
+
+```@example main-growth
+# Create Hamiltonian flow from OCP
+f = Flow(ocp)
+nothing # hide
+```
+
+!!! note
+
+    For more details about the flow construction, see [this page](@ref manual-flow-others).
+
+The shooting function enforces the transversality conditions $p(t_f) = 0$ and $p_\lambda(t_f) = 0$. Using `augment=true`, the flow automatically returns $(x(t_f), p(t_f), p_\lambda(t_f))$, with $p_\lambda(t_0) = 0$ by construction.
+
+```@example main-growth
+# Shooting function: S(p0, λ) = (p(tf), pλ(tf))
+# We want both components to be zero at tf
+function shoot!(s, p0, λ)
+    _, px_tf, pλ_tf = f(t0, x0, p0, tf, λ; augment=true)
+    s[1] = px_tf
+    s[2] = pλ_tf
+    return nothing
+end
+
+# Auxiliary in-place NLE function
+nle!(s, y, _) = shoot!(s, y...)
+nothing # hide
+```
+
+We use the direct solution to initialize the shooting method:
+
+```@example main-growth
+# Extract solution from direct method for initialization
+p_direct = costate(direct_sol)
+λ_direct = variable(direct_sol)
+
+# Initial guess
+p0_guess = p_direct(t0)
+λ_guess = λ_direct
+
+# NLE problem with initial guess (2 unknowns: p0, λ)
+prob_indirect = NonlinearProblem(nle!, [p0_guess, λ_guess])
+
+# Solve shooting equations
+shooting_sol = solve(prob_indirect; show_trace=Val(false))
+p0_sol, λ_sol = shooting_sol.u
+
+println("Indirect solution:")
+println("Initial costate: p0 = ", p0_sol)
+println("Parameter: λ = ", λ_sol)
+nothing # hide
+```
+
+Finally, we compute and plot the indirect solution:
+
+```@example main-growth
+# Compute and plot indirect solution
+indirect_sol = f((t0, tf), x0, p0_sol, λ_sol; saveat=range(t0, tf, 200))
+plot!(plt, indirect_sol; linestyle=:dash, lw=2, label="Indirect", color=2)
+```
+
+The direct and indirect solutions match closely, both fitting the perturbed observed data.
 
 ## Example 2: Harmonic oscillator pulsation optimization
 
+```@setup main-harmonic
+using OptimalControl
+using NLPModelsIpopt
+using Plots
+using OrdinaryDiffEq  # ODE solver
+using NonlinearSolve  # Nonlinear solver
+```
+
 Consider a harmonic oscillator:
 
 ```math
@@ -91,16 +202,18 @@ The analytical solution is $\omega = \pi/2 \approx 1.5708$, giving $q(t) = \cos(
 
 ### [Problem definition](@id example-control-free-problem-2)
 
-```@example main
+```@example main-harmonic
 # optimal control problem (pulsation optimization)
-ocp2 = @def begin
+q0 = 1; v0 = 0
+t0 = 0; tf = 1
+ocp = @def begin
     ω ∈ R, variable              # pulsation to optimize
-    t ∈ [0, 1], time
+    t ∈ [t0, tf], time
     x = (q, v) ∈ R², state
     
-    q(0) == 1.0
-    v(0) == 0.0
-    q(1) == 0.0                  # final condition
+    q(t0) == q0
+    v(t0) == v0
+    q(tf) == 0.0                  # final condition
     
     ẋ(t) == [v(t), -ω^2 * q(t)]
     
@@ -109,18 +222,18 @@ end
 nothing # hide
 ```
 
-### [Solution](@id example-control-free-solution-2)
+### [Direct method](@id example-control-free-direct-2)
 
-```@example main
-sol2 = solve(ocp2; display=false)
-println("Optimal pulsation: ω = ", variable(sol2))
-println("Objective value: ω² = ", objective(sol2))
+```@example main-harmonic
+direct_sol = solve(ocp; grid_size=20, display=false)
+println("Optimal pulsation: ω = ", variable(direct_sol))
+println("Objective value: ω² = ", objective(direct_sol))
 println("Expected: ω = π/2 ≈ 1.5708, ω² ≈ 2.4674")
 nothing # hide
 ```
 
-```@example main
-plot(sol2; size=(800, 400))
+```@example main-harmonic
+plot(direct_sol; size=(800, 400))
 ```
 
 The optimal pulsation should be close to $\omega = \pi/2 \approx 1.5708$, and the objective $\omega^2 \approx 2.4674$.
@@ -129,14 +242,14 @@ The optimal pulsation should be close to $\omega = \pi/2 \approx 1.5708$, and th
 
 For the harmonic oscillator, we can compare the numerical solution with the analytical one:
 
-```@example main
+```@example main-harmonic
 # analytical solution
 t_analytical = range(0, 1, 100)
 q_analytical = cos.(π * t_analytical / 2)
 v_analytical = -(π/2) * sin.(π * t_analytical / 2)
 
 # plot comparison
-plt = plot(sol2; size=(800, 600))
+plt = plot(direct_sol; size=(800, 600), label="Direct")
 plot!(plt, t_analytical, q_analytical; 
       label="q (analytical)", linestyle=:dash, linewidth=2, subplot=1)
 plot!(plt, t_analytical, v_analytical; 
@@ -145,6 +258,110 @@ plot!(plt, t_analytical, v_analytical;
 
 The numerical and analytical solutions should match very closely.
 
+### [Indirect method](@id example-control-free-indirect-2)
+
+We now solve the same problem using an indirect shooting method. For this control-free problem with a variable parameter, we use an **augmented Hamiltonian** approach. The Hamiltonian for this problem is:
+
+```math
+H(x, p, \omega) = p_1 v + p_2 (-\omega^2 q)
+```
+
+To handle the variable parameter $\omega$, we treat it as an additional state with zero dynamics. This gives us the augmented system with state $(q, v, \omega)$ and costate $(p_1, p_2, p_\omega)$, where:
+
+```math
+\begin{aligned}
+\frac{\mathrm{d}q}{\mathrm{d}t} &= \frac{\partial H}{\partial p_1} = v \\
+\frac{\mathrm{d}v}{\mathrm{d}t} &= \frac{\partial H}{\partial p_2} = -\omega^2 q \\
+\frac{\mathrm{d}\omega}{\mathrm{d}t} &= 0 \quad \text{(constant parameter)} \\
+\frac{\mathrm{d}p_1}{\mathrm{d}t} &= -\frac{\partial H}{\partial q} = \omega^2 p_2 \\
+\frac{\mathrm{d}p_2}{\mathrm{d}t} &= -\frac{\partial H}{\partial v} = -p_1 \\
+\frac{\mathrm{d}p_\omega}{\mathrm{d}t} &= -\frac{\partial H}{\partial \omega} = 2\omega q p_2
+\end{aligned}
+```
+
+For this problem with a Mayer cost $g(\omega) = \omega^2$, the transversality condition for the variable parameter is:
+
+```math
+p_\omega(t_f) - p_\omega(t_0)= -\frac{\partial g}{\partial \omega} = -2\omega
+```
+
+Assuming $p_\omega(t_0) = 0$, we have:
+
+```math
+p_\omega(t_f) = -\int_{t_0}^{t_f} \frac{\partial H}{\partial \omega}(t, x(t), p(t), \omega) \, \mathrm{d}t = -2\omega
+```
+
+We use CTFlows' `augment=true` feature to automatically compute $p_\omega(t_f)$ without manually constructing the augmented system:
+
+```@example main-harmonic
+# Create Hamiltonian flow from OCP
+f = Flow(ocp)
+nothing # hide
+```
+
+!!! note
+
+    For more details about the flow construction, see [this page](@ref manual-flow-others).
+
+The shooting function enforces the conditions:
+
+- Final condition: $q(t_f) = 0$
+- Free final velocity: $p_2(t_f) = 0$
+- Transversality condition for Mayer cost: $p_\omega(t_f) + 2\omega = 0$
+
+Using `augment=true`, the flow automatically returns $(x(t_f), p(t_f), p_\omega(t_f))$, with $p_\omega(t_0) = 0$ by construction.
+
+```@example main-harmonic
+# Shooting function: S(p0, ω)
+function shoot!(s, p0, ω)
+    x_tf, p_tf, pω_tf = f(t0, [q0, v0], p0, tf, ω; augment=true)
+    q_tf = x_tf[1]
+    pv_tf = p_tf[2]
+    s[1] = q_tf         # q(tf) = 0
+    s[2] = pv_tf        # p2(tf) = 0 (free final velocity)
+    s[3] = pω_tf + 2ω  # pω(tf) + 2ω = 0 (Mayer cost transversality)
+    return nothing
+end
+
+# Auxiliary in-place NLE function
+nle!(s, y, _) = shoot!(s, y[1:2], y[3])
+nothing # hide
+```
+
+We use the direct solution to initialize the shooting method:
+
+```@example main-harmonic
+# Extract solution from direct method for initialization
+p_direct = costate(direct_sol)
+ω_direct = variable(direct_sol)
+
+# Initial guess
+p0_guess = p_direct(t0)
+ω_guess = ω_direct
+
+# NLE problem with initial guess
+prob_indirect = NonlinearProblem(nle!, [p0_guess..., ω_guess])
+
+# Solve shooting equations
+shooting_sol = solve(prob_indirect; show_trace=Val(false))
+p0_sol, ω_sol = shooting_sol.u[1:2], shooting_sol.u[3]
+
+println("Indirect solution:")
+println("Initial costate: p0 = ", p0_sol)
+println("Parameter: ω = ", ω_sol)
+nothing # hide
+```
+
+Finally, we compute and plot the indirect solution:
+
+```@example main-harmonic
+# Compute and plot indirect solution
+indirect_sol = f((t0, tf), [q0, v0], p0_sol, ω_sol; saveat=range(t0, tf, 200))
+plot!(plt, indirect_sol; linestyle=:dash, lw=2, label="Indirect", color=2)
+```
+
+The direct and indirect solutions match closely, both finding the optimal pulsation $\omega \approx \pi/2$.
+
 !!! note "Applications"
 
     Control-free problems appear in many contexts:
diff --git a/docs/src/manual-differential-geometry.md b/docs/src/manual-differential-geometry.md
index 1875f619a..23e88329a 100644
--- a/docs/src/manual-differential-geometry.md
+++ b/docs/src/manual-differential-geometry.md
@@ -607,7 +607,7 @@ dg(3, [1, 2], [4, 5])
 For a non-autonomous Hamiltonian $H(t, x, p)$ and a function $G(t, x, p)$, the **total time derivative** along the Hamiltonian flow is:
 
 ```math
-\frac{d}{dt} G(t, x(t), p(t)) = \partial_t G + \{H, G\}.
+\frac{\mathrm{d}}{\mathrm{d}t} G(t, x(t), p(t)) = \partial_t G + \{H, G\}.
 ```
 
 This is the sum of:
diff --git a/docs/src/manual-flow-ocp.md b/docs/src/manual-flow-ocp.md
index bce982873..7df91da06 100644
--- a/docs/src/manual-flow-ocp.md
+++ b/docs/src/manual-flow-ocp.md
@@ -316,6 +316,126 @@ yf, pf = f(0, [x0, tf], [p0, 0], 1)
 
     In the [Goddard problem](https://control-toolbox.org/Tutorials.jl/stable/tutorial-goddard.html#tutorial-goddard-structure), you may find other constructions of flows, especially for singular and boundary arcs.
 
+## Augmented costate computation with `augment=true`
+
+When working with optimal control problems that have variables, it can be useful to compute the costate associated with the variable parameter. The `augment=true` keyword argument provides automatic computation of this costate without requiring manual construction of the augmented Hamiltonian system.
+
+### Mathematical background
+
+For an optimal control problem with Hamiltonian $H(t, x, p, v)$, where $x$ is the state, $p$ is the costate, and $v$ is a variable parameter, the **augmented system** treats the variable as an additional state with zero dynamics:
+
+```math
+\begin{aligned}
+\frac{\mathrm{d}x}{\mathrm{d}t} &= \frac{\partial H}{\partial p} \\
+\frac{\mathrm{d}v}{\mathrm{d}t} &= 0 \quad \text{(constant parameter)} \\
+\frac{\mathrm{d}p}{\mathrm{d}t} &= -\frac{\partial H}{\partial x} \\
+\frac{\mathrm{d}p_v}{\mathrm{d}t} &= -\frac{\partial H}{\partial v}
+\end{aligned}
+```
+
+With the initial condition $p_v(t_0) = 0$, the final costate $p_v(t_f)$ represents the accumulated sensitivity:
+
+```math
+p_v(t_f) = -\int_{t_0}^{t_f} \frac{\partial H}{\partial v}(t, x(t), p(t), v) \, \mathrm{d}t
+```
+
+### Usage
+
+Let us consider a harmonic oscillator problem where the pulsation $\omega$ is a variable parameter appearing in the dynamics:
+
+```@example main
+q0 = 1
+v0 = 0
+t0 = 0
+tf = 1
+
+ocp_aug = @def begin
+    ω ∈ R, variable              # pulsation to optimize
+    t ∈ [t0, tf], time
+    x = (q, v) ∈ R², state
+    u ∈ R, control
+    
+    q(t0) == q0
+    v(t0) == v0
+    q(tf) == 0.0
+    
+    ẋ(t) == [v(t), -ω^2 * q(t) + u(t)]
+    
+    ω^2 + 0.5∫(u(t)^2) → min
+end
+
+# Maximizing control from Pontryagin's principle
+u_aug(x, p, ω) = p[2]
+f_aug = Flow(ocp_aug, u_aug)
+nothing # hide
+```
+
+Without `augment=true`, the flow returns only the state and costate:
+
+```@example main
+ω_val = π/2
+p0_val = [1.0, 0.5]
+xf, pf = f_aug(t0, [q0, v0], p0_val, tf, ω_val)
+println("q(tf) = ", xf[1], ", v(tf) = ", xf[2])
+```
+
+With `augment=true`, the flow automatically computes and returns the costate associated with the variable `ω`:
+
+```@example main
+xf, pf, pω = f_aug(t0, [q0, v0], p0_val, tf, ω_val; augment=true)
+println("q(tf) = ", xf[1], ", v(tf) = ", xf[2], ", p_ω(tf) = ", pω)
+```
+
+The value `pω` represents the sensitivity of the Hamiltonian with respect to the pulsation parameter:
+
+```math
+p_{ω}(t_f) = -\int_{t_0}^{t_f} \frac{\partial H}{\partial \omega}(t, x(t), p(t), \omega) \, \mathrm{d}t
+```
+
+with $p_{\omega}(t_0) = 0$ by construction. This is particularly useful for computing transversality conditions in control-free problems.
+
+### Advantages
+
+The `augment=true` feature provides several benefits:
+
+- **No manual work**: No need to manually construct the augmented Hamiltonian or augmented ODEs
+- **Type-safe**: Automatic handling of scalar vs vector variables
+- **Robust**: Uses the existing, well-tested `Flow(Hamiltonian(...))` infrastructure
+- **Mathematical rigor**: Proper initial conditions and transversality handling
+
+### Error handling
+
+The `augment=true` option is only available for problems with variables:
+
+```julia
+# This will throw an error (no variable in the problem)
+ocp_no_var = @def begin
+    t ∈ [0, 1], time
+    x ∈ R, state
+    u ∈ R, control
+    x(0) == 0
+    ẋ(t) == u(t)
+    ∫(u(t)^2) → min
+end
+
+f_no_var = Flow(ocp_no_var, (x, p) -> p)
+f_no_var(0, 0, 1, 1; augment=true)  # ERROR: PreconditionError
+```
+
+Additionally, `augment=true` only works for point evaluation, not for trajectory computation:
+
+```julia
+# This works (point evaluation)
+xf, pf, pvf = f_aug(t0, x0, p0, tf, v; augment=true)
+
+# This will throw an error (trajectory call)
+sol = f_aug((t0, tf), x0, p0, v; augment=true)  # ERROR: PreconditionError
+```
+
+!!! note "Control-free problems"
+
+    The `augment=true` feature is particularly useful for control-free problems where the variable parameter appears in the dynamics. See the [control-free problems example](@ref example-control-free) for detailed applications with transversality conditions.
+
 ## Concatenation of arcs
 
 In this part, we present how to concatenate several flows. Let us consider the following problem.
diff --git a/test/suite/reexport/test_ctflows.jl b/test/suite/reexport/test_ctflows.jl
index 66ca73f31..64ad26777 100644
--- a/test/suite/reexport/test_ctflows.jl
+++ b/test/suite/reexport/test_ctflows.jl
@@ -10,6 +10,7 @@ module TestCtflows
 using Test: Test
 using OptimalControl # using is mandatory since we test exported symbols
 using CTFlows: CTFlows # needed for abstract type checks
+using OrdinaryDiffEq
 
 const VERBOSE = isdefined(Main, :TestOptions) ? Main.TestOptions.VERBOSE : true
 const SHOWTIMING = isdefined(Main, :TestOptions) ? Main.TestOptions.SHOWTIMING : true
@@ -357,6 +358,127 @@ function test_ctflows()
                 Test.@test dg(3, [1, 2], [4, 5]) ≈ 6
             end
         end
+
+        # ====================================================================
+        # FLOW FROM OCP AND AUGMENT TESTS
+        # ====================================================================
+        # These tests verify the Flow(ocp) construction for control-free problems
+        # and the augment=true feature for automatic costate computation.
+
+        Test.@testset "Flow from OCP and augment" begin
+            Test.@testset "Flow from Control-Free OCP" begin
+                # Define a simple control-free OCP (exponential growth)
+                t0 = 0
+                tf = 1
+                x0 = 1.0
+                
+                ocp = @def begin
+                    λ ∈ R, variable
+                    t ∈ [t0, tf], time
+                    x ∈ R, state
+                    x(t0) == x0
+                    ẋ(t) == λ * x(t)
+                    ∫(x(t)^2) → min
+                end
+                
+                # Test: Flow(ocp) works for control-free problems
+                f = Flow(ocp)
+                
+                # Test: basic call returns 2 values (state, costate)
+                λ_val = 0.5
+                p0 = 1.0
+                Test.@test applicable(f, t0, x0, p0, tf, λ_val)
+                xf, pf = f(t0, x0, p0, tf, λ_val)
+                Test.@test xf isa Real
+                Test.@test pf isa Real
+            end
+
+            Test.@testset "Flow with augment=true" begin
+                # Same OCP as above
+                t0 = 0
+                tf = 1
+                x0 = 1.0
+                
+                ocp = @def begin
+                    λ ∈ R, variable
+                    t ∈ [t0, tf], time
+                    x ∈ R, state
+                    x(t0) == x0
+                    ẋ(t) == λ * x(t)
+                    ∫(x(t)^2) → min
+                end
+                
+                f = Flow(ocp)
+                λ_val = 0.5
+                p0 = 1.0
+                
+                # Test: augment=true returns 3 values (state, costate, variable costate)
+                xf, pf, pλ = f(t0, x0, p0, tf, λ_val; augment=true)
+                Test.@test xf isa Real
+                Test.@test pf isa Real
+                Test.@test pλ isa Real  # The new one: costate of λ
+            end
+
+            Test.@testset "Manual vs Automatic Hamiltonian" begin
+                # Define OCP
+                t0 = 0
+                tf = 1
+                x0 = 1.0
+                λ_val = 0.5
+                p0 = 1.0
+                
+                ocp = @def begin
+                    λ ∈ R, variable
+                    t ∈ [t0, tf], time
+                    x ∈ R, state
+                    x(t0) == x0
+                    ẋ(t) == λ * x(t)
+                    ∫(x(t)^2) → min
+                end
+                
+                # Manual Hamiltonian construction
+                H(x, p, λ) = p * λ * x - x^2
+                function H_aug(x_, p_)
+                    x, λ = x_
+                    p, _ = p_
+                    return H(x, p, λ)
+                end
+                f_manual = Flow(OptimalControl.Hamiltonian(H_aug))
+                
+                # Automatic Flow from OCP
+                f_auto = Flow(ocp)
+                
+                # Test: both give similar results
+                xf_manual, pf_manual = f_manual(t0, [x0, λ_val], [p0, 0.0], tf)
+                xf_auto, pf_auto = f_auto(t0, x0, p0, tf, λ_val)
+                
+                Test.@test xf_manual[1] ≈ xf_auto rtol=1e-6
+                Test.@test pf_manual[1] ≈ pf_auto rtol=1e-6
+            end
+
+            Test.@testset "Analytical Solution Check" begin
+                # For λ=0, x(t) = x0 (constant)
+                t0 = 0
+                tf = 1
+                x0 = 1.0
+                
+                ocp = @def begin
+                    λ ∈ R, variable
+                    t ∈ [t0, tf], time
+                    x ∈ R, state
+                    x(t0) == x0
+                    ẋ(t) == λ * x(t)
+                    ∫(x(t)^2) → min
+                end
+                
+                f = Flow(ocp)
+                λ_zero = 0.0
+                p0 = 1.0
+                
+                xf, pf = f(t0, x0, p0, tf, λ_zero)
+                Test.@test xf ≈ x0 rtol=1e-10  # x remains constant
+            end
+        end
     end
 end