diff --git a/docs/src/assets/Manifest.toml b/docs/src/assets/Manifest.toml
index 2140cea3..0f77cd71 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 = "d89a9defa5901dd1d745042c8e1ea5af87977104"
+git-tree-sha1 = "536c34374ab27f3d4d6a51d9e92ed89f146ca989"
 uuid = "1c39547c-7794-42f7-af83-d98194f657c2"
-version = "0.8.19-beta"
+version = "0.8.20"
 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 = "c878ed38da4040f618a459c30cec8e53286ddc26"
+git-tree-sha1 = "3e9df6b6cb96ccb05051ded9a5d4f76f649f6d0c"
 uuid = "34c4fa32-2049-4079-8329-de33c2a22e2d"
-version = "0.9.11"
+version = "0.9.12-beta"
 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 = "311bad1284ba1ccb23a6efe1d0cf531bdb7c46b9"
+git-tree-sha1 = "1a7445af1c937c291371e0f7c0a51a80f5b5d0eb"
 uuid = "32681960-a1b1-40db-9bff-a1ca817385d1"
-version = "0.8.10-beta"
+version = "0.8.12-beta"
 
 [[deps.CTSolvers]]
 deps = ["ADNLPModels", "CTBase", "CTModels", "CommonSolve", "DocStringExtensions", "ExaModels", "KernelAbstractions", "NLPModels", "SolverCore"]
diff --git a/docs/src/example-control-free.md b/docs/src/example-control-free.md
index dc3b3852..d4023b3a 100644
--- a/docs/src/example-control-free.md
+++ b/docs/src/example-control-free.md
@@ -7,10 +7,6 @@ Control-free problems are optimal control problems without a control variable. T
 
 This page demonstrates two simple examples with known analytical solutions.
 
-!!! compat "Upcoming feature"
-
-    The syntax shown here uses a workaround with a dummy control variable (`u ∈ R, control`) and a small penalty term (`1e-5*u(t)^2`) in the cost to force `u ≈ 0`, because the parser does not yet support omitting the control declaration. These workaround lines are marked with `# TO REMOVE WHEN POSSIBLE` and will be removed once native control-free syntax is implemented.
-
 First, we import the necessary packages:
 
 ```@example main
@@ -46,13 +42,11 @@ ocp1 = @def begin
     p ∈ R, variable              # growth rate to estimate
     t ∈ [0, 10], time
     x ∈ R, state
-    u ∈ R, control               # TO REMOVE WHEN POSSIBLE
     
     x(0) == 2.0
-    
     ẋ(t) == p * x(t)
     
-    ∫((x(t) - data(t))^2 + 1e-5*u(t)^2) → min  # fit to observed data # TO REMOVE u(t) when possible
+    ∫((x(t) - data(t))^2) → min  # fit to observed data
 end
 nothing # hide
 ```
@@ -103,7 +97,6 @@ ocp2 = @def begin
     ω ∈ R, variable              # pulsation to optimize
     t ∈ [0, 1], time
     x = (q, v) ∈ R², state
-    u ∈ R, control               # TO REMOVE WHEN POSSIBLE
     
     q(0) == 1.0
     v(0) == 0.0
@@ -111,7 +104,7 @@ ocp2 = @def begin
     
     ẋ(t) == [v(t), -ω^2 * q(t)]
     
-    ω^2 + 1e-5*∫(u(t)^2) → min   # minimize pulsation # TO REMOVE u(t) when possible
+    ω^2 → min   # minimize pulsation
 end
 nothing # hide
 ```
diff --git a/docs/src/example-double-integrator-energy.md b/docs/src/example-double-integrator-energy.md
index 7d908e6d..54c88cbd 100644
--- a/docs/src/example-double-integrator-energy.md
+++ b/docs/src/example-double-integrator-energy.md
@@ -135,8 +135,10 @@ We are now ready to solve the shooting equations.
 # auxiliary in-place NLE function
 nle!(s, p0, _) = s[:] = S(p0)
 
-# initial guess for the Newton solver
-p0_guess = [1, 1]
+# initial guess for the Newton solver from the direct solution
+t = time_grid(direct_sol) # the time grid as a vector
+p = costate(direct_sol)   # the costate as a function of time
+p0_guess = p(t0)          # initial costate
 
 # NLE problem with initial guess
 prob = NonlinearProblem(nle!, p0_guess)
diff --git a/docs/src/example-double-integrator-time.md b/docs/src/example-double-integrator-time.md
index 06c6096c..5da17658 100644
--- a/docs/src/example-double-integrator-time.md
+++ b/docs/src/example-double-integrator-time.md
@@ -90,14 +90,14 @@ nothing # hide
 Let us solve it with a direct method (we set the number of time steps to 200):
 
 ```@example main
-sol = solve(ocp; grid_size=200)
+direct_sol = solve(ocp; grid_size=200)
 nothing # hide
 ```
 
 and plot the solution:
 
 ```@example main
-plt = plot(sol; label="Direct", size=(800, 600))
+plt = plot(direct_sol; label="Direct", size=(800, 600))
 ```
 
 !!! note "Nota bene"
@@ -161,15 +161,24 @@ We are now ready to solve the shooting equations:
 # in-place shooting function
 nle!(s, ξ, λ) = shoot!(s, ξ[1:2], ξ[3], ξ[4]) 
 
-# initial guess: costate and final time
-ξ_guess = [0.1, 0.1, 0.5, 1]
+# initial guess for the Newton solver from direct method
+t = time_grid(direct_sol) # the time grid as a vector
+p = costate(direct_sol)   # the costate as a function of time
+p0_guess = p(t0)          # initial costate
+tf_guess = variable(direct_sol) # final time
+
+# find switching time t1 where p2(t) changes sign
+p2_values = [p(ti)[2] for ti in t]
+t1_guess = t[findfirst(i -> i > 1 && p2_values[i] * p2_values[i-1] < 0, 2:length(p2_values))]
+
+ξ_guess = [p0_guess[1], p0_guess[2], t1_guess, tf_guess]
 
 # NLE problem
 prob = NonlinearProblem(nle!, ξ_guess)
 
 # resolution of the shooting equations
-sol = solve(prob; show_trace=Val(true))
-p0, t1, tf = sol.u[1:2], sol.u[3], sol.u[4]
+shoot_sol = solve(prob; show_trace=Val(true))
+p0, t1, tf = shoot_sol.u[1:2], shoot_sol.u[3], shoot_sol.u[4]
 
 # print the solution
 println("\np0 = ", p0, "\nt1 = ", t1, "\ntf = ", tf)
@@ -182,10 +191,10 @@ Finally, we reconstruct and plot the solution obtained by the indirect method:
 φ = f_max * (t1, f_min)
 
 # compute the solution: state, costate, control...
-flow_sol = φ((t0, tf), x0, p0; saveat=range(t0, tf, 200))
+indirect_sol = φ((t0, tf), x0, p0; saveat=range(t0, tf, 200))
 
 # plot the solution on the previous plot
-plot!(plt, flow_sol; label="Indirect", color=2, linestyle=:dash)
+plot!(plt, indirect_sol; label="Indirect", color=2, linestyle=:dash)
 ```
 
 !!! note
diff --git a/docs/src/example-singular-control.md b/docs/src/example-singular-control.md
index fb38d094..fdd32a47 100644
--- a/docs/src/example-singular-control.md
+++ b/docs/src/example-singular-control.md
@@ -77,7 +77,7 @@ Let's plot the solution:
 
 ```@example main
 opt = (state_bounds_style=:none, control_bounds_style=:none)
-plt = plot(direct_sol; label="direct", size=(800, 800), opt...)
+plt = plot(direct_sol; label="Direct", size=(800, 800), opt...)
 ```
 
 ## Singular control by hand
@@ -190,7 +190,7 @@ u_s(\theta) = \sin^2\theta.
 Let's overlay this on the numerical solution:
 
 ```@example main
-T = time_grid(direct_sol, :control)
+T = time_grid(direct_sol)
 θ(t) = state(direct_sol)(t)[3]
 us(t) = sin(θ(t))^2
 plot!(plt, T, us; subplot=7, line=:dash, lw=2, label="us (hand)")
@@ -340,7 +340,7 @@ nothing # hide
 Plot the indirect solution alongside the direct solution:
 
 ```@example main
-plot!(plt, indirect_sol; label="indirect", color=2, linestyle=:dash, opt...)
+plot!(plt, indirect_sol; label="Indirect", color=2, linestyle=:dash, opt...)
 ```
 
 The indirect and direct solutions match very well, confirming that our singular control computation is correct.
diff --git a/docs/src/manual-solve.md b/docs/src/manual-solve.md
index 94500297..a0e1aeb7 100644
--- a/docs/src/manual-solve.md
+++ b/docs/src/manual-solve.md
@@ -230,6 +230,22 @@ Each strategy declares its available options. You can inspect them using `descri
 
 ### Discretizer options
 
+The collocation discretizer supports multiple integration schemes:
+
+- `:trapeze` - Trapezoidal rule (second-order accurate)
+- `:midpoint` - Midpoint rule (second-order accurate)  
+- `:euler` or `:euler_explicit` or `:euler_forward` - Explicit Euler method (first-order accurate)
+- `:euler_implicit` or `:euler_backward` - Implicit Euler method (first-order accurate, more stable for stiff problems)
+
+!!! note "Additional schemes with ADNLP modeler"
+
+    When using the `:adnlp` modeler, two additional high-order collocation schemes are available:
+
+    - `:gauss_legendre_2` - 2-point Gauss-Legendre collocation (fourth-order accurate)
+    - `:gauss_legendre_3` - 3-point Gauss-Legendre collocation (sixth-order accurate)
+
+    These schemes provide higher accuracy but require more computational effort.
+
 ```@example main
 describe(:collocation)
 ```
diff --git a/test/problems/control_free.jl b/test/problems/control_free.jl
index 19ce173e..5f5ebe30 100644
--- a/test/problems/control_free.jl
+++ b/test/problems/control_free.jl
@@ -29,13 +29,12 @@ function ExponentialGrowth()
         p ∈ R, variable              # growth rate to estimate
         t ∈ [0, 10], time
         x ∈ R, state
-        u ∈ R, control               # TO REMOVE WHEN POSSIBLE
 
         x(0) == 2.0
 
         ẋ(t) == p * x(t)
 
-        ∫((x(t) - data(t))^2 + 1e-5*u(t)^2) → min  # fit to observed data # TO REMOVE u(t) when possible
+        ∫((x(t) - data(t))^2) → min  # fit to observed data
     end
 
     return (
@@ -73,7 +72,6 @@ function HarmonicOscillator()
         ω ∈ R, variable              # pulsation to optimize
         t ∈ [0, 1], time
         x = (q, v) ∈ R², state
-        u ∈ R, control               # TO REMOVE WHEN POSSIBLE
 
         q(0) == 1.0
         v(0) == 0.0
@@ -81,7 +79,7 @@ function HarmonicOscillator()
 
         ẋ(t) == [v(t), -ω^2 * q(t)]
 
-        ω^2 + 1e-5*∫(u(t)^2) → min   # minimize pulsation # TO REMOVE u(t) when possible
+        ω^2 → min   # minimize pulsation
     end
 
     return (