diff --git a/docs/make.jl b/docs/make.jl
index 5609daab..29ed025b 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -207,6 +207,8 @@ with_api_reference(src_dir, ext_dir) do api_pages
             "Basic Examples" => [
                 "Energy minimisation" => "example-double-integrator-energy.md",
                 "Time mininimisation" => "example-double-integrator-time.md",
+                "Control-free problems" => "example-control-free.md",
+                "Singular control" => "example-singular-control.md",
             ],
             "Manual" => [
                 "Define a problem" => "manual-abstract.md",
@@ -221,6 +223,7 @@ with_api_reference(src_dir, ext_dir) do api_pages
                 ],
                 "Solution characteristics" => "manual-solution.md",
                 "Plot a solution" => "manual-plot.md",
+                "Differential geometry tools" => "manual-differential-geometry.md",
                 "Compute flows" => [
                     "From optimal control problems" => "manual-flow-ocp.md",
                     "From Hamiltonians and others" => "manual-flow-others.md",
diff --git a/docs/src/assets/Manifest.toml b/docs/src/assets/Manifest.toml
index 64f46c30..2140cea3 100644
--- a/docs/src/assets/Manifest.toml
+++ b/docs/src/assets/Manifest.toml
@@ -1,6 +1,6 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.12.1"
+julia_version = "1.12.5"
 manifest_format = "2.0"
 project_hash = "d470f84ad521e691a87a863486162d9840a77916"
 
@@ -56,9 +56,9 @@ version = "0.4.5"
 
 [[deps.Accessors]]
 deps = ["CompositionsBase", "ConstructionBase", "Dates", "InverseFunctions", "MacroTools"]
-git-tree-sha1 = "856ecd7cebb68e5fc87abecd2326ad59f0f911f3"
+git-tree-sha1 = "2eeb2c9bef11013efc6f8f97f32ee59b146b09fb"
 uuid = "7d9f7c33-5ae7-4f3b-8dc6-eff91059b697"
-version = "0.1.43"
+version = "0.1.44"
 
     [deps.Accessors.extensions]
     AxisKeysExt = "AxisKeys"
@@ -141,9 +141,9 @@ version = "1.11.0"
 
 [[deps.Atomix]]
 deps = ["UnsafeAtomics"]
-git-tree-sha1 = "29bb0eb6f578a587a49da16564705968667f5fa8"
+git-tree-sha1 = "b8651b2eb5796a386b0398a20b519a6a6150f75c"
 uuid = "a9b6321e-bd34-4604-b9c9-b65b8de01458"
-version = "1.1.2"
+version = "1.1.3"
 
     [deps.Atomix.extensions]
     AtomixCUDAExt = "CUDA"
@@ -180,9 +180,9 @@ version = "0.1.6"
 
 [[deps.BracketingNonlinearSolve]]
 deps = ["CommonSolve", "ConcreteStructs", "NonlinearSolveBase", "PrecompileTools", "Reexport", "SciMLBase"]
-git-tree-sha1 = "4999dff8efd76814f6662519b985aeda975a1924"
+git-tree-sha1 = "d9b66401c1fa982c7ca984d0566af5a9b3551420"
 uuid = "70df07ce-3d50-431d-a3e7-ca6ddb60ac1e"
-version = "1.11.0"
+version = "1.12.0"
 weakdeps = ["ChainRulesCore", "ForwardDiff"]
 
     [deps.BracketingNonlinearSolve.extensions]
@@ -208,9 +208,9 @@ version = "0.2.7"
 
 [[deps.CTBase]]
 deps = ["DocStringExtensions"]
-git-tree-sha1 = "fac598a2d962f564824d37494a307652941c723d"
+git-tree-sha1 = "2f689cdb1c28d452f143f2c76897cd5273e89cf7"
 uuid = "54762871-cc72-4466-b8e8-f6c8b58076cd"
-version = "0.18.6-beta"
+version = "0.18.7"
 
     [deps.CTBase.extensions]
     CoveragePostprocessing = ["Coverage"]
@@ -225,16 +225,16 @@ version = "0.18.6-beta"
     Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [[deps.CTDirect]]
-deps = ["ADNLPModels", "CTModels", "CTSolvers", "DocStringExtensions", "ExaModels", "SolverCore", "SparseArrays"]
-git-tree-sha1 = "9d9c1c5c1d17cc02692378ae49d03e90b0e469e8"
+deps = ["ADNLPModels", "CTModels", "CTSolvers", "DocStringExtensions", "ExaModels", "SolverCore", "SparseArrays", "SparseConnectivityTracer"]
+git-tree-sha1 = "db1699ffab2f31aeda5dbb61dfa9015773cb70fc"
 uuid = "790bbbee-bee9-49ee-8912-a9de031322d5"
-version = "1.0.5"
+version = "1.0.6"
 
 [[deps.CTFlows]]
 deps = ["CTBase", "CTModels", "DocStringExtensions", "ForwardDiff", "LinearAlgebra", "MLStyle", "MacroTools"]
-git-tree-sha1 = "71ec5ebc9464b1d79e64ea2e53675fa0506e20c6"
+git-tree-sha1 = "d89a9defa5901dd1d745042c8e1ea5af87977104"
 uuid = "1c39547c-7794-42f7-af83-d98194f657c2"
-version = "0.8.16-beta"
+version = "0.8.19-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 = "c90116a8e74243ae258329325f0b7545cc5385f5"
+git-tree-sha1 = "c878ed38da4040f618a459c30cec8e53286ddc26"
 uuid = "34c4fa32-2049-4079-8329-de33c2a22e2d"
-version = "0.9.10-beta"
+version = "0.9.11"
 weakdeps = ["JLD2", "JSON3", "Plots"]
 
     [deps.CTModels.extensions]
@@ -260,9 +260,9 @@ version = "0.8.10-beta"
 
 [[deps.CTSolvers]]
 deps = ["ADNLPModels", "CTBase", "CTModels", "CommonSolve", "DocStringExtensions", "ExaModels", "KernelAbstractions", "NLPModels", "SolverCore"]
-git-tree-sha1 = "21092a83de84623f3d2e55b301ffa54d8cb10e3b"
+git-tree-sha1 = "46f28bf27a6158333527417d97787b8125b1a07c"
 uuid = "d3e8d392-8e4b-4d9b-8e92-d7d4e3650ef6"
-version = "0.4.10-beta"
+version = "0.4.12"
 
     [deps.CTSolvers.extensions]
     CTSolversCUDA = "CUDA"
@@ -341,16 +341,16 @@ uuid = "4889d778-9329-5762-9fec-0578a5d30366"
 version = "0.7.1+0"
 
 [[deps.Cairo_jll]]
-deps = ["Artifacts", "Bzip2_jll", "CompilerSupportLibraries_jll", "Fontconfig_jll", "FreeType2_jll", "Glib_jll", "JLLWrappers", "LZO_jll", "Libdl", "Pixman_jll", "Xorg_libXext_jll", "Xorg_libXrender_jll", "Zlib_jll", "libpng_jll"]
-git-tree-sha1 = "a21c5464519504e41e0cbc91f0188e8ca23d7440"
+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"
 uuid = "83423d85-b0ee-5818-9007-b63ccbeb887a"
-version = "1.18.5+1"
+version = "1.18.6+0"
 
 [[deps.ChainRulesCore]]
 deps = ["Compat", "LinearAlgebra"]
-git-tree-sha1 = "e4c6a16e77171a5f5e25e9646617ab1c276c5607"
+git-tree-sha1 = "12177ad6b3cad7fd50c8b3825ce24a99ad61c18f"
 uuid = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
-version = "1.26.0"
+version = "1.26.1"
 weakdeps = ["SparseArrays"]
 
     [deps.ChainRulesCore.extensions]
@@ -512,9 +512,9 @@ version = "1.8.1"
 
 [[deps.DataStructures]]
 deps = ["OrderedCollections"]
-git-tree-sha1 = "e357641bb3e0638d353c4b29ea0e40ea644066a6"
+git-tree-sha1 = "e86f4a2805f7f19bec5129bc9150c38208e5dc23"
 uuid = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
-version = "0.19.3"
+version = "0.19.4"
 
 [[deps.DataValueInterfaces]]
 git-tree-sha1 = "bfc1187b79289637fa0ef6d4436ebdfe6905cbd6"
@@ -540,13 +540,14 @@ version = "1.9.1"
 
 [[deps.DiffEqBase]]
 deps = ["ArrayInterface", "BracketingNonlinearSolve", "ConcreteStructs", "DocStringExtensions", "FastBroadcast", "FastClosures", "FastPower", "FunctionWrappers", "FunctionWrappersWrappers", "LinearAlgebra", "Logging", "Markdown", "MuladdMacro", "PrecompileTools", "Printf", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators", "SciMLStructures", "Setfield", "Static", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"]
-git-tree-sha1 = "1719cd1b0a12e01775dc6db1577dd6ace1798fee"
+git-tree-sha1 = "803d463050fae6121f5d15fd449e9538cb993ad9"
 uuid = "2b5f629d-d688-5b77-993f-72d75c75574e"
-version = "6.210.1"
+version = "6.214.1"
 
     [deps.DiffEqBase.extensions]
     DiffEqBaseCUDAExt = "CUDA"
     DiffEqBaseChainRulesCoreExt = "ChainRulesCore"
+    DiffEqBaseDynamicQuantitiesExt = "DynamicQuantities"
     DiffEqBaseEnzymeExt = ["ChainRulesCore", "Enzyme"]
     DiffEqBaseFlexUnitsExt = "FlexUnits"
     DiffEqBaseForwardDiffExt = ["ForwardDiff"]
@@ -565,6 +566,7 @@ version = "6.210.1"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
     ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
     Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+    DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
     Enzyme = "7da242da-08ed-463a-9acd-ee780be4f1d9"
     FlexUnits = "76e01b6b-c995-4ce6-8559-91e72a3d4e95"
     ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
@@ -678,7 +680,7 @@ version = "0.2.0"
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
 uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
-version = "1.6.0"
+version = "1.7.0"
 
 [[deps.EnumX]]
 git-tree-sha1 = "c49898e8438c828577f04b92fc9368c388ac783c"
@@ -686,9 +688,9 @@ uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 version = "1.0.7"
 
 [[deps.EnzymeCore]]
-git-tree-sha1 = "990991b8aa76d17693a98e3a915ac7aa49f08d1a"
+git-tree-sha1 = "24bbb6fc8fb87eb71c1f8d00184a60fc22c63903"
 uuid = "f151be2c-9106-41f4-ab19-57ee4f262869"
-version = "0.8.18"
+version = "0.8.19"
 weakdeps = ["Adapt", "ChainRulesCore"]
 
     [deps.EnzymeCore.extensions]
@@ -779,10 +781,15 @@ uuid = "b22a6f82-2f65-5046-a5b2-351ab43fb4e5"
 version = "8.0.1+1"
 
 [[deps.FastBroadcast]]
-deps = ["ArrayInterface", "LinearAlgebra", "Polyester", "Static", "StaticArrayInterface", "StrideArraysCore"]
-git-tree-sha1 = "ab1b34570bcdf272899062e1a56285a53ecaae08"
+deps = ["ArrayInterface", "LinearAlgebra"]
+git-tree-sha1 = "e3e64918b1604ba8b1734c4a27febdfe5d09e235"
 uuid = "7034ab61-46d4-4ed7-9d0f-46aef9175898"
-version = "0.3.5"
+version = "1.3.1"
+weakdeps = ["Polyester", "Static"]
+
+    [deps.FastBroadcast.extensions]
+    FastBroadcastPolyesterExt = "Polyester"
+    FastBroadcastStaticExt = "Static"
 
 [[deps.FastClosures]]
 git-tree-sha1 = "acebe244d53ee1b461970f8910c235b259e772ef"
@@ -887,9 +894,9 @@ version = "1.3.7"
 
 [[deps.ForwardDiff]]
 deps = ["CommonSubexpressions", "DiffResults", "DiffRules", "LinearAlgebra", "LogExpFunctions", "NaNMath", "Preferences", "Printf", "Random", "SpecialFunctions"]
-git-tree-sha1 = "eef4c86803f47dcb61e9b8790ecaa96956fdd8ae"
+git-tree-sha1 = "cddeab6487248a39dae1a960fff0ac17b2a28888"
 uuid = "f6369f11-7733-5829-9624-2563aa707210"
-version = "1.3.2"
+version = "1.3.3"
 weakdeps = ["StaticArrays"]
 
     [deps.ForwardDiff.extensions]
@@ -897,9 +904,9 @@ weakdeps = ["StaticArrays"]
 
 [[deps.FreeType2_jll]]
 deps = ["Artifacts", "Bzip2_jll", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "2c5512e11c791d1baed2049c5652441b28fc6a31"
+git-tree-sha1 = "70329abc09b886fd2c5d94ad2d9527639c421e3e"
 uuid = "d7e528f0-a631-5988-bf34-fe36492bcfd7"
-version = "2.13.4+0"
+version = "2.14.3+1"
 
 [[deps.FriBidi_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -913,10 +920,10 @@ uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
 version = "1.1.3"
 
 [[deps.FunctionWrappersWrappers]]
-deps = ["FunctionWrappers"]
-git-tree-sha1 = "b104d487b34566608f8b4e1c39fb0b10aa279ff8"
+deps = ["FunctionWrappers", "PrecompileTools", "TruncatedStacktraces"]
+git-tree-sha1 = "6874da243fb93e34201d7d4587ffa0e920682f64"
 uuid = "77dc65aa-8811-40c2-897b-53d922fa7daf"
-version = "0.1.3"
+version = "1.0.0"
 
 [[deps.Future]]
 deps = ["Random"]
@@ -947,15 +954,15 @@ version = "0.2.0"
 
 [[deps.GPUCompiler]]
 deps = ["ExprTools", "InteractiveUtils", "LLVM", "Libdl", "Logging", "PrecompileTools", "Preferences", "Scratch", "Serialization", "TOML", "Tracy", "UUIDs"]
-git-tree-sha1 = "966946d226e8b676ca6409454718accb18c34c54"
+git-tree-sha1 = "fedfe5e7db7035271c3f58359007f971da1dde87"
 uuid = "61eb1bfa-7361-4325-ad38-22787b887f55"
-version = "1.8.2"
+version = "1.9.1"
 
 [[deps.GPUToolbox]]
 deps = ["LLVM"]
-git-tree-sha1 = "9e9186b09a13b7f094f87d1a9bb266d8780e1b1c"
+git-tree-sha1 = "a589b6c1a0eff953571f5d8b0474f5020831114d"
 uuid = "096a3bc2-3ced-46d0-87f4-dd12716f4bfc"
-version = "1.0.0"
+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"]
@@ -1059,9 +1066,9 @@ version = "1.13.1+0"
 
 [[deps.Hwloc_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "XML2_jll", "Xorg_libpciaccess_jll"]
-git-tree-sha1 = "157e2e5838984449e44af851a52fe374d56b9ada"
+git-tree-sha1 = "baaaebd42ed9ee1bd9173cfd56910e55a8622ee1"
 uuid = "e33a78d0-f292-5ffc-b300-72abe9b543c8"
-version = "2.13.0+0"
+version = "2.13.0+1"
 
 [[deps.IOCapture]]
 deps = ["Logging", "Random"]
@@ -1143,9 +1150,9 @@ version = "1.0.0"
 
 [[deps.JLD2]]
 deps = ["ChunkCodecLibZlib", "ChunkCodecLibZstd", "FileIO", "MacroTools", "Mmap", "OrderedCollections", "PrecompileTools", "ScopedValues"]
-git-tree-sha1 = "8f8ff711442d1f4cfc0d86133e7ee03d62ec9b98"
+git-tree-sha1 = "941f87a0ae1b14d1ac2fa57245425b23a9d7a516"
 uuid = "033835bb-8acc-5ee8-8aae-3f567f8a3819"
-version = "0.6.3"
+version = "0.6.4"
 weakdeps = ["UnPack"]
 
     [deps.JLD2.extensions]
@@ -1165,9 +1172,9 @@ version = "1.7.1"
 
 [[deps.JSON]]
 deps = ["Dates", "Logging", "Parsers", "PrecompileTools", "StructUtils", "UUIDs", "Unicode"]
-git-tree-sha1 = "b3ad4a0255688dcb895a52fafbaae3023b588a90"
+git-tree-sha1 = "67c6f1f085cb2671c93fe34244c9cccde30f7a26"
 uuid = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
-version = "1.4.0"
+version = "1.5.0"
 
     [deps.JSON.extensions]
     JSONArrowExt = ["ArrowTypes"]
@@ -1230,9 +1237,9 @@ version = "15.1.0"
 
 [[deps.KernelAbstractions]]
 deps = ["Adapt", "Atomix", "InteractiveUtils", "MacroTools", "PrecompileTools", "Requires", "StaticArrays", "UUIDs"]
-git-tree-sha1 = "fb14a863240d62fbf5922bf9f8803d7df6c62dc8"
+git-tree-sha1 = "f2e76d3ced51a2a9e185abc0b97494c7273f649f"
 uuid = "63c18a36-062a-441e-b654-da1e3ab1ce7c"
-version = "0.9.40"
+version = "0.9.41"
 weakdeps = ["EnzymeCore", "LinearAlgebra", "SparseArrays"]
 
     [deps.KernelAbstractions.extensions]
@@ -1254,9 +1261,9 @@ version = "3.100.3+0"
 
 [[deps.LDLFactorizations]]
 deps = ["AMD", "LinearAlgebra", "SparseArrays", "Test"]
-git-tree-sha1 = "70f582b446a1c3ad82cf87e62b878668beef9d13"
+git-tree-sha1 = "d75c5cb8d6ac9c359ae9eb8e87e446ba9f221dd4"
 uuid = "40e66cde-538c-5869-a4ad-c39174c6795b"
-version = "0.10.1"
+version = "0.10.2"
 
 [[deps.LERC_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
@@ -1291,12 +1298,6 @@ git-tree-sha1 = "eb62a3deb62fc6d8822c0c4bef73e4412419c5d8"
 uuid = "1d63c593-3942-5779-bab2-d838dc0a180e"
 version = "18.1.8+0"
 
-[[deps.LZO_jll]]
-deps = ["Artifacts", "JLLWrappers", "Libdl"]
-git-tree-sha1 = "1c602b1127f4751facb671441ca72715cc95938a"
-uuid = "dd4b983a-f0e5-5f8d-a1b7-129d4a5fb1ac"
-version = "2.10.3+0"
-
 [[deps.LaTeXStrings]]
 git-tree-sha1 = "dda21b8cbd6a6c40d9d02a73230f9d70fed6918c"
 uuid = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
@@ -1344,7 +1345,7 @@ version = "0.6.4"
 [[deps.LibCURL_jll]]
 deps = ["Artifacts", "LibSSH2_jll", "Libdl", "OpenSSL_jll", "Zlib_jll", "nghttp2_jll"]
 uuid = "deac9b47-8bc7-5906-a0fe-35ac56dc84c0"
-version = "8.11.1+1"
+version = "8.15.0+0"
 
 [[deps.LibGit2]]
 deps = ["LibGit2_jll", "NetworkOptions", "Printf", "SHA"]
@@ -1451,10 +1452,10 @@ version = "2.13.0"
     Metal = "dde4c033-4e86-420c-a63e-0dd931031962"
 
 [[deps.LinearSolve]]
-deps = ["ArrayInterface", "ChainRulesCore", "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 = "57a7bea58da7de6906f2621294ea35656cb40c5f"
+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 = "1ddad5f2b0717f71f1588b3519e2f7d80fc2ce65"
 uuid = "7ed4a6bd-45f5-4d41-b270-4a48e9bafcae"
-version = "3.66.0"
+version = "3.69.0"
 
     [deps.LinearSolve.extensions]
     LinearSolveAMDGPUExt = "AMDGPU"
@@ -1465,6 +1466,7 @@ version = "3.66.0"
     LinearSolveCUDAExt = "CUDA"
     LinearSolveCUDSSExt = "CUDSS"
     LinearSolveCUSOLVERRFExt = ["CUSOLVERRF", "SparseArrays"]
+    LinearSolveChainRulesCoreExt = "ChainRulesCore"
     LinearSolveCliqueTreesExt = ["CliqueTrees", "SparseArrays"]
     LinearSolveEnzymeExt = ["EnzymeCore", "SparseArrays"]
     LinearSolveFastAlmostBandedMatricesExt = "FastAlmostBandedMatrices"
@@ -1492,6 +1494,7 @@ version = "3.66.0"
     CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
     CUDSS = "45b445bb-4962-46a0-9369-b4df9d0f772e"
     CUSOLVERRF = "a8cc9031-bad2-4722-94f5-40deabb4245c"
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
     CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
     EnzymeCore = "f151be2c-9106-41f4-ab19-57ee4f262869"
     FastAlmostBandedMatrices = "9d29842c-ecb8-4973-b1e9-a27b1157504e"
@@ -1569,9 +1572,9 @@ version = "0.5.16"
 
 [[deps.MadNCL]]
 deps = ["Atomix", "KernelAbstractions", "LinearAlgebra", "MadNLP", "NLPModels", "Printf", "Random", "SparseArrays"]
-git-tree-sha1 = "fd3e96ec2e760b48084345414b7943a7b607aa0c"
+git-tree-sha1 = "8cdb50494fa73f9af44aabadfac51d39413a707e"
 uuid = "434a0bcb-5a7c-42b2-a9d3-9e3f760e7af0"
-version = "0.2.0"
+version = "0.2.1"
 weakdeps = ["MadNLPGPU"]
 
     [deps.MadNCL.extensions]
@@ -1678,7 +1681,7 @@ version = "0.3.7"
 
 [[deps.MozillaCACerts_jll]]
 uuid = "14a3606d-f60d-562e-9121-12d972cd8159"
-version = "2025.5.20"
+version = "2025.11.4"
 
 [[deps.MuladdMacro]]
 git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab"
@@ -1699,9 +1702,9 @@ version = "0.11.2"
 
 [[deps.NLPModelsKnitro]]
 deps = ["KNITRO", "NLPModels", "NLPModelsModifiers", "SolverCore"]
-git-tree-sha1 = "f35898a07a02c4f46bc08eca0b19116e894ef464"
+git-tree-sha1 = "90e9e12bf859aa3b9453a852015a552592b8ea95"
 uuid = "bec4dd0d-7755-52d5-9a02-22f0ffc7efcb"
-version = "0.10.1"
+version = "0.10.2"
 
 [[deps.NLPModelsModifiers]]
 deps = ["FastClosures", "LinearAlgebra", "LinearOperators", "NLPModels", "Printf", "SparseArrays"]
@@ -1742,10 +1745,10 @@ 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", "SimpleNonlinearSolve", "StaticArraysCore", "SymbolicIndexingInterface"]
-git-tree-sha1 = "d27bcf0cebf8786edcc2eaa4455c959e680334e7"
+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 = "373b688150f160b68083733657246532389e4280"
 uuid = "8913a72c-1f9b-4ce2-8d82-65094dcecaec"
-version = "4.16.0"
+version = "4.17.0"
 
     [deps.NonlinearSolve.extensions]
     NonlinearSolveFastLevenbergMarquardtExt = "FastLevenbergMarquardt"
@@ -1775,10 +1778,10 @@ version = "4.16.0"
     Sundials = "c3572dad-4567-51f8-b174-8c6c989267f4"
 
 [[deps.NonlinearSolveBase]]
-deps = ["ADTypes", "Adapt", "ArrayInterface", "CommonSolve", "Compat", "ConcreteStructs", "DifferentiationInterface", "EnzymeCore", "FastClosures", "LinearAlgebra", "LogExpFunctions", "Markdown", "MaybeInplace", "PreallocationTools", "Preferences", "Printf", "RecursiveArrayTools", "SciMLBase", "SciMLJacobianOperators", "SciMLLogging", "SciMLOperators", "SciMLStructures", "Setfield", "StaticArraysCore", "SymbolicIndexingInterface", "TimerOutputs"]
-git-tree-sha1 = "4f595a0977d6e048fa1e3c382b088b950f8c7934"
+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"
 uuid = "be0214bd-f91f-a760-ac4e-3421ce2b2da0"
-version = "2.15.0"
+version = "2.20.0"
 
     [deps.NonlinearSolveBase.extensions]
     NonlinearSolveBaseBandedMatricesExt = "BandedMatrices"
@@ -1869,7 +1872,7 @@ version = "1.6.1"
 [[deps.OpenSSL_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "458c3c95-2e84-50aa-8efc-19380b2a3a95"
-version = "3.5.1+0"
+version = "3.5.4+0"
 
 [[deps.OpenSpecFun_jll]]
 deps = ["Artifacts", "CompilerSupportLibraries_jll", "JLLWrappers", "Libdl"]
@@ -1890,27 +1893,27 @@ version = "1.8.1"
 
 [[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"]
-git-tree-sha1 = "924e1db15095c7df6b844231c00c40d756a7553d"
+git-tree-sha1 = "47271adac597af08263b6ea2669e93040b45c4d0"
 uuid = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
-version = "6.108.0"
+version = "6.111.0"
 
 [[deps.OrdinaryDiffEqAdamsBashforthMoulton]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqLowOrderRK", "Polyester", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"]
-git-tree-sha1 = "8307937159c3aeec5f19f4b661d82d96d25a3ff1"
+git-tree-sha1 = "2e44acb684dfcdc2e41851a988733e30b28a8478"
 uuid = "89bda076-bce5-4f1c-845f-551c83cdda9a"
-version = "1.9.0"
+version = "1.11.0"
 
 [[deps.OrdinaryDiffEqBDF]]
 deps = ["ADTypes", "ArrayInterface", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqSDIRK", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "StaticArrays", "TruncatedStacktraces"]
-git-tree-sha1 = "22b0c4f7939af140b7f7f4ce2cce90d9f72fa515"
+git-tree-sha1 = "8ebd2eb20a8ebb12911b40282d5b41e2db6d107b"
 uuid = "6ad6398a-0878-4a85-9266-38940aa047c8"
-version = "1.22.0"
+version = "1.24.0"
 
 [[deps.OrdinaryDiffEqCore]]
-deps = ["ADTypes", "Accessors", "Adapt", "ArrayInterface", "ConcreteStructs", "DataStructures", "DiffEqBase", "DocStringExtensions", "EnumX", "EnzymeCore", "FastBroadcast", "FastClosures", "FastPower", "FillArrays", "FunctionWrappersWrappers", "InteractiveUtils", "LinearAlgebra", "Logging", "MacroTools", "MuladdMacro", "Polyester", "PrecompileTools", "Preferences", "Random", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLLogging", "SciMLOperators", "SciMLStructures", "Static", "StaticArrayInterface", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"]
-git-tree-sha1 = "b4a8d9b96931c2fc69126233bbe6d1a11b053d77"
+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", "StaticArraysCore", "SymbolicIndexingInterface", "TruncatedStacktraces"]
+git-tree-sha1 = "73e431f66fb9854474c0bec1c11053e8889629c1"
 uuid = "bbf590c4-e513-4bbe-9b18-05decba2e5d8"
-version = "3.22.0"
+version = "3.27.0"
 
     [deps.OrdinaryDiffEqCore.extensions]
     OrdinaryDiffEqCoreMooncakeExt = "Mooncake"
@@ -1922,15 +1925,15 @@ version = "3.22.0"
 
 [[deps.OrdinaryDiffEqDefault]]
 deps = ["ADTypes", "DiffEqBase", "EnumX", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqBDF", "OrdinaryDiffEqCore", "OrdinaryDiffEqRosenbrock", "OrdinaryDiffEqTsit5", "OrdinaryDiffEqVerner", "PrecompileTools", "Preferences", "Reexport", "SciMLBase"]
-git-tree-sha1 = "0f40d969dd10d1b226c864bf7dc4b4b8933bc130"
+git-tree-sha1 = "5fbd116f93790ce79c4beff8206e1b8e45c492a2"
 uuid = "50262376-6c5a-4cf5-baba-aaf4f84d72d7"
-version = "1.13.0"
+version = "1.14.0"
 
 [[deps.OrdinaryDiffEqDifferentiation]]
-deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqCore", "SciMLBase", "SciMLOperators", "SparseMatrixColorings", "StaticArrayInterface", "StaticArrays"]
-git-tree-sha1 = "c85968ea296acaff5de6ed0d9b64ddb00f4ea01f"
+deps = ["ADTypes", "ArrayInterface", "ConcreteStructs", "ConstructionBase", "DiffEqBase", "DifferentiationInterface", "FastBroadcast", "FiniteDiff", "ForwardDiff", "FunctionWrappersWrappers", "LinearAlgebra", "LinearSolve", "OrdinaryDiffEqCore", "SciMLBase", "SciMLOperators", "SparseMatrixColorings", "StaticArrays"]
+git-tree-sha1 = "c7d06493d3327cc45f0d8d7e641203d28779e837"
 uuid = "4302a76b-040a-498a-8c04-15b101fed76b"
-version = "2.2.1"
+version = "2.6.0"
 weakdeps = ["SparseArrays"]
 
     [deps.OrdinaryDiffEqDifferentiation.extensions]
@@ -1938,153 +1941,153 @@ weakdeps = ["SparseArrays"]
 
 [[deps.OrdinaryDiffEqExplicitRK]]
 deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"]
-git-tree-sha1 = "c5b900878b088776b8d1bd5a7b1d94e301e21c4e"
+git-tree-sha1 = "17d8bfeab9d6a7b31d116e950aa7204fd444ec19"
 uuid = "9286f039-9fbf-40e8-bf65-aa933bdc4db0"
-version = "1.9.0"
+version = "1.11.0"
 
 [[deps.OrdinaryDiffEqExponentialRK]]
 deps = ["ADTypes", "DiffEqBase", "ExponentialUtilities", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "RecursiveArrayTools", "Reexport", "SciMLBase"]
-git-tree-sha1 = "72156f954b199ff23dada0e8c0f12c44503b5cf9"
+git-tree-sha1 = "32fb54c8000c82fd32c55f7bf52f4651a3bd0e61"
 uuid = "e0540318-69ee-4070-8777-9e2de6de23de"
-version = "1.13.0"
+version = "1.15.0"
 
 [[deps.OrdinaryDiffEqExtrapolation]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "FastPower", "LinearSolve", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "Polyester", "RecursiveArrayTools", "Reexport", "SciMLBase"]
-git-tree-sha1 = "129730b7b6cb60cc9c18e0db5861f4a7ed2c30b9"
+git-tree-sha1 = "57bf89f5b144f6e5af1146fb7aee1544c6b76950"
 uuid = "becaefa8-8ca2-5cf9-886d-c06f3d2bd2c4"
-version = "1.16.0"
+version = "1.18.0"
 
 [[deps.OrdinaryDiffEqFIRK]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "FastGaussQuadrature", "FastPower", "LinearAlgebra", "LinearSolve", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "Polyester", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators"]
-git-tree-sha1 = "342c716e0c15ab44203f68a78f98800ec560df82"
+git-tree-sha1 = "b8af0f1b9925e5ae58a4adf5afc60269a3ac6c44"
 uuid = "5960d6e9-dd7a-4743-88e7-cf307b64f125"
-version = "1.23.0"
+version = "1.25.0"
 
 [[deps.OrdinaryDiffEqFeagin]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "Polyester", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"]
-git-tree-sha1 = "b123f64a8635a712ceb037a7d2ffe2a1875325d3"
+git-tree-sha1 = "302999c99bc454cf274d39b9ba3c2415f6fbb1cb"
 uuid = "101fe9f7-ebb6-4678-b671-3a81e7194747"
-version = "1.8.0"
+version = "1.10.0"
 
 [[deps.OrdinaryDiffEqFunctionMap]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"]
-git-tree-sha1 = "cbd291508808caf10cf455f974c2025e886ed2a3"
+git-tree-sha1 = "2d2d59f9e530bca2f99934ecfde5430f6ad54bfe"
 uuid = "d3585ca7-f5d3-4ba6-8057-292ed1abd90f"
-version = "1.9.0"
+version = "1.11.0"
 
 [[deps.OrdinaryDiffEqHighOrderRK]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"]
-git-tree-sha1 = "9584dcc90cf10216de7aa0f2a1edc0f54d254cf6"
+git-tree-sha1 = "46d3c3b871253c21b91fe30257eff55ce61c486a"
 uuid = "d28bc4f8-55e1-4f49-af69-84c1a99f0f58"
-version = "1.9.0"
+version = "1.11.0"
 
 [[deps.OrdinaryDiffEqIMEXMultistep]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "Reexport", "SciMLBase"]
-git-tree-sha1 = "9280abaf9ac36d60dd774113f7ce8a7f826d6e2e"
+git-tree-sha1 = "bd5440f7f09ce75df59e7f5d2f5b8e52d24aae3c"
 uuid = "9f002381-b378-40b7-97a6-27a27c83f129"
-version = "1.12.0"
+version = "1.14.0"
 
 [[deps.OrdinaryDiffEqLinear]]
 deps = ["DiffEqBase", "ExponentialUtilities", "LinearAlgebra", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "SciMLBase", "SciMLOperators"]
-git-tree-sha1 = "c92913fa5942ed9bc748f3e79a5c693c8ec0c3d7"
+git-tree-sha1 = "a1217ddfc979da3a2d9cf11b083688562e8ef333"
 uuid = "521117fe-8c41-49f8-b3b6-30780b3f0fb5"
-version = "1.10.0"
+version = "1.12.0"
 
 [[deps.OrdinaryDiffEqLowOrderRK]]
 deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"]
-git-tree-sha1 = "78223e34d4988070443465cd3f2bdc38d6bd14b0"
+git-tree-sha1 = "ee10d0fa457ee9805dcf1ab20e9a37958952c540"
 uuid = "1344f307-1e59-4825-a18e-ace9aa3fa4c6"
-version = "1.10.0"
+version = "1.12.0"
 
 [[deps.OrdinaryDiffEqLowStorageRK]]
 deps = ["Adapt", "DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static", "StaticArrays"]
-git-tree-sha1 = "bd032c73716bc538033af041ca8903df6c813bfd"
+git-tree-sha1 = "66ebd44d699322583305497eedd371877feb6679"
 uuid = "b0944070-b475-4768-8dec-fb6eb410534d"
-version = "1.12.0"
+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 = "a75727e93ffef0f0bc408372988f7bc0767b1781"
+git-tree-sha1 = "f6722d6a96c27263ac3ea6159a3174914567a807"
 uuid = "127b3ac7-2247-4354-8eb6-78cf4e7c58e8"
-version = "1.23.0"
+version = "1.25.0"
 
 [[deps.OrdinaryDiffEqNordsieck]]
 deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqTsit5", "Polyester", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"]
-git-tree-sha1 = "facea9aaf48eed5e9ba66d8b3246e51417c084d0"
+git-tree-sha1 = "21117475cb0d27942d31e4b8a27337392ae680b0"
 uuid = "c9986a66-5c92-4813-8696-a7ec84c806c8"
-version = "1.9.0"
+version = "1.11.0"
 
 [[deps.OrdinaryDiffEqPDIRK]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "Polyester", "Reexport", "SciMLBase", "StaticArrays"]
-git-tree-sha1 = "c95dd60623e11464e6079b77d2ce604fb399a02d"
+git-tree-sha1 = "274f4faaef2b137bdf05fb221a869d3d773aff22"
 uuid = "5dd0a6cf-3d4b-4314-aa06-06d4e299bc89"
-version = "1.11.0"
+version = "1.13.0"
 
 [[deps.OrdinaryDiffEqPRK]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "Polyester", "Reexport", "SciMLBase"]
-git-tree-sha1 = "baa77b7f874cda1f58f8c793fc7a9778e78a91c5"
+git-tree-sha1 = "1fe4f58f170c68b187ac42dc5f00932bb162ddfa"
 uuid = "5b33eab2-c0f1-4480-b2c3-94bc1e80bda1"
-version = "1.8.0"
+version = "1.10.0"
 
 [[deps.OrdinaryDiffEqQPRK]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static"]
-git-tree-sha1 = "9e351a8f923c843adb48945318437e051f6ee139"
+git-tree-sha1 = "275900d64ab7be5f1bca8c7bfd48166838b6ab44"
 uuid = "04162be5-8125-4266-98ed-640baecc6514"
-version = "1.8.0"
+version = "1.10.0"
 
 [[deps.OrdinaryDiffEqRKN]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "Polyester", "RecursiveArrayTools", "Reexport", "SciMLBase"]
-git-tree-sha1 = "b086c6d1b4153c9ff4b3f184a9ba7829413cc502"
+git-tree-sha1 = "edf3411246d6030e165a5fba513cbdcde5c72f79"
 uuid = "af6ede74-add8-4cfd-b1df-9a4dbb109d7a"
-version = "1.10.0"
+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 = "f11347f3f01a5b00dae2b565e73795ee138cdc68"
+git-tree-sha1 = "88f722a74df3dadaf5f31c454a280f1201a1b368"
 uuid = "43230ef6-c299-4910-a778-202eb28ce4ce"
-version = "1.25.0"
+version = "1.28.0"
 
 [[deps.OrdinaryDiffEqSDIRK]]
-deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"]
-git-tree-sha1 = "0b766d820e3b948881f1f246899de9ef3d329224"
+deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "LinearAlgebra", "MacroTools", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "TruncatedStacktraces"]
+git-tree-sha1 = "c5e26fbadbad137fa2b7e3161a290a1becfec502"
 uuid = "2d112036-d095-4a1e-ab9a-08536f3ecdbf"
-version = "1.12.0"
+version = "1.14.0"
 
 [[deps.OrdinaryDiffEqSSPRK]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static", "StaticArrays"]
-git-tree-sha1 = "8abc61382a0c6469aa9c3bff2d61c9925a088320"
+git-tree-sha1 = "66a8402e3e50eb07a8d4cb96c730eb8fc2bb372b"
 uuid = "669c94d9-1f4b-4b64-b377-1aa079aa2388"
-version = "1.11.0"
+version = "1.13.0"
 
 [[deps.OrdinaryDiffEqStabilizedIRK]]
 deps = ["ADTypes", "DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "OrdinaryDiffEqDifferentiation", "OrdinaryDiffEqNonlinearSolve", "OrdinaryDiffEqStabilizedRK", "RecursiveArrayTools", "Reexport", "SciMLBase", "StaticArrays"]
-git-tree-sha1 = "cf6856c731ddf9866e3e22612cce5e270f071545"
+git-tree-sha1 = "6688927070213dc11058ee2316027c64eb5052f7"
 uuid = "e3e12d00-db14-5390-b879-ac3dd2ef6296"
-version = "1.11.0"
+version = "1.13.0"
 
 [[deps.OrdinaryDiffEqStabilizedRK]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "RecursiveArrayTools", "Reexport", "SciMLBase", "StaticArrays"]
-git-tree-sha1 = "d156a972fa7bc37bf8377d33a7d51d152e354d4c"
+git-tree-sha1 = "bda2155c5abd7b1e4965d35f930f78e320d7bafe"
 uuid = "358294b1-0aab-51c3-aafe-ad5ab194a2ad"
-version = "1.8.0"
+version = "1.10.0"
 
 [[deps.OrdinaryDiffEqSymplecticRK]]
 deps = ["DiffEqBase", "FastBroadcast", "MuladdMacro", "OrdinaryDiffEqCore", "Polyester", "RecursiveArrayTools", "Reexport", "SciMLBase"]
-git-tree-sha1 = "9b783806fe2dc778649231cb3932cb71b63222d9"
+git-tree-sha1 = "b18a4e7973e73e84cfd79222e2389565b323b588"
 uuid = "fa646aed-7ef9-47eb-84c4-9443fc8cbfa8"
-version = "1.11.0"
+version = "1.13.0"
 
 [[deps.OrdinaryDiffEqTsit5]]
 deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static", "TruncatedStacktraces"]
-git-tree-sha1 = "8be4cba85586cd2efa6c76d1792c548758610901"
+git-tree-sha1 = "b305a813cc926ce31458e2b955612d9b719f1c38"
 uuid = "b1df2697-797e-41e3-8120-5422d3b24e4a"
-version = "1.9.0"
+version = "1.11.0"
 
 [[deps.OrdinaryDiffEqVerner]]
 deps = ["DiffEqBase", "FastBroadcast", "LinearAlgebra", "MuladdMacro", "OrdinaryDiffEqCore", "Polyester", "PrecompileTools", "Preferences", "RecursiveArrayTools", "Reexport", "SciMLBase", "Static", "TruncatedStacktraces"]
-git-tree-sha1 = "5ca5dbbfea89e14f283ce9fe2301c528ff4ec007"
+git-tree-sha1 = "918dbbaeea6684d4d65907994650d9614bdd2eb2"
 uuid = "79d7bb75-1356-48c1-b8c0-6832512096c2"
-version = "1.11.0"
+version = "1.13.0"
 
 [[deps.PCRE2_jll]]
 deps = ["Artifacts", "Libdl"]
@@ -2118,7 +2121,7 @@ version = "0.44.2+0"
 [[deps.Pkg]]
 deps = ["Artifacts", "Dates", "Downloads", "FileWatching", "LibGit2", "Libdl", "Logging", "Markdown", "Printf", "Random", "SHA", "TOML", "Tar", "UUIDs", "p7zip_jll"]
 uuid = "44cfe95a-1eb2-52ea-b672-e2afdf69b78f"
-version = "1.12.0"
+version = "1.12.1"
 weakdeps = ["REPL"]
 
     [deps.Pkg.extensions]
@@ -2176,9 +2179,9 @@ version = "1.4.3"
 
 [[deps.PreallocationTools]]
 deps = ["Adapt", "ArrayInterface", "PrecompileTools"]
-git-tree-sha1 = "dc8d6bde5005a0eac05ae8faf1eceaaca166cfa4"
+git-tree-sha1 = "e16b73bf892c55d16d53c9c0dbd0fb31cb7e25da"
 uuid = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
-version = "1.1.2"
+version = "1.2.0"
 weakdeps = ["ForwardDiff", "ReverseDiff", "SparseConnectivityTracer"]
 
     [deps.PreallocationTools.extensions]
@@ -2200,9 +2203,9 @@ version = "1.5.2"
 
 [[deps.PrettyTables]]
 deps = ["Crayons", "LaTeXStrings", "Markdown", "PrecompileTools", "Printf", "REPL", "Reexport", "StringManipulation", "Tables"]
-git-tree-sha1 = "211530a7dc76ab59087f4d4d1fc3f086fbe87594"
+git-tree-sha1 = "624de6279ab7d94fc9f672f0068107eb6619732c"
 uuid = "08abe8d2-0d0c-5749-adfa-8a2ac140af0d"
-version = "3.2.3"
+version = "3.3.2"
 
     [deps.PrettyTables.extensions]
     PrettyTablesTypstryExt = "Typstry"
@@ -2286,9 +2289,9 @@ version = "0.6.12"
 
 [[deps.RecursiveArrayTools]]
 deps = ["Adapt", "ArrayInterface", "DocStringExtensions", "GPUArraysCore", "LinearAlgebra", "PrecompileTools", "RecipesBase", "StaticArraysCore", "SymbolicIndexingInterface"]
-git-tree-sha1 = "18d2a6fd1ea9a8205cadb3a5704f8e51abdd748b"
+git-tree-sha1 = "e4fd3369c78666a65ccec25dba28a0b181434ab2"
 uuid = "731186ca-8d62-57ce-b412-fbd966d074cd"
-version = "3.48.0"
+version = "3.52.0"
 
     [deps.RecursiveArrayTools.extensions]
     RecursiveArrayToolsFastBroadcastExt = "FastBroadcast"
@@ -2370,9 +2373,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 = "0be0208add9b6836a701e0ac3ad30bda72fee51d"
+git-tree-sha1 = "908c0bf271604d09393a21c142116ab26f66f67c"
 uuid = "0bca4576-84f4-4d90-8ffe-ffa030f20462"
-version = "2.150.0"
+version = "2.154.0"
 
     [deps.SciMLBase.extensions]
     SciMLBaseChainRulesCoreExt = "ChainRulesCore"
@@ -2553,9 +2556,9 @@ version = "1.2.1"
 
 [[deps.SparseMatrixColorings]]
 deps = ["ADTypes", "DocStringExtensions", "LinearAlgebra", "PrecompileTools", "Random", "SparseArrays"]
-git-tree-sha1 = "fa43a02c01e3e3cb065c89bf9b648b89e3c06f18"
+git-tree-sha1 = "1c1be8c6fdfaf9b6c9e156c509e672953b8e6af7"
 uuid = "0a514795-09f3-496d-8182-132a7b665d35"
-version = "0.4.25"
+version = "0.4.26"
 
     [deps.SparseMatrixColorings.extensions]
     SparseMatrixColoringsCUDAExt = "CUDA"
@@ -2572,9 +2575,9 @@ version = "0.4.25"
 
 [[deps.SpecialFunctions]]
 deps = ["IrrationalConstants", "LogExpFunctions", "OpenLibm_jll", "OpenSpecFun_jll"]
-git-tree-sha1 = "5acc6a41b3082920f79ca3c759acbcecf18a8d78"
+git-tree-sha1 = "2700b235561b0335d5bef7097a111dc513b8655e"
 uuid = "276daf66-3868-5448-9aa4-cd146d93841b"
-version = "2.7.1"
+version = "2.7.2"
 weakdeps = ["ChainRulesCore"]
 
     [deps.SpecialFunctions.extensions]
@@ -2664,16 +2667,18 @@ version = "1.11.0"
 
 [[deps.StructUtils]]
 deps = ["Dates", "UUIDs"]
-git-tree-sha1 = "28145feabf717c5d65c1d5e09747ee7b1ff3ed13"
+git-tree-sha1 = "fa95b3b097bcef5845c142ea2e085f1b2591e92c"
 uuid = "ec057cc2-7a8d-4b58-b3b3-92acb9f63b42"
-version = "2.6.3"
+version = "2.7.1"
 
     [deps.StructUtils.extensions]
     StructUtilsMeasurementsExt = ["Measurements"]
+    StructUtilsStaticArraysCoreExt = ["StaticArraysCore"]
     StructUtilsTablesExt = ["Tables"]
 
     [deps.StructUtils.weakdeps]
     Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
+    StaticArraysCore = "1e83bf80-4336-4d27-bf5d-d5a4f845583c"
     Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 
 [[deps.StyledStrings]]
@@ -2800,9 +2805,9 @@ version = "0.4.1"
 
 [[deps.UnoSolver]]
 deps = ["LinearAlgebra", "OpenBLAS32_jll", "Uno_jll"]
-git-tree-sha1 = "cf39af9c8be02fd5f5027bc3184e8582b8f6d2a0"
+git-tree-sha1 = "a947341afb633f09cdcf875090981b80b37ec31f"
 uuid = "1baa60ac-02f7-4b39-a7a8-2f4f58486b05"
-version = "0.2.7"
+version = "0.2.10"
 
     [deps.UnoSolver.extensions]
     UnoSolverMathOptInterfaceExt = "MathOptInterface"
@@ -2814,14 +2819,14 @@ version = "0.2.7"
 
 [[deps.Uno_jll]]
 deps = ["ASL_jll", "Artifacts", "CompilerSupportLibraries_jll", "HSL_jll", "HiGHS_jll", "JLLWrappers", "LLVMOpenMP_jll", "Libdl", "METIS_jll", "MUMPS_seq_jll", "SPRAL_jll", "libblastrampoline_jll"]
-git-tree-sha1 = "30b1deeaeb5de7c0e0e4a7f9a253195147ed0e9e"
+git-tree-sha1 = "3530444c365f78f2d84af6224e4fffbd9e115e49"
 uuid = "396d5378-14f1-5ab1-981d-48acd51740ed"
-version = "2.5.0+0"
+version = "2.7.0+0"
 
 [[deps.UnsafeAtomics]]
-git-tree-sha1 = "b13c4edda90890e5b04ba24e20a310fbe6f249ff"
+git-tree-sha1 = "0f30765c32d66d58e41f4cb5624d4fc8a82ec13b"
 uuid = "013be700-e6cd-48c3-b4a1-df204f14c38f"
-version = "0.3.0"
+version = "0.3.1"
 weakdeps = ["LLVM"]
 
     [deps.UnsafeAtomics.extensions]
@@ -3078,9 +3083,9 @@ version = "1.28.1+0"
 
 [[deps.libpng_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Zlib_jll"]
-git-tree-sha1 = "e015f211ebb898c8180887012b938f3851e719ac"
+git-tree-sha1 = "e2a7072fc0cdd7949528c1455a3e5da4122e1153"
 uuid = "b53b4c65-9356-5827-b1ea-8c7a1a84506f"
-version = "1.6.55+0"
+version = "1.6.56+0"
 
 [[deps.libva_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Xorg_libX11_jll", "Xorg_libXext_jll", "Xorg_libXfixes_jll", "libdrm_jll"]
@@ -3112,9 +3117,9 @@ uuid = "1317d2d5-d96f-522e-a858-c73665f53c3e"
 version = "2022.0.0+1"
 
 [[deps.p7zip_jll]]
-deps = ["Artifacts", "Libdl"]
+deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "3f19e933-33d8-53b3-aaab-bd5110c3b7a0"
-version = "17.5.0+2"
+version = "17.7.0+0"
 
 [[deps.x264_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl"]
diff --git a/docs/src/example-control-free.md b/docs/src/example-control-free.md
new file mode 100644
index 00000000..dc3b3852
--- /dev/null
+++ b/docs/src/example-control-free.md
@@ -0,0 +1,162 @@
+# [Control-free problems](@id example-control-free)
+
+Control-free problems are optimal control problems without a control variable. They are used for **optimizing constant parameters in dynamical systems**, such as:
+
+- Identifying unknown parameters from observed data (parameter estimation)
+- Finding optimal parameters for a given performance criterion
+
+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
+using OptimalControl
+using NLPModelsIpopt
+using Plots
+```
+
+## Example 1: Exponential growth rate estimation
+
+Consider a system with exponential growth:
+
+```math
+\dot{x}(t) = p \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:
+
+```math
+\min_{p} \int_0^{10} (x(t) - x_{\text{obs}}(t))^2 \, dt
+```
+
+The analytical solution is $p = 0.5$.
+
+### [Problem definition](@id example-control-free-problem-1)
+
+```@example main
+# observed data (analytical solution)
+data(t) = 2.0 * exp(0.5 * t)
+
+# optimal control problem (parameter estimation)
+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
+end
+nothing # hide
+```
+
+### [Solution](@id example-control-free-solution-1)
+
+```@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")
+nothing # hide
+```
+
+```@example main
+plot(sol1; size=(800, 400))
+```
+
+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.
+
+## Example 2: Harmonic oscillator pulsation optimization
+
+Consider a harmonic oscillator:
+
+```math
+\ddot{q}(t) = -\omega^2 q(t)
+```
+
+with initial conditions $q(0) = 1$, $\dot{q}(0) = 0$ and final condition $q(1) = 0$. We want to find the minimal pulsation $\omega$ satisfying these constraints:
+
+```math
+    \begin{aligned}
+        & \text{Minimise} && \omega^2 \\
+        & \text{subject to} \\
+        & && \ddot{q}(t) = -\omega^2 q(t), \\[1.0em]
+        & && q(0) = 1, \quad \dot{q}(0) = 0, \\[0.5em]
+        & && q(1) = 0.
+    \end{aligned}
+```
+
+The analytical solution is $\omega = \pi/2 \approx 1.5708$, giving $q(t) = \cos(\pi t / 2)$.
+
+### [Problem definition](@id example-control-free-problem-2)
+
+```@example main
+# optimal control problem (pulsation optimization)
+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
+    q(1) == 0.0                  # final condition
+    
+    ẋ(t) == [v(t), -ω^2 * q(t)]
+    
+    ω^2 + 1e-5*∫(u(t)^2) → min   # minimize pulsation # TO REMOVE u(t) when possible
+end
+nothing # hide
+```
+
+### [Solution](@id example-control-free-solution-2)
+
+```@example main
+sol2 = solve(ocp2; display=false)
+println("Optimal pulsation: ω = ", variable(sol2))
+println("Objective value: ω² = ", objective(sol2))
+println("Expected: ω = π/2 ≈ 1.5708, ω² ≈ 2.4674")
+nothing # hide
+```
+
+```@example main
+plot(sol2; size=(800, 400))
+```
+
+The optimal pulsation should be close to $\omega = \pi/2 \approx 1.5708$, and the objective $\omega^2 \approx 2.4674$.
+
+## Comparison with analytical solutions
+
+For the harmonic oscillator, we can compare the numerical solution with the analytical one:
+
+```@example main
+# 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))
+plot!(plt, t_analytical, q_analytical; 
+      label="q (analytical)", linestyle=:dash, linewidth=2, subplot=1)
+plot!(plt, t_analytical, v_analytical; 
+      label="v (analytical)", linestyle=:dash, linewidth=2, subplot=2)
+```
+
+The numerical and analytical solutions should match very closely.
+
+!!! note "Applications"
+
+    Control-free problems appear in many contexts:
+    - **System identification**: estimating physical parameters (mass, damping, stiffness) from experimental data
+    - **Optimal design**: finding optimal geometric or physical parameters (length, stiffness, etc.)
+    - **Inverse problems**: reconstructing unknown inputs or initial conditions from partial observations
+    
+    See the [syntax documentation](@ref manual-abstract-control-free) for more details on defining control-free problems.
diff --git a/docs/src/example-double-integrator-energy.md b/docs/src/example-double-integrator-energy.md
index d35b3976..7d908e6d 100644
--- a/docs/src/example-double-integrator-energy.md
+++ b/docs/src/example-double-integrator-energy.md
@@ -49,8 +49,7 @@ ocp = @def begin
     x(t0) == x0
     x(tf) == xf
 
-    ∂(q)(t) == v(t)
-    ∂(v)(t) == u(t)
+    ẋ(t) == [v(t), u(t)]
 
     0.5∫( u(t)^2 ) → min
 end
@@ -68,8 +67,7 @@ nothing # hide
     \begin{aligned}
         & \text{Minimise} && \frac{1}{2}\int_0^1 u^2(t) \,\mathrm{d}t \\
         & \text{subject to} \\
-        & && \dot{q}(t) = v(t), \\[0.5em]
-        & && \dot{v}(t) = u(t), \\[1.0em]
+        & && \dot{x}(t) = [v(t), u(t)], \\[1.0em]
         & && x(0) = (-1,0), \\[0.5em] 
         & && x(1) = (0,0).
     \end{aligned}
@@ -200,8 +198,7 @@ ocp = @def begin
     x(t0) == x0
     x(tf) == xf
 
-    ∂(q)(t) == v(t)
-    ∂(v)(t) == u(t)
+    ẋ(t) == [v(t), u(t)]
 
     0.5∫( u(t)^2 ) → min
 end
diff --git a/docs/src/example-double-integrator-time.md b/docs/src/example-double-integrator-time.md
index 984fd081..06c6096c 100644
--- a/docs/src/example-double-integrator-time.md
+++ b/docs/src/example-double-integrator-time.md
@@ -18,9 +18,7 @@ This problem can be interpreted as a simple model for a wagon with constant mass
 <img src="./assets/chariot.svg" style="display: block; margin: 0 auto 20px auto;" width="400px">
 ```
 
-First, we need to import the [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl) package to define the 
-optimal control problem and [NLPModelsIpopt.jl](https://jso.dev/NLPModelsIpopt.jl) to solve it. 
-We also need to import the [Plots.jl](https://docs.juliaplots.org) package to plot the solution.
+First, we need to import the [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl) package to define the optimal control problem and [NLPModelsIpopt.jl](https://jso.dev/NLPModelsIpopt.jl) to solve it. We also need to import the [Plots.jl](https://docs.juliaplots.org) package to plot the solution.
 
 ```@example main
 using OptimalControl
@@ -38,24 +36,22 @@ Let us define the problem:
 ```
 
 ```@example main
-ocp = @def begin
+t0 = 0; x0 = [-1, 0]; xf = [0, 0]
 
+ocp = @def begin
     tf ∈ R,          variable
-    t ∈ [0, tf],     time
+    t ∈ [t0, tf],    time
     x = (q, v) ∈ R², state
     u ∈ R,           control
 
     -1 ≤ u(t) ≤ 1
 
-    q(0)  == -1
-    v(0)  == 0
-    q(tf) == 0
-    v(tf) == 0
+    x(t0) == x0
+    x(tf) == xf
 
-    ẋ(t) == [v(t), u(t)]
+    ẋ(t) == [v(t), u(t)]
 
     tf → min
-
 end
 nothing # hide
 ```
@@ -71,13 +67,10 @@ nothing # hide
     \begin{aligned}
         & \text{Minimise} && t_f \\[0.5em]
         & \text{subject to} \\[0.5em]
-        & && \dot q(t) = v(t), \\
-        & && \dot v(t) = u(t), \\[0.5em]
+        & && \dot x(t) = [v(t), u(t)], \\[0.5em]
         & && -1 \le u(t) \le 1, \\[0.5em]
-        & && q(0) = -1, \\[0.5em]
-        & && v(0) = 0, \\[0.5em]
-        & && q(t_f) = 0, \\[0.5em]
-        & && v(t_f) = 0.
+        & && x(0) = (-1, 0), \\[0.5em]
+        & && x(t_f) = (0, 0).
     \end{aligned}
 ```
 
@@ -151,9 +144,6 @@ nothing # hide
 The shooting function enforces the conditions:
 
 ```@example main
-t0 = 0
-x0 = [-1, 0]
-xf = [ 0, 0]
 function shoot!(s, p0, t1, tf) 
     x_t0, p_t0 = x0, p0
     x_t1, p_t1 = f_max(t0, x_t0, p_t0, t1)
diff --git a/docs/src/example-singular-control.md b/docs/src/example-singular-control.md
new file mode 100644
index 00000000..fb38d094
--- /dev/null
+++ b/docs/src/example-singular-control.md
@@ -0,0 +1,352 @@
+# [Singular control](@id example-singular-control)
+
+For control-affine systems of the form
+
+```math
+\dot{q}(t) = f_0(q(t)) + u(t) f_1(q(t)), \quad u(t) \in [u_{\min}, u_{\max}],
+```
+
+the pseudo-Hamiltonian is $H = H_0 + u H_1$, where $H_i(q, p) = \langle p, f_i(q) \rangle$ are the Hamiltonian lifts of the vector fields $f_0$ and $f_1$.
+
+When the **switching function** $H_1$ vanishes on a time interval (i.e., $H_1(q(t), p(t)) = 0$ for $t \in [t_1, t_2]$), the arc is called **singular**. On such arcs, the control cannot be determined directly from the maximization condition and must be computed by successive differentiation of $H_1$ along the flow.
+
+This page demonstrates how to compute singular controls both by hand and using differential geometry tools from OptimalControl.jl, then verifies the result numerically using direct and indirect methods.
+
+First, we import the necessary packages:
+
+```@example main
+using OptimalControl
+using NLPModelsIpopt
+using Plots
+```
+
+## Problem definition
+
+We consider a vehicle moving in the plane with drift. The state is $q = (x, y, \theta)$ where $(x, y)$ is the position and $\theta$ is the orientation. The dynamics are:
+
+```math
+\dot{x}(t) = \cos\theta(t), \quad \dot{y}(t) = \sin\theta(t) + x(t), \quad \dot{\theta}(t) = u(t),
+```
+
+with control constraint $u(t) \in [-1, 1]$.
+
+We want to find the time-optimal transfer from the origin $(0, 0)$ with free initial orientation to the target position $(1, 0)$ with free final orientation:
+
+```@example main
+ocp = @def begin
+    
+    tf ∈ R, variable
+    t ∈ [0, tf], time
+    q = (x, y, θ) ∈ R³, state
+    u ∈ R, control
+
+    -1 ≤ u(t) ≤ 1                     # Control bounds
+    -π/2 ≤ θ(t) ≤ π/2                 # State bounds (helps direct method convergence)
+
+    x(0) == 0
+    y(0) == 0
+    x(tf) == 1
+    y(tf) == 0
+
+    ∂(q)(t) == [cos(θ(t)), sin(θ(t)) + x(t), u(t)]
+
+    tf → min
+
+end
+nothing # hide
+```
+
+This is a control-affine system with:
+
+```math
+f_0(q) = \begin{pmatrix} \cos\theta \\ \sin\theta + x \\ 0 \end{pmatrix}, \quad
+f_1(q) = \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}.
+```
+
+## Direct method
+
+We solve the problem using a direct method:
+
+```@example main
+direct_sol = solve(ocp; display=false)
+println("Optimal time: tf = ", variable(direct_sol))
+nothing # hide
+```
+
+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...)
+```
+
+## Singular control by hand
+
+The pseudo-Hamiltonian for this time-optimal problem is:
+
+```math
+H(q, p, u) = p_1 \cos\theta + p_2(\sin\theta + x) + p_3 u.
+```
+
+This is control-affine: $H = H_0 + u H_1$ with:
+
+```math
+H_0(q, p) = p_1 \cos\theta + p_2(\sin\theta + x), \quad H_1(q, p) = p_3.
+```
+
+The switching function is $H_1 = p_3$. On a singular arc, we have $H_1 = 0$ and all its time derivatives must vanish.
+
+**First derivative:**
+
+```math
+\dot{H}_1 = \{H, H_1\} = \{H_0, H_1\} =: H_{01}.
+```
+
+Computing the Poisson bracket:
+
+```math
+H_{01} = \frac{\partial H_0}{\partial p_1} \frac{\partial H_1}{\partial x} - \frac{\partial H_0}{\partial x} \frac{\partial H_1}{\partial p_1}
+       + \frac{\partial H_0}{\partial p_2} \frac{\partial H_1}{\partial y} - \frac{\partial H_0}{\partial y} \frac{\partial H_1}{\partial p_2}
+       + \frac{\partial H_0}{\partial p_3} \frac{\partial H_1}{\partial \theta} - \frac{\partial H_0}{\partial \theta} \frac{\partial H_1}{\partial p_3}.
+```
+
+Since $H_1 = p_3$ depends only on $p_3$, the only non-zero contribution comes from the $(\theta, p_3)$ pair:
+
+```math
+H_{01} = \frac{\partial H_0}{\partial \theta} \frac{\partial H_1}{\partial p_3} - \frac{\partial H_0}{\partial p_3} \frac{\partial H_1}{\partial \theta} = (-p_1 \sin\theta + p_2 \cos\theta) \cdot 1 - 0 = -p_1 \sin\theta + p_2 \cos\theta.
+```
+
+On the singular arc, $H_{01} = 0$, which gives the constraint:
+
+```math
+p_2 \cos\theta = p_1 \sin\theta.
+```
+
+**Second derivative:**
+
+```math
+\dot{H}_{01} = \{H, H_{01}\} = \{H_0, H_{01}\} + u \{H_1, H_{01}\} =: H_{001} + u H_{101}.
+```
+
+For the arc to remain singular, $\dot{H}_{01} = 0$, which gives:
+
+```math
+u_s = -\frac{H_{001}}{H_{101}},
+```
+
+whenever $H_{101} \neq 0$. Computing $H_{001} = \{H_0, H_{01}\}$ with $H_{01} = -p_1 \sin\theta + p_2 \cos\theta$:
+
+The only non-zero contribution comes from the $(x, p_1)$ pair:
+
+```math
+H_{001} = \frac{\partial H_0}{\partial x} \frac{\partial H_{01}}{\partial p_1} - \frac{\partial H_0}{\partial p_1} \frac{\partial H_{01}}{\partial x} = p_2 \cdot (-\sin\theta) - \cos\theta \cdot 0 = -p_2 \sin\theta.
+```
+
+Computing $H_{101} = \{H_1, H_{01}\}$ with $H_1 = p_3$ and $H_{01} = -p_1 \sin\theta + p_2 \cos\theta$:
+
+The only non-zero contribution comes from the $(\theta, p_3)$ pair:
+
+```math
+H_{101} = \frac{\partial H_1}{\partial \theta} \frac{\partial H_{01}}{\partial p_3} - \frac{\partial H_1}{\partial p_3} \frac{\partial H_{01}}{\partial \theta} = 0 - 1 \cdot (-p_1 \cos\theta - p_2 \sin\theta) = p_1 \cos\theta + p_2 \sin\theta.
+```
+
+Therefore:
+
+```math
+u_s = -\frac{H_{001}}{H_{101}} = \frac{p_2 \sin\theta}{p_1 \cos\theta + p_2 \sin\theta}.
+```
+
+!!! note "Non-degeneracy condition"
+
+    We can show that $H_{101} \neq 0$ on the singular arc. From the constraint $p_1 \sin\theta = p_2 \cos\theta$, if we had $H_{101} = p_1 \cos\theta + p_2 \sin\theta = 0$, then:
+
+    ```math
+    \begin{pmatrix} \cos\theta & \sin\theta \\ -\sin\theta & \cos\theta \end{pmatrix}
+    \begin{pmatrix} p_1 \\ p_2 \end{pmatrix} = \begin{pmatrix} 0 \\ 0 \end{pmatrix}.
+    ```
+
+    Since this matrix has determinant 1 (hence is invertible), we would have $p_1 = p_2 = 0$. Combined with $p_3 = 0$ (from $H_1 = 0$), this gives $p = 0$, which is impossible for a time-minimization problem.
+
+**Simplification using the constraint:**
+
+Multiply numerator and denominator by $\sin\theta$:
+
+```math
+u_s = \frac{p_2 \sin^2\theta}{p_1 \cos\theta \sin\theta + p_2 \sin^2\theta}.
+```
+
+From the constraint $p_1 \sin\theta = p_2 \cos\theta$, we have $p_1 \cos\theta \sin\theta = p_2 \cos^2\theta$. Substituting in the denominator:
+
+```math
+u_s = \frac{p_2 \sin^2\theta}{p_2 \cos^2\theta + p_2 \sin^2\theta} = \frac{p_2 \sin^2\theta}{p_2(\cos^2\theta + \sin^2\theta)} = \sin^2\theta.
+```
+
+So the singular control is:
+
+```math
+u_s(\theta) = \sin^2\theta.
+```
+
+Let's overlay this on the numerical solution:
+
+```@example main
+T = time_grid(direct_sol, :control)
+θ(t) = state(direct_sol)(t)[3]
+us(t) = sin(θ(t))^2
+plot!(plt, T, us; subplot=7, line=:dash, lw=2, label="us (hand)")
+plot(plt[7]; size=(800, 400))
+```
+
+## Singular control via Poisson brackets
+
+We can compute the same result using the differential geometry tools from OptimalControl.jl. See the [differential geometry tools manual](@ref manual-differential-geometry) for detailed explanations.
+
+First, define the vector fields:
+
+```@example main
+F0(q) = [cos(q[3]), sin(q[3]) + q[1], 0]
+F1(q) = [0, 0, 1]
+nothing # hide
+```
+
+Compute their Hamiltonian lifts:
+
+```@example main
+H0 = Lift(F0)
+H1 = Lift(F1)
+nothing # hide
+```
+
+Compute the iterated Poisson brackets:
+
+```@example main
+H01 = @Lie {H0, H1}
+H001 = @Lie {H0, H01}
+H101 = @Lie {H1, H01}
+nothing # hide
+```
+
+The singular control is:
+
+```@example main
+us_bracket(q, p) = -H001(q, p) / H101(q, p)
+nothing # hide
+```
+
+Let's verify this gives the same result:
+
+```@example main
+q(t) = state(direct_sol)(t)
+p(t) = costate(direct_sol)(t)
+us_b(t) = us_bracket(q(t), p(t))
+plot!(plt, T, us_b; subplot=7, line=:dashdot, lw=2, label="us (brackets)")
+plot(plt[7]; size=(800, 400))
+```
+
+Both methods give the same singular control, which matches the numerical solution from the direct method.
+
+## Indirect shooting method
+
+We now solve the problem using an indirect shooting method based on the singular control we computed. This approach is similar to the one used in the [double integrator example](@ref example-double-integrator-energy).
+
+First, import the necessary packages:
+
+```@example main
+using OrdinaryDiffEq
+using NonlinearSolve
+```
+
+Define the singular control in feedback form:
+
+```@example main
+u_indirect(x) = sin(x[3])^2
+nothing # hide
+```
+
+Build the flow for the singular arc:
+
+```@example main
+f = Flow(ocp, (x, p, tf) -> u_indirect(x))
+nothing # hide
+```
+
+Define the shooting function. We have 5 unknowns: the initial costate $p_0 \in \mathbb{R}^3$, the initial orientation $\theta_0$, and the final time $t_f$. We must define 5 equations to solve for these unknowns.
+
+```@example main
+t0 = 0
+
+function shoot!(s, p0, θ0, tf)
+    q_t0, p_t0 = [0, 0, θ0], p0
+    q_tf, p_tf = f(t0, q_t0, p_t0, tf)
+    
+    s[1] = q_tf[1] - 1      # x(tf) = 1 (boundary condition)
+    s[2] = q_tf[2]          # y(tf) = 0 (boundary condition)
+    s[3] = p_t0[3]          # pθ(0) = 0 (transversality condition)
+    s[4] = p_tf[3]          # pθ(tf) = 0 (transversality condition)
+    
+    # H(tf) = 1 (for time-optimal with p^0 = -1)
+    pxf = p_tf[1]
+    pyf = p_tf[2]
+    θf = q_tf[3]
+    s[5] = pxf * cos(θf) + pyf * (sin(θf) + 1) - 1
+end
+nothing # hide
+```
+
+Use the direct solution to provide an initial guess:
+
+```@example main
+p0 = costate(direct_sol)(t0)
+θ0 = state(direct_sol)(t0)[3]
+tf = variable(direct_sol)
+
+println("Initial guess:")
+println("p0 = ", p0)
+println("θ0 = ", θ0)
+println("tf = ", tf)
+nothing # hide
+```
+
+Set up and solve the nonlinear system:
+
+```@example main
+# Auxiliary in-place NLE function
+nle!(s, ξ, _) = shoot!(s, ξ[1:3], ξ[4], ξ[5])
+
+# Initial guess for the Newton solver
+ξ_guess = [p0..., θ0, tf]
+
+# NLE problem with initial guess
+prob = NonlinearProblem(nle!, ξ_guess)
+
+# Resolution of the shooting equations
+shooting_sol = solve(prob; show_trace=Val(false))
+p0_sol, θ0_sol, tf_sol = shooting_sol.u[1:3], shooting_sol.u[4], shooting_sol.u[5]
+
+println("Shooting solution:")
+println("p0 = ", p0_sol)
+println("θ0 = ", θ0_sol)
+println("tf = ", tf_sol)
+nothing # hide
+```
+
+Reconstruct the indirect solution:
+
+```@example main
+indirect_sol = f((t0, tf_sol), [0, 0, θ0_sol], p0_sol; saveat=range(t0, tf_sol, 100))
+nothing # hide
+```
+
+Plot the indirect solution alongside the direct solution:
+
+```@example main
+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.
+
+## See also
+
+- [Differential geometry tools](@ref manual-differential-geometry) — Mathematical definitions and usage of `Lift`, `Poisson`, `@Lie`
+- [Goddard tutorial](@extref Tutorials tutorial-goddard) — More complex example with bang, singular, and boundary arcs
+- [Compute flows from optimal control problems](@ref manual-flow-ocp) — Using flows for indirect methods
diff --git a/docs/src/manual-abstract.md b/docs/src/manual-abstract.md
index 1986c883..93454693 100644
--- a/docs/src/manual-abstract.md
+++ b/docs/src/manual-abstract.md
@@ -24,7 +24,8 @@ end
 ```
 
 !!! warning
-    Note that the full code of the definition above is not provided (hence the `...`) The same is true for most examples below (only those without `...` are indeed complete). Also note that problem definitions must at least include definitions for time, state, control, dynamics and cost.
+    - Note that the full code of the definition above is not provided (hence the `...`) The same is true for most examples below (only those without `...` are indeed complete).
+    - Also note that problem definitions must at least include definitions for time, state, dynamics and cost. The control declaration is optional (see [Control-free problems](@ref manual-abstract-control-free)).
 
 Aliases `v₁`, `v₂` (and `v1`, `v2`) are automatically defined and can be used in subsequent expressions instead of `v[1]` and `v[2]`. The user can also define her own aliases for the components (one alias per dimension):
 
@@ -127,6 +128,45 @@ end
 !!! note
     One dimensional variable, state or control are treated as scalars (`Real`), not vectors (`Vector`). In Julia, for `x::Real`, it is possible to write `x[1]` (and `x[1][1]`...) so it is OK (though useless) to write `x₁`, `x1` or `x[1]` instead of simply `x` to access the corresponding value. Conversely it is *not* OK to use such an `x` as a vector, for instance as in `...f(x)...` where `f(x::Vector{T}) where {T <: Real}`.
 
+## [Control-free problems](@id manual-abstract-control-free)
+
+The control declaration is **optional**. You can define problems without control for:
+
+- **Parameter estimation**: Identify unknown parameters in the dynamics from observed data
+- **Dynamic optimization**: Optimize constant parameters subject to ODE constraints
+
+For example, to estimate a growth rate parameter:
+
+```julia
+@def begin
+    p ∈ R, variable              # parameter to estimate
+    t ∈ [0, 10], time
+    x ∈ R, state
+    x(0) == 2.0
+    ẋ(t) == p * x(t)             # dynamics depends on p
+    ∫(x(t) - data(t))² → min     # fit to observed data
+end
+```
+
+Or to optimize the pulsation of a harmonic oscillator:
+
+```julia
+@def begin
+    ω ∈ R, variable              # pulsation to optimize
+    t ∈ [0, 1], time
+    x = (q, v) ∈ R², state
+    q(0) == 1.0
+    v(0) == 0.0
+    q(1) == 0.0                  # final condition
+    ẋ(t) == [v(t), -ω²*q(t)]    # harmonic oscillator
+    ω² → min                     # minimize pulsation
+end
+```
+
+!!! compat "Upcoming feature"
+
+    Control-free problem syntax (omitting the control declaration) is currently being implemented. For now, use a dummy control with `u ∈ R, control` and `u(t) == 0` as a workaround. See the [Control-free problems example](@ref example-control-free) for executable examples.
+
 ## [Dynamics](@id manual-abstract-dynamics)
 
 ```julia
@@ -136,7 +176,7 @@ end
 The dynamics is given in the standard vectorial ODE form:
 
 ```math
-    \dot{x}(t) = f([t, ]x(t), u(t)[, v])
+    \dot{x}(t) = f([t, ]x(t)[, u(t)][, v])
 ```
 
 depending on whether it is autonomous / with a variable or not (the parser will detect time and variable dependences,
@@ -534,7 +574,7 @@ end
 
 ## Misc
 
-- Declarations (of variable - if any -, time, state and control) must be done first. Then, dynamics, constraints and cost can be introduced in an arbitrary order.
+- Declarations (of variable - if any -, time, state and control - if any -) must be done first. Then, dynamics, constraints and cost can be introduced in an arbitrary order.
 - It is possible to provide numbers / labels (as in math equations) for the constraints to improve readability (this is mostly for future use, typically to retrieve the Lagrange multiplier associated with the discretisation of a given constraint):
 
 ```julia
diff --git a/docs/src/manual-differential-geometry.md b/docs/src/manual-differential-geometry.md
new file mode 100644
index 00000000..1875f619
--- /dev/null
+++ b/docs/src/manual-differential-geometry.md
@@ -0,0 +1,634 @@
+# [Differential geometry tools](@id manual-differential-geometry)
+
+Optimal control theory relies on differential geometry tools to analyze Hamiltonian systems, compute singular controls, study controllability, and more. This page introduces the main operators available in OptimalControl.jl: Hamiltonian lift, Lie derivatives, Poisson brackets, Lie brackets, and partial time derivatives.
+
+!!! note "Type qualification"
+
+    Types like `Hamiltonian`, `HamiltonianLift`, `VectorField`, and `HamiltonianVectorField` are **not exported** by OptimalControl.jl. You must qualify them with `OptimalControl.` when using them (e.g., `OptimalControl.VectorField`). Functions and operators (`Lift`, `⋅`, `Lie`, `Poisson`, `@Lie`, `∂ₜ`) are exported and can be used directly.
+
+First, import the package:
+
+```@example main
+using OptimalControl
+```
+
+## Hamiltonian lift
+
+Given a vector field $X: \mathbb{R}^n \to \mathbb{R}^n$, its **Hamiltonian lift** is the function $H_X: \mathbb{R}^n \times (\mathbb{R}^n)^* \to \mathbb{R}$ defined by
+
+```math
+H_X(x, p) = \langle p, X(x) \rangle = \sum_{i=1}^n p_i X_i(x).
+```
+
+### From plain Julia functions
+
+The simplest way to compute a Hamiltonian lift is from a plain Julia function. By default, the function is treated as **autonomous** (time-independent) and **non-variable** (no extra parameter):
+
+```@example main
+# Define a vector field as a Julia function
+X(x) = [x[2], -x[1]]
+
+# Compute its Hamiltonian lift
+H = Lift(X)
+
+# Evaluate at a point (x, p)
+x = [1, 2]
+p = [3, 4]
+H(x, p)
+```
+
+The result is $H(x, p) = p_1 x_2 + p_2 (-x_1) = 3 \times 2 + 4 \times (-1) = 2$.
+
+### From VectorField type
+
+You can also use the `OptimalControl.VectorField` type, which allows more control over the function's properties:
+
+```@example main-1
+using OptimalControl # hide
+# Wrap in VectorField (autonomous, non-variable by default)
+X = OptimalControl.VectorField(x -> [x[2], -x[1]])
+H = Lift(X)
+
+# This returns a HamiltonianLift object
+H([1, 2], [3, 4])
+```
+
+### Non-autonomous case
+
+For time-dependent vector fields, use `autonomous=false`:
+
+```@example main-2
+using OptimalControl # hide
+# Non-autonomous vector field: X(t, x) = [t*x[2], -x[1]]
+X(t, x) = [t * x[2], -x[1]]
+H = Lift(X; autonomous=false)
+
+# Signature is now H(t, x, p)
+H(2, [1, 2], [3, 4])
+```
+
+### Variable case
+
+For vector fields depending on an additional parameter $v$, use `variable=true`:
+
+```@example main-3
+using OptimalControl # hide
+# Variable vector field: X(x, v) = [x[2] + v, -x[1]]
+X(x, v) = [x[2] + v, -x[1]]
+H = Lift(X; variable=true)
+
+# Signature is now H(x, p, v)
+H([1, 2], [3, 4], 1)
+```
+
+## Lie derivative
+
+The **Lie derivative** of a function $f: \mathbb{R}^n \to \mathbb{R}$ along a vector field $X$ is defined by
+
+```math
+(X \cdot f)(x) = f'(x) \cdot X(x) = \sum_{i=1}^n \frac{\partial f}{\partial x_i}(x) X_i(x).
+```
+
+This represents the directional derivative of $f$ along $X$.
+
+### [From plain Julia functions](@id lie-from-functions)
+
+When using plain Julia functions, they are treated as autonomous and non-variable:
+
+```@example main-4
+using OptimalControl # hide
+# Vector field and scalar function
+X(x) = [x[2], -x[1]]
+f(x) = x[1]^2 + x[2]^2
+
+# Lie derivative (using dot operator)
+Xf = X ⋅ f
+
+# Evaluate at a point
+Xf([1, 2])
+```
+
+For the harmonic oscillator with $X(x) = (x_2, -x_1)$ and energy $f(x) = x_1^2 + x_2^2$:
+
+```math
+(X \cdot f)(x) = 2x_1 x_2 + 2x_2(-x_1) = 0,
+```
+
+which confirms that energy is conserved along trajectories.
+
+### [From VectorField type](@id lie-from-vectorfield)
+
+```@example main-5
+using OptimalControl # hide
+# Using VectorField type
+X = OptimalControl.VectorField(x -> [x[2], -x[1]])
+g(x) = x[1]^2 + x[2]^2
+
+# Lie derivative
+Xg = X ⋅ g
+Xg([1, 2])
+```
+
+### Alternative syntax
+
+The `Lie` function is equivalent to the `⋅` operator:
+
+```@example main-5
+# These are equivalent
+Xg1 = X ⋅ g
+Xg2 = Lie(X, g)
+
+Xg1([1, 2]) == Xg2([1, 2])
+```
+
+### With keyword arguments
+
+For non-autonomous or variable cases, use the `Lie` function with keyword arguments:
+
+```@example main-6
+using OptimalControl # hide
+# Non-autonomous case
+X(t, x) = [t + x[2], -x[1]]
+f(t, x) = t + x[1]^2 + x[2]^2
+
+Xf = Lie(X, f; autonomous=false)
+Xf(1, [1, 2])
+```
+
+```@example main-7
+using OptimalControl # hide
+# Variable case
+X(x, v) = [x[2] + v, -x[1]]
+f(x, v) = x[1]^2 + x[2]^2 + v
+
+Xf = Lie(X, f; variable=true)
+Xf([1, 2], 1)
+```
+
+### With VectorField type
+
+You can also create the VectorField explicitly with the keywords, then use it without keywords in the Lie function:
+
+```@example main-7a
+using OptimalControl # hide
+# Non-autonomous VectorField created with keywords
+X = OptimalControl.VectorField((t, x) -> [t + x[2], -x[1]]; autonomous=false)
+f(t, x) = t + x[1]^2 + x[2]^2
+
+# No keywords needed here - the VectorField already knows its properties
+Xf = Lie(X, f)
+Xf(1, [1, 2])
+```
+
+```@example main-7b
+using OptimalControl # hide
+# Variable VectorField created with keywords
+X = OptimalControl.VectorField((x, v) -> [x[2] + v, -x[1]]; variable=true)
+f(x, v) = x[1]^2 + x[2]^2 + v
+
+# No keywords needed here
+Xf = Lie(X, f)
+Xf([1, 2], 1)
+```
+
+## Poisson bracket
+
+For two functions $f, g: \mathbb{R}^n \times (\mathbb{R}^n)^* \to \mathbb{R}$, the **Poisson bracket** is defined by
+
+```math
+\{f, g\}(x, p) = \sum_{i=1}^n \left( \frac{\partial f}{\partial p_i} \frac{\partial g}{\partial x_i} - \frac{\partial f}{\partial x_i} \frac{\partial g}{\partial p_i} \right).
+```
+
+### Properties
+
+The Poisson bracket satisfies:
+
+- **Bilinearity**: $\{af + bg, h\} = a\{f, h\} + b\{g, h\}$ for scalars $a, b$
+- **Antisymmetry**: $\{f, g\} = -\{g, f\}$
+- **Leibniz rule**: $\{fg, h\} = f\{g, h\} + g\{f, h\}$
+- **Jacobi identity**: $\{\{f, g\}, h\} + \{\{h, f\}, g\} + \{\{g, h\}, f\} = 0$
+
+### [From plain Julia functions](@id poisson-from-functions)
+
+```@example main-8
+using OptimalControl # hide
+# Define two Hamiltonian functions
+f(x, p) = p[1] * x[2] + p[2] * x[1]
+g(x, p) = x[1]^2 + p[2]^2
+
+# Compute the Poisson bracket
+H = Poisson(f, g)
+
+# Evaluate at a point
+x = [1, 2]
+p = [3, 4]
+H(x, p)
+```
+
+### Verify antisymmetry
+
+```@example main-8
+Hfg = Poisson(f, g)
+Hgf = Poisson(g, f)
+
+println("Poisson(f, g) = ", Hfg(x, p))
+println("Poisson(g, f) = ", Hgf(x, p))
+println("Sum = ", Hfg(x, p) + Hgf(x, p))
+```
+
+### From Hamiltonian type
+
+```@example main-8
+# Wrap in Hamiltonian type
+F = OptimalControl.Hamiltonian(f)
+G = OptimalControl.Hamiltonian(g)
+
+H = Poisson(F, G)
+H(x, p)
+```
+
+### [With keyword arguments](@id poisson-kwargs)
+
+```@example main-9
+using OptimalControl # hide
+# Non-autonomous case
+f(t, x, p) = t + p[1] * x[2] + p[2] * x[1]
+g(t, x, p) = t^2 + x[1]^2 + p[2]^2
+
+H = Poisson(f, g; autonomous=false)
+H(1, [1, 2], [3, 4])
+```
+
+### With Hamiltonian type and keywords
+
+You can also create Hamiltonian objects explicitly with keywords, then use them without keywords in the Poisson function. **Important**: both Hamiltonians must have the same time and variable dependencies:
+
+```@example main-9a
+using OptimalControl # hide
+# Non-autonomous Hamiltonians created with keywords
+f(t, x, p) = t + p[1] * x[2] + p[2] * x[1]
+g(t, x, p) = t^2 + x[1]^2 + p[2]^2
+
+F = OptimalControl.Hamiltonian(f; autonomous=false)
+G = OptimalControl.Hamiltonian(g; autonomous=false)
+
+# No keywords needed here - both Hamiltonians are already non-autonomous
+H = Poisson(F, G)
+H(1, [1, 2], [3, 4])
+```
+
+```@example main-9b
+using OptimalControl # hide
+# Variable Hamiltonians created with keywords
+f(x, p, v) = x[1]^2 + p[2]^2 + v
+g(x, p, v) = x[2]^2 + p[1]^2 + 2*v
+
+F = OptimalControl.Hamiltonian(f; variable=true)
+G = OptimalControl.Hamiltonian(g; variable=true)
+
+# Both are variable, so the Poisson bracket is also variable
+H = Poisson(F, G)
+H([1, 2], [3, 4], 1)
+```
+
+### Relation to Hamiltonian vector fields
+
+The Poisson bracket is closely related to the Lie derivative. If $\vec{H} = (\nabla_p H, -\nabla_x H)$ denotes the Hamiltonian vector field associated to $H$, then
+
+```math
+\{H, G\} = \vec{H} \cdot G.
+```
+
+This means the Poisson bracket of $H$ and $G$ equals the Lie derivative of $G$ along the Hamiltonian vector field of $H$.
+
+### Poisson bracket of Hamiltonian lifts
+
+When computing the Poisson bracket of two Hamiltonian lifts, the result is the Hamiltonian lift of the Lie bracket of the underlying vector fields:
+
+```@example main-10
+using OptimalControl # hide
+# Two vector fields
+X(x) = [x[1]^2, x[2]^2]
+Y(x) = [x[2], -x[1]]
+
+# Their Hamiltonian lifts
+HX = Lift(X)
+HY = Lift(Y)
+
+# Poisson bracket of lifts
+H = Poisson(HX, HY)
+H([1, 2], [3, 4])
+```
+
+This satisfies: $\{H_X, H_Y\} = H_{[X,Y]}$ where $[X,Y]$ is the Lie bracket of vector fields (see next section).
+
+## Lie bracket of vector fields
+
+For two vector fields $X, Y: \mathbb{R}^n \to \mathbb{R}^n$, the **Lie bracket** is the vector field $[X, Y]$ defined by
+
+```math
+[X, Y](x) = Y'(x) \cdot X(x) - X'(x) \cdot Y(x),
+```
+
+where $X'(x)$ denotes the Jacobian matrix of $X$ at $x$.
+
+### [From VectorField type](@id bracket-from-vectorfield)
+
+```@example main-11
+using OptimalControl # hide
+# Define two vector fields
+X = OptimalControl.VectorField(x -> [x[2], -x[1]])
+Y = OptimalControl.VectorField(x -> [x[1], x[2]])
+
+# Compute the Lie bracket
+Z = Lie(X, Y)
+
+# Evaluate at a point
+Z([1, 2])
+```
+
+### Relation to Poisson brackets
+
+If $H_X = \text{Lift}(X)$ and $H_Y = \text{Lift}(Y)$ are the Hamiltonian lifts, then:
+
+```math
+\{H_X, H_Y\} = H_{[X,Y]}.
+```
+
+Let's verify this numerically:
+
+```@example main-11
+# Hamiltonian lifts
+HX = Lift(x -> X(x))
+HY = Lift(x -> Y(x))
+
+# Poisson bracket of the lifts
+HXY = Poisson(HX, HY)
+
+# Lift of the Lie bracket
+HZ = Lift(x -> Z(x))
+
+# Compare at a point
+x = [1, 2]
+p = [3, 4]
+
+println("Poisson bracket: ", HXY(x, p))
+println("Lift of Lie bracket: ", HZ(x, p))
+```
+
+## The `@Lie` macro
+
+The `@Lie` macro provides a convenient syntax for computing Lie brackets (for vector fields) and Poisson brackets (for Hamiltonians).
+
+!!! warning "Important distinction"
+
+    - **Square brackets `[...]`** denote **Lie brackets** and work with:
+        - `VectorField` objects
+        - Plain Julia functions (automatically wrapped as `VectorField`)
+    - **Curly braces `{...}`** denote **Poisson brackets** and work with:
+        - Plain Julia functions (automatically wrapped as `Hamiltonian`)
+        - `Hamiltonian` objects
+
+    When using **only plain functions** (no `VectorField` or `Hamiltonian` objects), specify `autonomous` and/or `variable` keywords as needed to match your function signature. Keywords are optional - if not specified, they use the default values (`autonomous=true`, `variable=false`). If you mix plain functions with `VectorField` or `Hamiltonian` objects, the keywords are inferred from the `VectorField` or `Hamiltonian` objects.
+
+### Lie brackets with VectorField
+
+```@example main-12
+using OptimalControl # hide
+# Define vector fields
+F1 = OptimalControl.VectorField(x -> [0, -x[3], x[2]])
+F2 = OptimalControl.VectorField(x -> [x[3], 0, -x[1]])
+
+# Compute Lie bracket using macro
+F12 = @Lie [F1, F2]
+
+# Evaluate
+F12([1, 2, 3])
+```
+
+### Nested Lie brackets
+
+```@example main-12
+F3 = OptimalControl.VectorField(x -> [x[1], x[2], x[3]])
+F123 = @Lie [[F1, F2], F3]
+F123([1, 2, 3])
+```
+
+### Lie brackets with plain Julia functions
+
+You can also use plain Julia functions directly with the `@Lie` macro. The functions will be automatically wrapped in `VectorField` objects:
+
+```@example main-12a
+using OptimalControl # hide
+# Define plain Julia functions
+X(x) = [x[2], -x[1]]
+Y(x) = [x[1], x[2]]
+
+# Compute Lie bracket using macro with plain functions
+Z = @Lie [X, Y]
+
+# Evaluate
+Z([1, 2])
+```
+
+### With keyword arguments for plain functions
+
+For non-autonomous or variable cases, specify the keywords:
+
+```@example main-12b
+using OptimalControl # hide
+# Non-autonomous plain functions
+X(t, x) = [t + x[2], -x[1]]
+Y(t, x) = [x[1], t*x[2]]
+
+# Use autonomous=false keyword
+Z = @Lie [X, Y] autonomous=false
+Z(1, [1, 2])
+```
+
+```@example main-12c
+using OptimalControl # hide
+# Variable plain functions
+X(x, v) = [x[2] + v, -x[1]]
+Y(x, v) = [x[1], x[2] + v]
+
+# Use variable=true keyword
+Z = @Lie [X, Y] variable=true
+Z([1, 2], 1)
+```
+
+### Nested brackets with plain functions
+
+```@example main-12d
+using OptimalControl # hide
+X(x) = [0, -x[3], x[2]]
+Y(x) = [x[3], 0, -x[1]]
+Z(x) = [x[1], x[2], x[3]]
+
+# Nested Lie brackets
+nested = @Lie [[X, Y], Z]
+nested([1, 2, 3])
+```
+
+!!! tip "Plain functions vs VectorField"
+
+    Using plain functions with `@Lie [X, Y]` is convenient for quick computations. However, if you need to reuse the same vector field multiple times or want explicit control over the autonomy/variability, consider creating `VectorField` objects explicitly:
+
+    ```julia
+    # Explicit VectorField (keywords at creation)
+    X = OptimalControl.VectorField((t, x) -> [t + x[2], -x[1]]; autonomous=false)
+    Y = OptimalControl.VectorField((t, x) -> [x[1], t*x[2]]; autonomous=false)
+    Z = @Lie [X, Y]  # No keywords needed
+
+    # Plain functions (keywords at macro call)
+    X(t, x) = [t + x[2], -x[1]]
+    Y(t, x) = [x[1], t*x[2]]
+    Z = @Lie [X, Y] autonomous=false
+    ```
+
+### Poisson brackets from plain functions
+
+For Hamiltonian functions (plain Julia functions), use curly braces `{_, _}`:
+
+```@example main-13
+using OptimalControl # hide
+# Define Hamiltonian functions
+H0(x, p) = p[1] * x[2] + p[2] * (-x[1])
+H1(x, p) = p[2]
+
+# Compute Poisson bracket
+H01 = @Lie {H0, H1}
+
+# Evaluate
+H01([1, 2], [3, 4])
+```
+
+### Iterated Poisson brackets
+
+The macro is particularly useful for computing iterated brackets, which appear in singular control analysis:
+
+```@example main-13
+# First-order bracket
+H01 = @Lie {H0, H1}
+
+# Second-order brackets
+H001 = @Lie {H0, H01}
+H101 = @Lie {H1, H01}
+
+# Evaluate
+x = [1, 2]
+p = [3, 4]
+
+println("H01(x, p) = ", H01(x, p))
+println("H001(x, p) = ", H001(x, p))
+println("H101(x, p) = ", H101(x, p))
+```
+
+These iterated brackets are used to compute singular controls. For a pseudo-Hamiltonian of the form $H = H_0 + u H_1$, if the switching function $H_1$ vanishes on an interval (singular arc), the control is given by
+
+```math
+u_s = -\frac{H_{001}}{H_{101}},
+```
+
+provided $H_{101} \neq 0$. See the [singular control example](@ref example-singular-control) for a complete application.
+
+### [With keyword arguments](@id macro-kwargs)
+
+For non-autonomous functions, specify `autonomous=false`:
+
+```@example main-14
+using OptimalControl # hide
+# Non-autonomous Hamiltonians
+H0(t, x, p) = t + p[1] * x[2] + p[2] * (-x[1])
+H1(t, x, p) = p[2]
+
+# Poisson bracket with keyword
+H01 = @Lie {H0, H1} autonomous=false
+
+# Evaluate
+H01(1, [1, 2], [3, 4])
+```
+
+### Poisson brackets from Hamiltonian type
+
+```@example main-15
+using OptimalControl # hide
+# Using Hamiltonian type
+H1 = OptimalControl.Hamiltonian((x, p) -> x[1]^2 + p[2]^2)
+H2 = OptimalControl.Hamiltonian((x, p) -> x[2]^2 + p[1]^2)
+
+# Macro works with Hamiltonian objects too
+H12 = @Lie {H1, H2}
+H12([1, 1], [3, 2])
+```
+
+## Partial time derivative
+
+For non-autonomous functions $f(t, x, \ldots)$, the **partial derivative with respect to time** is computed using the `∂ₜ` operator:
+
+```math
+(\partial_t f)(t, x, \ldots) = \frac{\partial f}{\partial t}(t, x, \ldots).
+```
+
+### Basic usage
+
+```@example main-16
+using OptimalControl # hide
+# Define a time-dependent function
+f(t, x) = t * x
+
+# Compute partial time derivative
+df = ∂ₜ(f)
+
+# Evaluate
+df(0, 8)
+```
+
+```@example main-16
+df(2, 3)
+```
+
+### More complex example
+
+```@example main-17
+using OptimalControl # hide
+# Function with multiple arguments
+g(t, x, p) = t^2 + x[1] * p[1] + x[2] * p[2]
+
+# Partial derivative
+dg = ∂ₜ(g)
+
+# At t=3, ∂g/∂t = 2t = 6
+dg(3, [1, 2], [4, 5])
+```
+
+### Relation to total time derivative
+
+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\}.
+```
+
+This is the sum of:
+
+- The **partial time derivative** $\partial_t G$ (explicit time dependence)
+- The **Poisson bracket** $\{H, G\}$ (evolution along the flow)
+
+This relation is fundamental in non-autonomous optimal control theory.
+
+## Summary
+
+| Function/Operator | Mathematical notation | Julia syntax |
+| ----------------- | --------------------- | ------------ |
+| Hamiltonian lift | $H_X(x,p) = \langle p, X(x) \rangle$ | `H = Lift(X)` |
+| Lie derivative | $(X \cdot f)(x) = f'(x) \cdot X(x)$ | `X ⋅ f` or `Lie(X, f)` |
+| Poisson bracket | $\{f,g\}(x,p) = \langle \nabla_p f, \nabla_x g \rangle - \langle \nabla_x f, \nabla_p g \rangle$ | `Poisson(f, g)` or `@Lie {f, g}` |
+| Lie bracket | $[X,Y](x) = Y'(x) X(x) - X'(x) Y(x)$ | `Lie(X, Y)` or `@Lie [X, Y]` |
+| Partial time derivative | $\partial_t f(t, x, \ldots)$ | `∂ₜ(f)` |
+
+## See also
+
+- [Compute flows from Hamiltonians and others](@ref manual-flow-others) — Using flows with Hamiltonian vector fields
+- [Singular control example](@ref example-singular-control) — Application to computing singular controls
+- [Goddard tutorial](@extref Tutorials tutorial-goddard) — Complex example with bang, singular, and boundary arcs
diff --git a/docs/src/manual-initial-guess.md b/docs/src/manual-initial-guess.md
index 43c76255..e970eca2 100644
--- a/docs/src/manual-initial-guess.md
+++ b/docs/src/manual-initial-guess.md
@@ -63,6 +63,20 @@ nothing # hide
 We first solve the problem without giving an initial guess.
 This will default to initialize all variables to 0.1.
 
+To visualize the default initial guess before solving, we can run the solver with `max_iter=0`:
+
+```@example main
+# visualize the default initial guess (no iterations)
+sol_init = solve(ocp1; init=nothing, max_iter=0, display=false)
+plot(sol_init; size=(600, 450))
+```
+
+!!! tip "Visualizing any initial guess"
+
+    This technique works with any initial guess specification. By setting `max_iter=0`, the solver stops immediately after initialization, allowing you to visualize the initial guess before the optimization process begins.
+
+Now let us solve the problem completely:
+
 ```@example main
 # solve the optimal control problem without initial guess
 sol = solve(ocp1; display=false)
diff --git a/docs/src/manual-model.md b/docs/src/manual-model.md
index a36a6dd9..43eb3173 100644
--- a/docs/src/manual-model.md
+++ b/docs/src/manual-model.md
@@ -70,7 +70,7 @@ where $p^0 = -1$ since we are in the normal case. From the Pontryagin maximum pr
 u(x, p) = p_v
 ```
 
-since $\partial^2_{uu} H = p^0 = - 1 < 0$. 
+since $\partial^2_{uu} H = p^0 = - 1 < 0$.
 
 ```@example main
 u = (x, p) -> p[2]          # control law in feedback form
@@ -231,7 +231,7 @@ is_autonomous(ocp)  # false if dynamics or cost depend on time
 
 For more details on autonomy, see the [Time dependence](@ref manual-model-time-dependence) section below.
 
-#### Summary table
+#### [Summary table](@id manual-model-summary-time)
 
 | Method | Returns | Description |
 | -------- | --------- | --------- |
@@ -293,7 +293,7 @@ Get the dimension of state box constraints:
 dim_state_constraints_box(ocp)  # returns number of box constraints on state
 ```
 
-#### Summary table
+#### [Summary table](@id manual-model-summary-state)
 
 | Method | Returns | Description |
 | -------- | --------- | --------- |
@@ -344,7 +344,7 @@ Get the dimension of control box constraints:
 dim_control_constraints_box(ocp)  # returns number of box constraints on control
 ```
 
-#### Summary table
+#### [Summary table](@id manual-model-summary-control)
 
 | Method | Returns | Description |
 | -------- | --------- | --------- |
@@ -395,7 +395,7 @@ Get the dimension of variable box constraints:
 dim_variable_constraints_box(ocp)  # returns number of box constraints on variable
 ```
 
-#### Summary table
+#### [Summary table](@id manual-model-summary-variable)
 
 | Method | Returns | Description |
 | -------- | --------- | --------- |
@@ -431,7 +431,7 @@ The first argument `dx` is mutated upon call and contains the state derivative.
 * `u`: control
 * `v`: variable
 
-#### Summary table
+#### [Summary table](@id manual-model-summary-dynamics)
 
 | Method | Returns | Description |
 | -------- | --------- | --------- |
@@ -493,7 +493,7 @@ v = [1.0, 2.0]
 f⁰(s, q, u, v)  # returns the integrand value
 ```
 
-#### Summary table
+#### [Summary table](@id manual-model-summary-objective)
 
 | Method | Returns | Description |
 | -------- | --------- | --------- |
@@ -595,7 +595,7 @@ dim_boundary_constraints_nl(ocp)  # number of nonlinear boundary constraints
 
     To get the dual variable (or Lagrange multiplier) associated to a constraint, use the [`dual`](@ref) method on a solution.
 
-#### Summary table
+#### [Summary table](@id manual-model-summary-constraints)
 
 | Method | Returns | Description |
 | -------- | --------- | --------- |
diff --git a/docs/src/manual-plot.md b/docs/src/manual-plot.md
index dc02e0c0..796ddc9d 100644
--- a/docs/src/manual-plot.md
+++ b/docs/src/manual-plot.md
@@ -22,12 +22,12 @@ The table below summarizes the main plotting arguments and links to the correspo
 
 | Section                                             | Relevant Arguments                                                                            |
 | :---------------------------------------------------| :-------------------------------------------------------------------------------------------- |
-| [Basic concepts](@ref manual-plot-basic)          | `size`, `state_style`, `costate_style`, `control_style`, `time_style`, `kwargs...`            |
-| [Split vs. group layout](@ref manual-plot-layout) | `layout`                                                                                      |
-| [Plotting control norm](@ref manual-plot-control) | `control`                                                                                     |
-| [Normalised time](@ref manual-plot-time)          | `time`                                                                                        |
-| [Constraints](@ref manual-plot-constraints)       | `state_bounds_style`, `control_bounds_style`, `path_style`, `path_bounds_style`, `dual_style` |
-| [What to plot](@ref manual-plot-select)           | `description...`                                                                              |
+| [Basic concepts](@ref manual-plot-basic)            | `size`, `state_style`, `costate_style`, `control_style`, `time_style`, `kwargs...`            |
+| [Split vs. group layout](@ref manual-plot-layout)   | `layout`                                                                                      |
+| [Plotting control norm](@ref manual-plot-control)   | `control`                                                                                     |
+| [Normalised time](@ref manual-plot-time)            | `time`                                                                                        |
+| [Constraints](@ref manual-plot-constraints)         | `state_bounds_style`, `control_bounds_style`, `path_style`, `path_bounds_style`, `dual_style` |
+| [What to plot](@ref manual-plot-select)             | `description...`                                                                              |
 
 You can plot solutions obtained from the `solve` function or from a flow computed using an optimal control problem and a control law. See the [Basic Concepts](@ref manual-plot-basic) and [From Flow function](@ref manual-plot-flow) sections for details.
 
@@ -104,6 +104,7 @@ plotattr("color") # Specific Attribute Example
 You can also visit the Plot documentation online to get the descriptions of the attributes:
 
 - To pass attributes to the plot, see the [attributes plot](https://docs.juliaplots.org/latest/generated/attributes_plot/) documentation. For instance, you can specify the size of the figure.
+
 ```@raw html
 <details style="margin-left:3em"><summary>List of plot attributes.</summary>
 ```
@@ -117,7 +118,9 @@ end # hide
 ```@raw html
 </details>
 ```
+
 - You can pass attributes to all subplots at once by referring to the [attributes subplot](https://docs.juliaplots.org/latest/generated/attributes_subplot/) documentation. For example, you can specify the location of the legends.
+
 ```@raw html
 <details style="margin-left:3em"><summary>List of subplot attributes.</summary>
 ```
@@ -131,7 +134,9 @@ end # hide
 ```@raw html
 </details>
 ```
+
 - Similarly, you can pass axis attributes to all subplots. See the [attributes axis](https://docs.juliaplots.org/latest/generated/attributes_axis/) documentation. For example, you can remove the grid from every subplot.
+
 ```@raw html
 <details style="margin-left:3em"><summary>List of axis attributes.</summary>
 ```
@@ -145,7 +150,9 @@ end # hide
 ```@raw html
 </details>
 ```
+
 - Finally, you can pass series attributes to all subplots. Refer to the [attributes series](https://docs.juliaplots.org/latest/generated/attributes_series/) documentation. For instance, you can set the width of the curves using `linewidth`.
+
 ```@raw html
 <details style="margin-left:3em"><summary>List of series attributes.</summary>
 ```
@@ -174,7 +181,7 @@ plot(sol;
      control_style = (color=:red, linewidth=2))       # style: control trajectory
 ```
 
-Vertical axes at the initial and final times are automatically plotted. The style can me modified with the `time_style` keyword argument. 
+Vertical axes at the initial and final times are automatically plotted. The style can me modified with the `time_style` keyword argument.
 Additionally, you can choose not to display for instance the state and the costate trajectories by setting their styles to `:none`. You can set to `:none` any style.
 
 ```@example main
@@ -235,7 +242,7 @@ If you prefer to get a more compact figure, you can use the `layout` optional ke
 ```@example main
 plot(sol; layout=:group)
 ```
-    
+
 The default layout value is `:split` which corresponds to the grid of subplots presented above.
 
 ```@example main
@@ -316,6 +323,7 @@ plot(plt[1]) # x₁
 plt = plot(sol)
 plot(plt[2]) # x₂
 ```
+
 ```@example main
 plt = plot(sol)
 plot(plt[3]) # p₁
diff --git a/docs/src/manual-solution.md b/docs/src/manual-solution.md
index 770dd669..2714f562 100644
--- a/docs/src/manual-solution.md
+++ b/docs/src/manual-solution.md
@@ -9,7 +9,7 @@ After covering these core functionalities, we'll delve into the **structure of a
 
 ---
 
-**Content**
+## Content
 
 * [Main functionalities](@ref manual-solution-main-functionalities)
 * [Solution struct](@ref manual-solution-struct)
@@ -195,7 +195,7 @@ times(sol)  # returns the TimesModel struct containing time information
 !!! note "Trajectory data formats"
     Trajectories (`state`, `control`, `costate`, `path_constraints_dual`) can be provided either as matrices (rows = time points, columns = components) or as functions `t -> vector` for interpolated or analytical data.
 
-#### Summary table
+#### [Summary table](@id manual-solution-summary-trajectories)
 
 | Method | Returns | Description |
 | -------- | --------- | ------------- |
@@ -218,7 +218,7 @@ Get the optimal objective value:
 objective(sol)  # returns the objective value
 ```
 
-#### Summary table
+#### [Summary table](@id manual-solution-summary-objective)
 
 | Method | Returns | Description |
 | -------- | --------- | ------------- |
@@ -375,7 +375,7 @@ Get additional solver-specific information:
 infos(sol)  # returns dictionary of solver info
 ```
 
-#### Summary table
+#### [Summary table](@id manual-solution-summary-solver)
 
 | Method | Returns | Description |
 | -------- | --------- | ------------- |
diff --git a/docs/src/manual-solve-explicit.md b/docs/src/manual-solve-explicit.md
index 7faa5d1f..665f4922 100644
--- a/docs/src/manual-solve-explicit.md
+++ b/docs/src/manual-solve-explicit.md
@@ -42,8 +42,7 @@ The mode is **automatically detected**: if any of `discretizer`, `modeler`, or `
 
 ### Creating strategy instances
 
-Each strategy is constructed with its options as keyword arguments. 
-First, load the required solver packages:
+Each strategy is constructed with its options as keyword arguments. First, load the required solver packages:
 
 ```julia
 # Load solver packages (only what you need)
diff --git a/src/imports/ctflows.jl b/src/imports/ctflows.jl
index 8ddeedc3..7a61d6e5 100644
--- a/src/imports/ctflows.jl
+++ b/src/imports/ctflows.jl
@@ -4,7 +4,7 @@
 @reexport import CTFlows: CTFlows # for generated code (prefix)
 
 # Types
-import CTFlows: Hamiltonian, HamiltonianLift, HamiltonianVectorField
+import CTFlows: Hamiltonian, HamiltonianLift, HamiltonianVectorField, VectorField
 
 # Methods
-@reexport import CTFlows: Lift, Flow, ⋅, Lie, Poisson, @Lie, *
+@reexport import CTFlows: Lift, Flow, ⋅, Lie, Poisson, @Lie, *, ∂ₜ
diff --git a/test/helpers/print_utils.jl b/test/helpers/print_utils.jl
index 537ac9ea..7db896da 100644
--- a/test/helpers/print_utils.jl
+++ b/test/helpers/print_utils.jl
@@ -182,7 +182,9 @@ end
         ref_obj::Union{Real, Nothing},
         iterations::Union{Int, Nothing} = nothing,
         memory_bytes::Union{Int, Nothing} = nothing,
-        show_memory::Bool = false
+        show_memory::Bool = false,
+        rtol::Real = 1e-2,
+        atol::Real = 1e-3
     )
 
 Display a formatted table line for a test result.
@@ -215,6 +217,8 @@ Format inspired by print_benchmark_line() in CTBenchmarks.jl.
 - `iterations`: Number of iterations (optional)
 - `memory_bytes`: Allocated memory in bytes (optional)
 - `show_memory`: Show memory (default: false)
+- `rtol`: Relative tolerance for color thresholds (default: 1e-2)
+- `atol`: Absolute tolerance for color thresholds (default: 1e-3)
 """
 function print_test_line(
     test_type::String,
@@ -229,9 +233,17 @@ function print_test_line(
     iterations::Union{Int,Nothing}=nothing,
     memory_bytes::Union{Int,Nothing}=nothing,
     show_memory::Bool=false,
+    rtol::Real=1e-2,
+    atol::Real=1e-3,
 )
-    # Relative error calculation
-    rel_error = ref_obj === nothing ? missing : abs(obj - ref_obj) / abs(ref_obj) * 100
+    # Error calculation: use absolute error when ref is near zero, otherwise relative error
+    error_value = if ref_obj === nothing
+        missing
+    elseif abs(ref_obj) < 1e-6
+        abs(obj - ref_obj)  # Absolute error for near-zero reference
+    else
+        abs(obj - ref_obj) / abs(ref_obj) * 100  # Relative error in %
+    end
 
     # Colored status (✓ green or ✗ red)
     if success
@@ -289,12 +301,27 @@ function print_test_line(
     print(lpad(ref_str, 14))
     print(" │ ")
 
-    # Error: scientific notation with 2 decimal places
-    err_str = rel_error === missing ? "N/A" : @sprintf("%.2e", rel_error / 100)  # Convert to fraction then scientific format
-    err_color = if rel_error === missing
+    # Error: display format depends on whether it's absolute or relative
+    err_str = if error_value === missing
+        "N/A"
+    elseif ref_obj !== nothing && abs(ref_obj) < 1e-6
+        # Absolute error: display as scientific notation
+        @sprintf("%.2e", error_value)
+    else
+        # Relative error: display as scientific notation (convert % to fraction)
+        @sprintf("%.2e", error_value / 100)
+    end
+    
+    # Color thresholds based on provided tolerances
+    err_color = if error_value === missing
         :white
+    elseif ref_obj !== nothing && abs(ref_obj) < 1e-6
+        # Absolute error: green if < atol, yellow if < 5*atol, red otherwise
+        (error_value < atol ? :green : (error_value < 5 * atol ? :yellow : :red))
     else
-        (rel_error < 1 ? :green : (rel_error < 5 ? :yellow : :red))
+        # Relative error (in %): green if < rtol*100, yellow if < 5*rtol*100, red otherwise
+        rtol_percent = rtol * 100
+        (error_value < rtol_percent ? :green : (error_value < 5 * rtol_percent ? :yellow : :red))
     end
     printstyled(lpad(err_str, 7); color=err_color)
 
diff --git a/test/problems/TestProblems.jl b/test/problems/TestProblems.jl
index 810755ce..a60069f9 100644
--- a/test/problems/TestProblems.jl
+++ b/test/problems/TestProblems.jl
@@ -7,9 +7,11 @@ include("goddard.jl")
 include("double_integrator.jl")
 include("quadrotor.jl")
 include("transfer.jl")
+include("control_free.jl")
 
 export Beam, Goddard
 export DoubleIntegratorTime, DoubleIntegratorEnergy, DoubleIntegratorEnergyConstrained
 export Quadrotor, Transfer
+export ExponentialGrowth, HarmonicOscillator
 
 end
diff --git a/test/problems/control_free.jl b/test/problems/control_free.jl
new file mode 100644
index 00000000..19ce173e
--- /dev/null
+++ b/test/problems/control_free.jl
@@ -0,0 +1,94 @@
+# Control-free optimal control problems for testing parameter estimation
+# and dynamic optimization capabilities.
+
+using OptimalControl
+
+"""
+    ExponentialGrowth()
+
+Return data for the exponential growth rate estimation problem.
+
+The problem consists in estimating the growth rate parameter `p` for:
+- ẋ(t) = p * x(t)
+- x(0) = 2.0
+
+by minimizing the squared error with observed data x_obs(t) = 2.0 * exp(0.5 * t):
+- ∫₀¹⁰ (x(t) - x_obs(t))² dt → min
+
+The function returns a NamedTuple with fields:
+  * `ocp`  – CTParser/@def optimal control problem
+  * `obj`  – reference optimal objective value (≈0.0, perfect fit)
+  * `name` – short problem name
+  * `p_expected` – expected parameter value (0.5)
+"""
+function ExponentialGrowth()
+    # observed data (analytical solution)
+    data(t) = 2.0 * exp(0.5 * t)
+
+    @def ocp 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
+    end
+
+    return (
+        ocp=ocp,
+        obj=0.0,  # perfect fit expected
+        name="exponential_growth",
+        init=nothing,
+        p_expected=0.5,
+    )
+end
+
+"""
+    HarmonicOscillator()
+
+Return data for the harmonic oscillator pulsation optimization problem.
+
+The problem consists in finding the minimal pulsation ω for:
+- q̈(t) = -ω² * q(t)
+- q(0) = 1.0, v(0) = 0.0
+- q(1) = 0.0
+
+by minimizing ω²:
+- ω² → min
+
+The analytical solution is ω = π/2 ≈ 1.5708.
+
+The function returns a NamedTuple with fields:
+  * `ocp`  – CTParser/@def optimal control problem
+  * `obj`  – reference optimal objective value (π²/4 ≈ 2.4674)
+  * `name` – short problem name
+  * `ω_expected` – expected pulsation value (π/2)
+"""
+function HarmonicOscillator()
+    @def ocp 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
+        q(1) == 0.0                  # final condition
+
+        ẋ(t) == [v(t), -ω^2 * q(t)]
+
+        ω^2 + 1e-5*∫(u(t)^2) → min   # minimize pulsation # TO REMOVE u(t) when possible
+    end
+
+    return (
+        ocp=ocp,
+        obj=π^2 / 4,  # ω² = (π/2)² = π²/4 ≈ 2.4674
+        name="harmonic_oscillator",
+        init=nothing,
+        ω_expected=π / 2,  # ≈ 1.5708
+    )
+end
diff --git a/test/suite/reexport/test_ctflows.jl b/test/suite/reexport/test_ctflows.jl
index 3e5fcfbf..66ca73f3 100644
--- a/test/suite/reexport/test_ctflows.jl
+++ b/test/suite/reexport/test_ctflows.jl
@@ -23,6 +23,7 @@ function test_ctflows()
                 OptimalControl.Hamiltonian,
                 OptimalControl.HamiltonianLift,
                 OptimalControl.HamiltonianVectorField,
+                OptimalControl.VectorField,
             )
                 Test.@test isdefined(OptimalControl, nameof(T))
                 Test.@test !isdefined(CurrentModule, nameof(T))
@@ -37,7 +38,7 @@ function test_ctflows()
             end
         end
         Test.@testset "Operators" begin
-            for op in (:⋅, :Lie, :Poisson, :*)
+            for op in (:⋅, :Lie, :Poisson, :*, :∂ₜ)
                 Test.@test isdefined(OptimalControl, op)
                 Test.@test isdefined(CurrentModule, op)
             end
@@ -128,7 +129,7 @@ function test_ctflows()
             end
 
             Test.@testset "@Lie Macro" begin
-                # Test basic macro usage
+                # Test basic macro usage with VectorField
                 X1 = CTFlows.VectorField(x -> [x[2], -x[1]])
                 X2 = CTFlows.VectorField(x -> [x[1], x[2]])
 
@@ -140,6 +141,108 @@ function test_ctflows()
                 Test.@test lie_macro_result([1, 2]) isa Vector
             end
 
+            Test.@testset "@Lie Macro with Plain Functions" begin
+                # ================================================================
+                # Autonomous plain functions
+                # ================================================================
+                Test.@testset "Autonomous plain functions" begin
+                    X(x) = [x[2], -x[1]]
+                    Y(x) = [x[1], x[2]]
+
+                    # Lie bracket with macro
+                    Z = @Lie [X, Y]
+                    Test.@test Z isa CTFlows.VectorField
+                    Test.@test Z([1, 2]) isa Vector
+
+                    # Should give same result as with VectorField objects
+                    X_vf = CTFlows.VectorField(X)
+                    Y_vf = CTFlows.VectorField(Y)
+                    Z_vf = @Lie [X_vf, Y_vf]
+                    Test.@test Z([1, 2]) ≈ Z_vf([1, 2])
+                end
+
+                # ================================================================
+                # Non-autonomous plain functions
+                # ================================================================
+                Test.@testset "Non-autonomous plain functions" begin
+                    X(t, x) = [t + x[2], -x[1]]
+                    Y(t, x) = [x[1], t * x[2]]
+
+                    Z = @Lie [X, Y] autonomous = false
+                    Test.@test Z isa CTFlows.VectorField
+                    Test.@test Z(1, [1, 2]) isa Vector
+
+                    # Verify against VectorField version
+                    X_vf = CTFlows.VectorField(X; autonomous=false)
+                    Y_vf = CTFlows.VectorField(Y; autonomous=false)
+                    Z_vf = @Lie [X_vf, Y_vf]
+                    Test.@test Z(1, [1, 2]) ≈ Z_vf(1, [1, 2])
+                end
+
+                # ================================================================
+                # Variable plain functions
+                # ================================================================
+                Test.@testset "Variable plain functions" begin
+                    X(x, v) = [x[2] + v, -x[1]]
+                    Y(x, v) = [x[1], x[2] + v]
+
+                    Z = @Lie [X, Y] variable = true
+                    Test.@test Z isa CTFlows.VectorField
+                    Test.@test Z([1, 2], 1) isa Vector
+
+                    # Verify against VectorField version
+                    X_vf = CTFlows.VectorField(X; variable=true)
+                    Y_vf = CTFlows.VectorField(Y; variable=true)
+                    Z_vf = @Lie [X_vf, Y_vf]
+                    Test.@test Z([1, 2], 1) ≈ Z_vf([1, 2], 1)
+                end
+
+                # ================================================================
+                # Non-autonomous + variable plain functions
+                # ================================================================
+                Test.@testset "Non-autonomous variable plain functions" begin
+                    X(t, x, v) = [t + x[2] + v, -x[1]]
+                    Y(t, x, v) = [x[1], t * x[2] + v]
+
+                    Z = @Lie [X, Y] autonomous = false variable = true
+                    Test.@test Z isa CTFlows.VectorField
+                    Test.@test Z(1, [1, 2], 1) isa Vector
+
+                    # Verify against VectorField version
+                    X_vf = CTFlows.VectorField(X; autonomous=false, variable=true)
+                    Y_vf = CTFlows.VectorField(Y; autonomous=false, variable=true)
+                    Z_vf = @Lie [X_vf, Y_vf]
+                    Test.@test Z(1, [1, 2], 1) ≈ Z_vf(1, [1, 2], 1)
+                end
+
+                # ================================================================
+                # Nested Lie brackets with plain functions
+                # ================================================================
+                Test.@testset "Nested brackets with plain functions" begin
+                    X(x) = [x[2], -x[1]]
+                    Y(x) = [x[1], x[2]]
+                    Z_func(x) = [0, x[1]]
+
+                    # [[X, Y], Z]
+                    nested = @Lie [[X, Y], Z_func]
+                    Test.@test nested isa CTFlows.VectorField
+                    Test.@test nested([1, 2]) isa Vector
+                end
+
+                # ================================================================
+                # Mixed: plain function + VectorField
+                # ================================================================
+                Test.@testset "Mixed plain function and VectorField" begin
+                    X(x) = [x[2], -x[1]]
+                    Y_vf = CTFlows.VectorField(x -> [x[1], x[2]])
+
+                    # Should work with one plain function and one VectorField
+                    Z = @Lie [X, Y_vf]
+                    Test.@test Z isa CTFlows.VectorField
+                    Test.@test Z([1, 2]) isa Vector
+                end
+            end
+
             Test.@testset "Complex Signature Tests" begin
                 # Test with different arities and keyword arguments
 
@@ -222,6 +325,37 @@ function test_ctflows()
                 result3 = @Lie {H1, H2}
                 Test.@test result3 isa CTFlows.Hamiltonian
             end
+
+            Test.@testset "VectorField Construction" begin
+                # Test VectorField construction patterns
+                X1 = OptimalControl.VectorField(x -> [x[2], -x[1]])
+                Test.@test X1 isa OptimalControl.VectorField
+                Test.@test X1([1, 2]) isa Vector
+
+                # Non-autonomous
+                X2 = OptimalControl.VectorField((t, x) -> [t + x[2], -x[1]]; autonomous=false)
+                Test.@test X2 isa OptimalControl.VectorField
+                Test.@test X2(1.0, [1, 2]) isa Vector
+
+                # Variable
+                X3 = OptimalControl.VectorField((x, v) -> [x[2] + v, -x[1]]; variable=true)
+                Test.@test X3 isa OptimalControl.VectorField
+                Test.@test X3([1, 2], 1.0) isa Vector
+            end
+
+            Test.@testset "∂ₜ Operator" begin
+                # Test partial time derivative
+                f = (t, x) -> t * x
+                df = ∂ₜ(f)
+                Test.@test df isa Function
+                Test.@test df(0, 8) ≈ 8
+                Test.@test df(2, 3) ≈ 3
+
+                # More complex function
+                g = (t, x, p) -> t^2 + x[1]*p[1]
+                dg = ∂ₜ(g)
+                Test.@test dg(3, [1, 2], [4, 5]) ≈ 6
+            end
         end
     end
 end
diff --git a/test/suite/solve/test_canonical.jl b/test/suite/solve/test_canonical.jl
index 57e3af21..414bdac0 100644
--- a/test/suite/solve/test_canonical.jl
+++ b/test/suite/solve/test_canonical.jl
@@ -32,6 +32,7 @@ const SHOWTIMING = isdefined(Main, :TestOptions) ? Main.TestOptions.SHOWTIMING :
 
 # Objective tolerance for comparison with reference values
 const OBJ_RTOL = 1e-2
+const OBJ_ATOL = 1e-3  # Absolute tolerance for small objectives
 
 # CUDA availability check
 is_cuda_on() = CUDA.functional()
@@ -97,6 +98,8 @@ function run_test(
             iters,
             memory_bytes > 0 ? memory_bytes : nothing,
             false,  # show_memory = false
+            OBJ_RTOL,  # rtol for color thresholds
+            OBJ_ATOL,  # atol for color thresholds
         )
     end
 
@@ -111,7 +114,12 @@ function run_test(
         Test.@test success
         if success
             Test.@test solve_result isa OptimalControl.AbstractSolution
-            Test.@test OptimalControl.objective(solve_result) ≈ pb.obj rtol = OBJ_RTOL
+            # Use absolute tolerance when reference objective is near zero
+            if abs(pb.obj) < 1e-6
+                Test.@test OptimalControl.objective(solve_result) ≈ pb.obj atol = OBJ_ATOL
+            else
+                Test.@test OptimalControl.objective(solve_result) ≈ pb.obj rtol = OBJ_RTOL
+            end
         end
     end
 end
@@ -191,13 +199,15 @@ function test_canonical()
         ]
 
         problems = [
+            # ("Quadrotor", Quadrotor()), # some DomainError some times
+            # ("Transfer", Transfer()), # debug: add later (currently issue with :exa)
             ("Beam", Beam()),
             ("Goddard", Goddard()),
-            # ("Quadrotor", Quadrotor()), # some DomainError some times
             ("DI_Time", DoubleIntegratorTime()),
             ("DI_Energy", DoubleIntegratorEnergy()),
             ("DI_EnergyCons", DoubleIntegratorEnergyConstrained()),
-            # ("Transfer", Transfer()), # debug: add later (currently issue with :exa)
+            ("ExpGrowth", ExponentialGrowth()),
+            ("HarmOsc", HarmonicOscillator()),
         ]
 
         # Use Refs for statistics to pass to helper functions