From 5d85c86a570ef96df88756726bd2f0f377a7f624 Mon Sep 17 00:00:00 2001
From: abavoil <17646791+abavoil@users.noreply.github.com>
Date: Wed, 15 Apr 2026 11:40:59 +0200
Subject: [PATCH 01/21] symbolics first commit

---
 docs/Project.toml              |   2 +
 docs/make.jl                   |   3 +-
 docs/src/assets/Manifest.toml  | 241 ++++++++++++++++++++++++++++++---
 docs/src/assets/Project.toml   |   2 +
 docs/src/tutorial-symbolics.md |   6 +
 5 files changed, 235 insertions(+), 19 deletions(-)
 create mode 100644 docs/src/tutorial-symbolics.md

diff --git a/docs/Project.toml b/docs/Project.toml
index 7966b8a..0d45e28 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -16,6 +16,7 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 BenchmarkTools = "1"
@@ -35,4 +36,5 @@ OrdinaryDiffEq = "6"
 Plots = "1"
 Printf = "1"
 Suppressor = "0.2"
+Symbolics = "7"
 julia = "1.10"
diff --git a/docs/make.jl b/docs/make.jl
index 47fc79f..041ed0c 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -48,7 +48,7 @@ repo_url = "github.com/control-toolbox/Tutorials.jl"
 Draft = false
 ```
 =#
-draft = false  # Draft mode: if true, @example blocks in markdown are not executed
+draft = true  # Draft mode: if true, @example blocks in markdown are not executed
 
 # ═══════════════════════════════════════════════════════════════════════════════
 # Build documentation
@@ -85,6 +85,7 @@ makedocs(;
             "Linear–quadratic regulator" => "tutorial-lqr.md",
             "Minimal action" => "tutorial-mam.md",
             "Model Predictive Control" => "tutorial-mpc.md",
+            "Symbolics Mechanics" => "tutorial-symbolics.md",
         ],
     ],
     plugins=[links],
diff --git a/docs/src/assets/Manifest.toml b/docs/src/assets/Manifest.toml
index c65e39b..65c19ad 100644
--- a/docs/src/assets/Manifest.toml
+++ b/docs/src/assets/Manifest.toml
@@ -1,8 +1,8 @@
 # This file is machine-generated - editing it directly is not advised
 
-julia_version = "1.12.5"
+julia_version = "1.12.6"
 manifest_format = "2.0"
-project_hash = "a947cf68efaf3e1dbdfceb29223045209d5bb275"
+project_hash = "1d73fc8b7306dbb42eb2a1fdd2ae0c8be50a6d19"
 
 [[deps.ADNLPModels]]
 deps = ["ADTypes", "ForwardDiff", "LinearAlgebra", "NLPModels", "Requires", "ReverseDiff", "SparseArrays", "SparseConnectivityTracer", "SparseMatrixColorings"]
@@ -38,6 +38,12 @@ git-tree-sha1 = "6252039f98492252f9e47c312c8ffda0e3b9e78d"
 uuid = "ae81ac8f-d209-56e5-92de-9978fef736f9"
 version = "0.1.3+0"
 
+[[deps.AbstractPlutoDingetjes]]
+deps = ["Pkg"]
+git-tree-sha1 = "6e1d2a35f2f90a4bc7c2ed98079b2ba09c35b83a"
+uuid = "6e696c72-6542-2067-7265-42206c756150"
+version = "1.3.2"
+
 [[deps.AbstractTrees]]
 git-tree-sha1 = "2d9c9a55f9c93e8887ad391fbae72f8ef55e1177"
 uuid = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
@@ -88,6 +94,12 @@ version = "1.1.3"
 uuid = "0dad84c5-d112-42e6-8d28-ef12dabb789f"
 version = "1.1.2"
 
+[[deps.ArnoldiMethod]]
+deps = ["LinearAlgebra", "Random", "StaticArrays"]
+git-tree-sha1 = "d57bd3762d308bded22c3b82d033bff85f6195c6"
+uuid = "ec485272-7323-5ecc-a04f-4719b315124d"
+version = "0.4.0"
+
 [[deps.ArrayInterface]]
 deps = ["Adapt", "LinearAlgebra"]
 git-tree-sha1 = "78b3a7a536b4b0a747a0f296ea77091ca0a9f9a3"
@@ -151,10 +163,15 @@ uuid = "2a0f44e3-6c83-55bd-87e4-b1978d98bd5f"
 version = "1.11.0"
 
 [[deps.BenchmarkTools]]
-deps = ["Compat", "JSON", "Logging", "Printf", "Profile", "Statistics", "UUIDs"]
-git-tree-sha1 = "6876e30dc02dc69f0613cb6ece242144f2ca9e56"
+deps = ["Compat", "JSON", "Logging", "PrecompileTools", "Printf", "Profile", "Statistics", "UUIDs"]
+git-tree-sha1 = "9670d3febc2b6da60a0ae57846ba74670290653f"
 uuid = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"
-version = "1.7.0"
+version = "1.8.0"
+
+[[deps.Bijections]]
+git-tree-sha1 = "a2d308fcd4c2fb90e943cf9cd2fbfa9c32b69733"
+uuid = "e2ed5e7c-b2de-5872-ae92-c73ca462fb04"
+version = "0.2.2"
 
 [[deps.BitFlags]]
 git-tree-sha1 = "0691e34b3bb8be9307330f88d1a3c3f25466c24d"
@@ -210,9 +227,9 @@ version = "0.18.7"
 
 [[deps.CTDirect]]
 deps = ["ADNLPModels", "CTModels", "CTSolvers", "DocStringExtensions", "ExaModels", "SolverCore", "SparseArrays", "SparseConnectivityTracer"]
-git-tree-sha1 = "a4d812f60412f47bf29e685a921d0d675f581d55"
+git-tree-sha1 = "db1699ffab2f31aeda5dbb61dfa9015773cb70fc"
 uuid = "790bbbee-bee9-49ee-8912-a9de031322d5"
-version = "1.0.7-beta"
+version = "1.0.6"
 
 [[deps.CTFlows]]
 deps = ["CTBase", "CTModels", "DocStringExtensions", "ForwardDiff", "LinearAlgebra", "MLStyle", "MacroTools"]
@@ -248,9 +265,9 @@ version = "0.8.13"
 
 [[deps.CTSolvers]]
 deps = ["ADNLPModels", "CTBase", "CTModels", "CommonSolve", "DocStringExtensions", "ExaModels", "KernelAbstractions", "NLPModels", "SolverCore"]
-git-tree-sha1 = "37a9c44ca8fcc6948854fb796b75e545a4efbb9e"
+git-tree-sha1 = "22c283a24bd1b51cf2795074a6ec3fe0bff78adb"
 uuid = "d3e8d392-8e4b-4d9b-8e92-d7d4e3650ef6"
-version = "0.4.14"
+version = "0.4.15"
 
     [deps.CTSolvers.extensions]
     CTSolversCUDA = "CUDA"
@@ -334,6 +351,11 @@ git-tree-sha1 = "37ea44092930b1811e666c3bc38065d7d87fcc74"
 uuid = "5ae59095-9a9b-59fe-a467-6f913c188581"
 version = "0.13.1"
 
+[[deps.Combinatorics]]
+git-tree-sha1 = "08c8b6831dc00bfea825826be0bc8336fc369860"
+uuid = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
+version = "1.0.2"
+
 [[deps.CommonSolve]]
 git-tree-sha1 = "78ea4ddbcf9c241827e7035c3a03e2e456711470"
 uuid = "38540f10-b2f7-11e9-35d8-d573e4eb0ff2"
@@ -365,6 +387,11 @@ deps = ["Artifacts", "Libdl"]
 uuid = "e66e0078-7015-5450-92f7-15fbd957f2ae"
 version = "1.3.0+1"
 
+[[deps.CompositeTypes]]
+git-tree-sha1 = "bce26c3dab336582805503bed209faab1c279768"
+uuid = "b152e2b5-7a66-4b01-a709-34e65c35f657"
+version = "0.1.4"
+
 [[deps.CompositionsBase]]
 git-tree-sha1 = "802bb88cd69dfd1509f6670416bd4434015693ad"
 uuid = "a33af91c-f02d-484b-be07-31d278c5ca2b"
@@ -389,17 +416,13 @@ version = "2.5.1"
 git-tree-sha1 = "b4b092499347b18a015186eae3042f72267106cb"
 uuid = "187b0558-2788-49d3-abe0-74a17ed4e7c9"
 version = "1.6.0"
+weakdeps = ["IntervalSets", "LinearAlgebra", "StaticArrays"]
 
     [deps.ConstructionBase.extensions]
     ConstructionBaseIntervalSetsExt = "IntervalSets"
     ConstructionBaseLinearAlgebraExt = "LinearAlgebra"
     ConstructionBaseStaticArraysExt = "StaticArrays"
 
-    [deps.ConstructionBase.weakdeps]
-    IntervalSets = "8197267c-284f-5f27-9208-e0e47529a953"
-    LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
-    StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-
 [[deps.Contour]]
 git-tree-sha1 = "439e35b0b36e2e5881738abc8857bd92ad6ff9a8"
 uuid = "d38c429a-6771-53c6-b99e-75d170b6e991"
@@ -455,6 +478,12 @@ git-tree-sha1 = "9e2f36d3c96a820c678f2f1f1782582fcf685bae"
 uuid = "8bb1440f-4735-579b-a4ab-409b98df4dab"
 version = "1.9.1"
 
+[[deps.Dictionaries]]
+deps = ["Indexing", "Random", "Serialization"]
+git-tree-sha1 = "a55766a9c8f66cf19ffcdbdb1444e249bb4ace33"
+uuid = "85a47980-9c8c-11e8-2b9f-f7ca1fa99fb4"
+version = "0.4.6"
+
 [[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 = "87e2ad6d4ae98505218e2f97cafcfa296dc97d37"
@@ -588,11 +617,31 @@ git-tree-sha1 = "d8a8cb2d5b0181fbbd41861016b221b0202c9bc5"
 uuid = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
 version = "1.1.0"
 
+[[deps.DomainSets]]
+deps = ["CompositeTypes", "FunctionMaps", "IntervalSets", "LinearAlgebra", "StaticArrays"]
+git-tree-sha1 = "1af39efbaf76fb648432b5efaac0d73af6760407"
+uuid = "5b8099bc-c8ec-5219-889f-1d9e522a28bf"
+version = "0.8.0"
+
+    [deps.DomainSets.extensions]
+    DomainSetsMakieExt = "Makie"
+    DomainSetsRandomExt = "Random"
+
+    [deps.DomainSets.weakdeps]
+    Makie = "ee78f7c6-11fb-53f2-987a-cfe4a2b5a57a"
+    Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+
 [[deps.Downloads]]
 deps = ["ArgTools", "FileWatching", "LibCURL", "NetworkOptions"]
 uuid = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
 version = "1.7.0"
 
+[[deps.DynamicPolynomials]]
+deps = ["LinearAlgebra", "MultivariatePolynomials", "MutableArithmetics", "Reexport", "StarAlgebras", "Test"]
+git-tree-sha1 = "5bfabc3827dfdd164359bad6800c115a81280c00"
+uuid = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
+version = "0.6.6"
+
 [[deps.EnumX]]
 git-tree-sha1 = "c49898e8438c828577f04b92fc9368c388ac783c"
 uuid = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
@@ -815,6 +864,12 @@ git-tree-sha1 = "7a214fdac5ed5f59a22c2d9a885a16da1c74bbc7"
 uuid = "559328eb-81f9-559d-9380-de523a88c83c"
 version = "1.0.17+0"
 
+[[deps.FunctionMaps]]
+deps = ["CompositeTypes", "LinearAlgebra", "StaticArrays"]
+git-tree-sha1 = "31bd99a57edf98990d1c21486032963955450e8d"
+uuid = "a85aefff-f8ca-4649-a888-c8e5398bc76c"
+version = "0.1.2"
+
 [[deps.FunctionWrappers]]
 git-tree-sha1 = "d62485945ce5ae9c0c48f124a84998d755bae00e"
 uuid = "069b7b12-0de2-55c6-9aab-29f3d0a68a2e"
@@ -918,6 +973,19 @@ git-tree-sha1 = "8a6dbda1fd736d60cc477d99f2e7a042acfa46e8"
 uuid = "3b182d85-2403-5c21-9c21-1e1f0cc25472"
 version = "1.3.15+0"
 
+[[deps.Graphs]]
+deps = ["ArnoldiMethod", "DataStructures", "Inflate", "LinearAlgebra", "Random", "SimpleTraits", "SparseArrays", "Statistics"]
+git-tree-sha1 = "7eb45fe833a5b7c51cf6d89c5a841d5967e44be3"
+uuid = "86223c79-3864-5bf0-83f7-82e725a168b6"
+version = "1.14.0"
+
+    [deps.Graphs.extensions]
+    GraphsSharedArraysExt = "SharedArrays"
+
+    [deps.Graphs.weakdeps]
+    Distributed = "8ba89e20-285c-5b6f-9357-94700520ee1b"
+    SharedArrays = "1a1011a3-84de-559e-8e89-a11a2f7dc383"
+
 [[deps.Grisu]]
 git-tree-sha1 = "53bb909d1151e57e2484c3d1b53e19552b887fb2"
 uuid = "42e2da0e-8278-4e71-bc24-59509adca0fe"
@@ -952,6 +1020,16 @@ git-tree-sha1 = "debdd00ffef04665ccbb3e150747a77560e8fad1"
 uuid = "615f187c-cbe4-4ef1-ba3b-2fcf58d6d173"
 version = "0.1.1"
 
+[[deps.Indexing]]
+git-tree-sha1 = "ce1566720fd6b19ff3411404d4b977acd4814f9f"
+uuid = "313cdc1a-70c2-5d6a-ae34-0150d3930a38"
+version = "1.1.1"
+
+[[deps.Inflate]]
+git-tree-sha1 = "d1b1b796e47d94588b3757fe84fbf65a5ec4a80d"
+uuid = "d25df0c9-e2be-5dd7-82c8-3ad0b3e990b9"
+version = "0.1.5"
+
 [[deps.InlineStrings]]
 git-tree-sha1 = "8f3d257792a522b4601c24a577954b0a8cd7334d"
 uuid = "842dd82b-1e85-43dc-bf29-5d0ee9dffc48"
@@ -965,6 +1043,11 @@ version = "1.4.5"
     ArrowTypes = "31f734f8-188a-4ce0-8406-c8a06bd891cd"
     Parsers = "69de0a69-1ddd-5017-9359-2bf0b02dc9f0"
 
+[[deps.IntegerMathUtils]]
+git-tree-sha1 = "4c1acff2dc6b6967e7e750633c50bc3b8d83e617"
+uuid = "18e54dd8-cb9d-406c-a71d-865a43cbb235"
+version = "0.1.3"
+
 [[deps.IntelOpenMP_jll]]
 deps = ["Artifacts", "JLLWrappers", "LazyArtifacts", "Libdl"]
 git-tree-sha1 = "ec1debd61c300961f98064cfb21287613ad7f303"
@@ -976,6 +1059,17 @@ deps = ["Markdown"]
 uuid = "b77e0a4c-d291-57a0-90e8-8db25a27a240"
 version = "1.11.0"
 
+[[deps.IntervalSets]]
+git-tree-sha1 = "79d6bd28c8d9bccc2229784f1bd637689b256377"
+uuid = "8197267c-284f-5f27-9208-e0e47529a953"
+version = "0.7.14"
+weakdeps = ["Random", "RecipesBase", "Statistics"]
+
+    [deps.IntervalSets.extensions]
+    IntervalSetsRandomExt = "Random"
+    IntervalSetsRecipesBaseExt = "RecipesBase"
+    IntervalSetsStatisticsExt = "Statistics"
+
 [[deps.InverseFunctions]]
 git-tree-sha1 = "a779299d77cd080bf77b97535acecd73e1c5e5cb"
 uuid = "3587e190-3f89-42d0-90ee-14403ec27112"
@@ -1218,9 +1312,9 @@ weakdeps = ["LineSearches"]
 
 [[deps.LineSearches]]
 deps = ["LinearAlgebra", "NLSolversBase", "NaNMath", "Printf"]
-git-tree-sha1 = "738bdcacfef25b3a9e4a39c28613717a6b23751e"
+git-tree-sha1 = "a666999510c794fe1d9dfa629ef33366f11103aa"
 uuid = "d3d80556-e9d4-5f37-9878-2ab0fcc64255"
-version = "7.6.0"
+version = "7.6.1"
 
 [[deps.LinearAlgebra]]
 deps = ["Libdl", "OpenBLAS_jll", "libblastrampoline_jll"]
@@ -1460,6 +1554,22 @@ git-tree-sha1 = "cac9cc5499c25554cba55cd3c30543cff5ca4fab"
 uuid = "46d2c3a1-f734-5fdb-9937-b9b9aeba4221"
 version = "0.2.4"
 
+[[deps.MultivariatePolynomials]]
+deps = ["DataStructures", "LinearAlgebra", "MutableArithmetics", "StarAlgebras"]
+git-tree-sha1 = "b45f1ed8448ea20885cb4c5029c2a462fe2682bf"
+uuid = "102ac46a-7ee4-5c85-9060-abc95bfdeaa3"
+version = "0.5.17"
+weakdeps = ["ChainRulesCore"]
+
+    [deps.MultivariatePolynomials.extensions]
+    MultivariatePolynomialsChainRulesCoreExt = "ChainRulesCore"
+
+[[deps.MutableArithmetics]]
+deps = ["LinearAlgebra", "SparseArrays", "Test"]
+git-tree-sha1 = "7c25249fc13a070f5ba433c50e21e22bb33c6fb0"
+uuid = "d8a4904e-b15c-11e9-3269-09a3773c0cb0"
+version = "1.7.1"
+
 [[deps.NLPModels]]
 deps = ["FastClosures", "LinearAlgebra", "LinearOperators", "Printf", "SparseArrays"]
 git-tree-sha1 = "5919b1be82ca6ebfc9754cb22ddfb54f56f1d9ac"
@@ -1969,6 +2079,12 @@ version = "3.3.2"
     [deps.PrettyTables.weakdeps]
     Typstry = "f0ed7684-a786-439e-b1e3-3b82803b501e"
 
+[[deps.Primes]]
+deps = ["IntegerMathUtils"]
+git-tree-sha1 = "25cdd1d20cd005b52fc12cb6be3f75faaf59bb9b"
+uuid = "27ebfcd6-29c5-5fa9-bf4b-fb8fc14df3ae"
+version = "0.5.7"
+
 [[deps.Printf]]
 deps = ["Unicode"]
 uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
@@ -2024,6 +2140,11 @@ deps = ["SHA"]
 uuid = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 version = "1.11.0"
 
+[[deps.ReadOnlyArrays]]
+git-tree-sha1 = "e6f7ddf48cf141cb312b078ca21cb2d29d0dc11d"
+uuid = "988b38a3-91fc-5605-94a2-ee2116b3bd83"
+version = "0.2.0"
+
 [[deps.RecipesBase]]
 deps = ["PrecompileTools"]
 git-tree-sha1 = "5c3d09cc4f31f5fc6af001c250bf1278733100ff"
@@ -2256,6 +2377,12 @@ version = "2.11.1"
     ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
     Tracker = "9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c"
 
+[[deps.SimpleTraits]]
+deps = ["InteractiveUtils", "MacroTools"]
+git-tree-sha1 = "be8eeac05ec97d379347584fa9fe2f5f76795bcb"
+uuid = "699a6c99-e7fa-54fc-8d76-47d257e15c1d"
+version = "0.9.5"
+
 [[deps.Sockets]]
 uuid = "6462fe0b-24de-5631-8697-dd941f90decc"
 version = "1.11.0"
@@ -2336,6 +2463,12 @@ git-tree-sha1 = "4f96c596b8c8258cc7d3b19797854d368f243ddc"
 uuid = "860ef19b-820b-49d6-a774-d7a799459cd3"
 version = "1.0.4"
 
+[[deps.StarAlgebras]]
+deps = ["LinearAlgebra", "MutableArithmetics", "SparseArrays"]
+git-tree-sha1 = "235b1f9d287bbf34083b3d0829343a7942c0ad1c"
+uuid = "0c0c59c1-dc5f-42e9-9a8b-b5dc384a6cd1"
+version = "0.3.0"
+
 [[deps.Static]]
 deps = ["CommonWorldInvalidations", "IfElse", "PrecompileTools", "SciMLPublic"]
 git-tree-sha1 = "49440414711eddc7227724ae6e570c7d5559a086"
@@ -2408,9 +2541,9 @@ version = "0.4.4"
 
 [[deps.StructUtils]]
 deps = ["Dates", "UUIDs"]
-git-tree-sha1 = "fa95b3b097bcef5845c142ea2e085f1b2591e92c"
+git-tree-sha1 = "aab80fbf866600f3299dd7f6656d80e7be177cfe"
 uuid = "ec057cc2-7a8d-4b58-b3b3-92acb9f63b42"
-version = "2.7.1"
+version = "2.7.2"
 
     [deps.StructUtils.extensions]
     StructUtilsMeasurementsExt = ["Measurements"]
@@ -2451,6 +2584,63 @@ weakdeps = ["PrettyTables"]
     [deps.SymbolicIndexingInterface.extensions]
     SymbolicIndexingInterfacePrettyTablesExt = "PrettyTables"
 
+[[deps.SymbolicLimits]]
+deps = ["SymbolicUtils", "TermInterface"]
+git-tree-sha1 = "5085671d2cba1eb02136a3d6661c583e801984c1"
+uuid = "19f23fe9-fdab-4a78-91af-e7b7767979c3"
+version = "1.1.0"
+
+[[deps.SymbolicUtils]]
+deps = ["AbstractTrees", "ArrayInterface", "Combinatorics", "ConstructionBase", "DataStructures", "Dictionaries", "DocStringExtensions", "DynamicPolynomials", "EnumX", "ExproniconLite", "Graphs", "LinearAlgebra", "MacroTools", "Moshi", "MultivariatePolynomials", "MutableArithmetics", "NaNMath", "PrecompileTools", "ReadOnlyArrays", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArraysCore", "SymbolicIndexingInterface", "TaskLocalValues", "TermInterface", "WeakCacheSets"]
+git-tree-sha1 = "552eaeac659f5802d45c4d936272a712ca13a969"
+uuid = "d1185830-fcd6-423d-90d6-eec64667417b"
+version = "4.24.1"
+
+    [deps.SymbolicUtils.extensions]
+    SymbolicUtilsChainRulesCoreExt = "ChainRulesCore"
+    SymbolicUtilsDistributionsExt = "Distributions"
+    SymbolicUtilsLabelledArraysExt = "LabelledArrays"
+    SymbolicUtilsReverseDiffExt = "ReverseDiff"
+
+    [deps.SymbolicUtils.weakdeps]
+    ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"
+    Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+    LabelledArrays = "2ee39098-c373-598a-b85f-a56591580800"
+    ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267"
+
+[[deps.Symbolics]]
+deps = ["ADTypes", "AbstractPlutoDingetjes", "ArrayInterface", "Bijections", "CommonWorldInvalidations", "ConstructionBase", "DataStructures", "DiffRules", "DocStringExtensions", "DomainSets", "DynamicPolynomials", "Libdl", "LinearAlgebra", "LogExpFunctions", "MacroTools", "Markdown", "Moshi", "MultivariatePolynomials", "MutableArithmetics", "NaNMath", "PrecompileTools", "Preferences", "Primes", "RecipesBase", "Reexport", "RuntimeGeneratedFunctions", "SciMLPublic", "Setfield", "SparseArrays", "SpecialFunctions", "StaticArraysCore", "SymbolicIndexingInterface", "SymbolicLimits", "SymbolicUtils", "TermInterface"]
+git-tree-sha1 = "5c5db34512d5349b5e68d75cd707a3190fb26c81"
+uuid = "0c5d862f-8b57-4792-8d23-62f2024744c7"
+version = "7.18.1"
+
+    [deps.Symbolics.extensions]
+    SymbolicsD3TreesExt = "D3Trees"
+    SymbolicsDistributionsExt = "Distributions"
+    SymbolicsForwardDiffExt = "ForwardDiff"
+    SymbolicsGroebnerExt = "Groebner"
+    SymbolicsHypergeometricFunctionsExt = "HypergeometricFunctions"
+    SymbolicsLatexifyExt = ["Latexify", "LaTeXStrings"]
+    SymbolicsLuxExt = "Lux"
+    SymbolicsNemoExt = "Nemo"
+    SymbolicsPreallocationToolsExt = ["PreallocationTools", "ForwardDiff"]
+    SymbolicsSymPyExt = "SymPy"
+    SymbolicsSymPyPythonCallExt = "SymPyPythonCall"
+
+    [deps.Symbolics.weakdeps]
+    D3Trees = "e3df1716-f71e-5df9-9e2d-98e193103c45"
+    Distributions = "31c24e10-a181-5473-b8eb-7969acd0382f"
+    ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
+    Groebner = "0b43b601-686d-58a3-8a1c-6623616c7cd4"
+    HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
+    LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
+    Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
+    Lux = "b2108857-7c20-44ae-9111-449ecde12c47"
+    Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
+    PreallocationTools = "d236fae5-4411-538c-8e31-a6e3d9e00b46"
+    SymPy = "24249f21-da20-56a4-8eb1-6a02cf4ae2e6"
+    SymPyPythonCall = "bc8888f7-b21e-4b7c-a06a-5d9c9496438c"
+
 [[deps.TOML]]
 deps = ["Dates"]
 uuid = "fa267f1f-6049-4f14-aa54-33bafae1ed76"
@@ -2473,12 +2663,22 @@ deps = ["ArgTools", "SHA"]
 uuid = "a4e569a6-e804-4fa4-b0f3-eef7a1d5b13e"
 version = "1.10.0"
 
+[[deps.TaskLocalValues]]
+git-tree-sha1 = "67e469338d9ce74fc578f7db1736a74d93a49eb8"
+uuid = "ed4db957-447d-4319-bfb6-7fa9ae7ecf34"
+version = "0.1.3"
+
 [[deps.TensorCore]]
 deps = ["LinearAlgebra"]
 git-tree-sha1 = "1feb45f88d133a655e001435632f019a9a1bcdb6"
 uuid = "62fd8b95-f654-4bbd-a8a5-9c27f68ccd50"
 version = "0.1.1"
 
+[[deps.TermInterface]]
+git-tree-sha1 = "d673e0aca9e46a2f63720201f55cc7b3e7169b16"
+uuid = "8ea1fca8-c5ef-4a55-8b96-4e9afe9c9a3c"
+version = "2.0.0"
+
 [[deps.Test]]
 deps = ["InteractiveUtils", "Logging", "Random", "Serialization"]
 uuid = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
@@ -2566,6 +2766,11 @@ git-tree-sha1 = "96478df35bbc2f3e1e791bc7a3d0eeee559e60e9"
 uuid = "a2964d1f-97da-50d4-b82a-358c7fce9d89"
 version = "1.24.0+0"
 
+[[deps.WeakCacheSets]]
+git-tree-sha1 = "386050ae4353310d8ff9c228f83b1affca2f7f38"
+uuid = "d30d5f5c-d141-4870-aa07-aabb0f5fe7d5"
+version = "0.1.0"
+
 [[deps.XML2_jll]]
 deps = ["Artifacts", "JLLWrappers", "Libdl", "Libiconv_jll", "Zlib_jll"]
 git-tree-sha1 = "80d3930c6347cfce7ccf96bd3bafdf079d9c0390"
diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml
index 7966b8a..0d45e28 100644
--- a/docs/src/assets/Project.toml
+++ b/docs/src/assets/Project.toml
@@ -16,6 +16,7 @@ OrdinaryDiffEq = "1dea7af3-3e70-54e6-95c3-0bf5283fa5ed"
 Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb"
+Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
 BenchmarkTools = "1"
@@ -35,4 +36,5 @@ OrdinaryDiffEq = "6"
 Plots = "1"
 Printf = "1"
 Suppressor = "0.2"
+Symbolics = "7"
 julia = "1.10"
diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
new file mode 100644
index 0000000..5d0c72a
--- /dev/null
+++ b/docs/src/tutorial-symbolics.md
@@ -0,0 +1,6 @@
+# Titre
+
+```@meta
+Draft = false
+```
+

From 45ec95f5e36d42b82a5a10b7c6f8bb40ab47f55e Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 12:01:20 +0200
Subject: [PATCH 02/21] fix trigger

---
 .github/workflows/CI.yml            | 2 +-
 .github/workflows/Documentation.yml | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index 4f8bd47..c8b19ff 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -9,7 +9,7 @@ on:
     
 jobs:
   call:
-    if: ${{ ! contains(github.event.pull_request.labels.*.name, 'run ci') }}
+    if: contains(github.event.pull_request.labels.*.name, 'run ci')
     uses: control-toolbox/CTActions/.github/workflows/ci.yml@main
     with:
       runs_on: '["ubuntu-latest"]'
diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
index bdbd632..9bffc84 100644
--- a/.github/workflows/Documentation.yml
+++ b/.github/workflows/Documentation.yml
@@ -9,5 +9,5 @@ on:
   
 jobs:
   call:
-    if: ${{ ! contains(github.event.pull_request.labels.*.name, 'run documentation') }}
+    if: contains(github.event.pull_request.labels.*.name, 'run documentation')
     uses: control-toolbox/CTActions/.github/workflows/documentation.yml@main

From e1d4189bea0d2763158b19e2a893ac124690ee2e Mon Sep 17 00:00:00 2001
From: abavoil <17646791+abavoil@users.noreply.github.com>
Date: Wed, 15 Apr 2026 12:02:35 +0200
Subject: [PATCH 03/21] add script and dependencies for symbolics

---
 Project.toml                   |   4 +
 docs/src/tutorial-symbolics.md | 129 +++++++++++++++++++++++++++++++++
 2 files changed, 133 insertions(+)

diff --git a/Project.toml b/Project.toml
index 26db0c8..9139201 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,5 +3,9 @@ uuid = "cb10daa6-a5e5-4c25-a171-ae181b8ea3c9"
 version = "0.4.1"
 authors = ["Olivier Cots <olivier.cots@toulouse-inp.fr>"]
 
+[deps]
+StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
+
 [compat]
+StaticArrays = "1.9.18"
 julia = "1.10"
diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index 5d0c72a..0099859 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -4,3 +4,132 @@
 Draft = false
 ```
 
+```@example main
+# ==========================================
+# Setup & Imports
+# ==========================================
+import Pkg
+Pkg.activate("Control.jl/ModelingToolkit.jl")
+
+using OptimalControl
+using Plots
+using StaticArrays
+
+import NLPModelsIpopt
+
+using Symbolics
+
+# ==========================================
+# 1. Variables and Parameters
+# ==========================================
+# Physical Constants
+const m_c_val = 5.0
+const m_p_val = 1.0
+const l_val = 2.0
+const g_val = 9.81
+const tf_val = 2.0
+
+# Symbolic Variables
+@variables t
+D = Differential(t)
+@variables m_c m_p l g u
+@variables x(t) θ(t)
+@variables v ω dv dω # Static variables to solve for accelerations
+
+q = [x, θ]
+
+# ==========================================
+# 2. Automated Kinematics & Jacobians
+# ==========================================
+p_c = [x, 0.0]
+p_p = [x + l * sin(θ), l * cos(θ)]
+F = [u, 0.0]
+
+T = 0.5 * m_c * sum(D.(p_c) .^ 2) + 0.5 * m_p * sum(D.(p_p) .^ 2)
+V = g * (m_p * p_p[2])
+P_non_conservative = transpose(D.(p_c)) * F
+
+# ==========================================
+# 3. Euler-Lagrange & Matrix Inversion
+# ==========================================
+L = T - V
+
+A = D.(Symbolics.gradient(L, D.(q)))
+B = Symbolics.gradient(L, q)
+Q = Symbolics.gradient(P_non_conservative, D.(q))
+
+# Euler-Lagrange: d/dt( dL/dq_dot ) - dL/dq = Forces
+el_eqs = expand_derivatives.(A - B - Q)
+
+# Freeze time derivatives into static algebraic variables
+sub_rules = Dict(D(x) => v, D(θ) => ω, D(D(x)) => dv, D(D(θ)) => dω)
+res = Symbolics.substitute.(el_eqs, (sub_rules,))
+
+# Extract Mass Matrix and Bias vector
+Mass = Symbolics.jacobian(res, [dv, dω])
+bias = Symbolics.substitute.(res, (Dict(dv => 0.0, dω => 0.0),))
+
+# Analytically invert the Mass Matrix (Symbolics handles 2x2 natively)
+accel = Mass \ (-bias)
+accel = Symbolics.simplify_fractions.(accel)
+
+# The fully explicit state derivatives: X_dot = [v, ω, v_dot, ω_dot]
+dx_dt = [v, ω, accel[1], accel[2]]
+
+# ==========================================
+# 4. Extract Explicit Dynamics into a Julia Function
+# ==========================================
+f_expr = build_function(dx_dt, [x, θ, v, ω], u, [m_c, m_p, l, g];
+    expression=Val{false}, force_SA=true)
+f_mtk = f_expr[1]
+
+const p_vals = [m_c_val, m_p_val, l_val, g_val]
+
+cartpole_dynamics(X, U) = f_mtk(X, U, p_vals)
+
+# ==========================================
+# 5. OptimalControl.jl dynamic optimization problem
+# ==========================================
+@def cartpole_ocp begin
+    t ∈ [0, tf_val], time
+    X = (x_ocp, θ_ocp, v_ocp, ω_ocp) ∈ R⁴, state
+    F_ctrl ∈ R, control
+
+    X(0) == [0, 0, 0, 0.1]
+    X(tf_val) - X(0) == [0, 0, 0, 0]
+
+    # Inject our explicitly inverted Symbolic dynamics!
+    Ẋ(t) == cartpole_dynamics(X(t), F_ctrl(t))
+
+    ∫(F_ctrl(t)^2) → min
+end
+
+# ==========================================
+# 6. Solvers and Collocation
+# ==========================================
+initial_guess = (state=[0.0, 0.0, 0.0, 0.1], control=0.0)
+
+@time sol = solve(cartpole_ocp; display=true, grid_size=100, init=initial_guess)
+
+# ==========================================
+# 7. Extracting Results
+# ==========================================
+println("--- Optimal Limit Cycle Found ---")
+println("Total Control Energy (∫F² dt): ", sol.objective)
+
+tsol = time_grid(sol)
+# tsol = range(sol.model.times.initial.time, sol.model.times.final.time, 1001)
+Xsol = state(sol).(tsol)
+Fsol = control(sol).(tsol)
+
+xsol = [X[1] for X in Xsol]
+θsol = [X[2] for X in Xsol]
+vsol = [X[3] for X in Xsol]
+ωsol = [X[4] for X in Xsol]
+
+qplot = plot(tsol, [xsol θsol], label=["x" "θ"], title="Configuration")
+dqplot = plot(tsol, [vsol ωsol], label=["v" "ω"], title="Velocities")
+uplot = plot(tsol, Fsol, label="u", title="Control", linetype=:steppost)
+
+plot(qplot, dqplot, uplot, layout=3, size=(800, 600))
+```

From 44c5c1e6a25fa3f9bf7a2860e078227637c21504 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 12:05:35 +0200
Subject: [PATCH 04/21] fix trigger

---
 .github/workflows/CI.yml            | 1 +
 .github/workflows/Documentation.yml | 1 +
 2 files changed, 2 insertions(+)

diff --git a/.github/workflows/CI.yml b/.github/workflows/CI.yml
index c8b19ff..e742fcc 100644
--- a/.github/workflows/CI.yml
+++ b/.github/workflows/CI.yml
@@ -6,6 +6,7 @@ on:
       - main
     tags: '*'
   pull_request:
+    types: [labeled, opened, synchronize, reopened]
     
 jobs:
   call:
diff --git a/.github/workflows/Documentation.yml b/.github/workflows/Documentation.yml
index 9bffc84..96b813c 100644
--- a/.github/workflows/Documentation.yml
+++ b/.github/workflows/Documentation.yml
@@ -6,6 +6,7 @@ on:
       - main
     tags: '*'
   pull_request:
+    types: [labeled, opened, synchronize, reopened]
   
 jobs:
   call:

From 209cddcc1133f2c1f6eb43eed3636012a1e23d3f Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 12:14:26 +0200
Subject: [PATCH 05/21] Add detailed explanations and structure to
 tutorial-symbolics.md

---
 docs/src/tutorial-symbolics.md | 176 ++++++++++++++++++++++++++++++++-
 1 file changed, 175 insertions(+), 1 deletion(-)

diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index 0099859..515278d 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -1,9 +1,127 @@
-# Titre
+# Cart-Pole Limit Cycle via Symbolic Lagrangian Mechanics
 
 ```@meta
 Draft = false
 ```
 
+This tutorial demonstrates how to combine **symbolic derivation of equations of motion** (via
+[Symbolics.jl](https://symbolics.juliasymbolics.org/)) with **direct optimal control**
+(via [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl/)) to find a
+**periodic orbit** (limit cycle) for a cart-pole system.
+
+The key idea is to let the computer do the hard mechanics: we write the Lagrangian in a few
+lines, and the Euler–Lagrange equations — including mass-matrix inversion — are derived
+automatically.
+
+## The Cart-Pole System
+
+The system consists of a cart of mass ``m_c`` sliding on a frictionless horizontal rail, with
+a rigid pendulum of mass ``m_p`` and length ``l`` attached to it. A horizontal force ``u``
+(the control input) acts on the cart.
+
+```@raw html
+<figure>
+  <img src="cartpole_sketch.svg" alt="Cart-pole diagram" width="350"/>
+  <figcaption>Fig. 1 — Cart-pole system. The angle θ is measured from the upright position.</figcaption>
+</figure>
+```
+
+The configuration vector is ``q = (x,\, \theta)^\top``, where ``x`` is the cart position and
+``\theta = 0`` corresponds to the **upright** (unstable) equilibrium of the pendulum.
+
+### Positions
+
+The Cartesian positions of the two bodies are:
+
+```math
+p_c = \begin{pmatrix} x \\ 0 \end{pmatrix}, \qquad
+p_p = \begin{pmatrix} x + l\sin\theta \\ l\cos\theta \end{pmatrix}.
+```
+
+### Lagrangian
+
+The kinetic and potential energies are:
+
+```math
+T = \tfrac{1}{2}m_c\,\|\dot{p}_c\|^2 + \tfrac{1}{2}m_p\,\|\dot{p}_p\|^2
+  = \tfrac{1}{2}(m_c+m_p)\dot{x}^2
+    + m_p l\,\dot{x}\dot{\theta}\cos\theta
+    + \tfrac{1}{2}m_p l^2\dot{\theta}^2,
+```
+
+```math
+V = m_p\,g\,l\cos\theta.
+```
+
+The Lagrangian is ``\mathcal{L} = T - V``, and the virtual work of the control force gives the
+generalised force vector ``Q = (u,\, 0)^\top``.
+
+### Euler–Lagrange Equations
+
+The equations of motion follow from:
+
+```math
+\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{q}_i}
+- \frac{\partial \mathcal{L}}{\partial q_i} = Q_i, \qquad i = 1,2.
+```
+
+They can be written in the standard **manipulator form**:
+
+```math
+M(q)\,\ddot{q} = -C(q,\dot{q}) + \tau(q, u),
+```
+
+where the symmetric positive-definite **mass matrix** is:
+
+```math
+M(q) =
+\begin{pmatrix}
+  m_c + m_p & m_p l \cos\theta \\
+  m_p l \cos\theta & m_p l^2
+\end{pmatrix},
+```
+
+and the right-hand side collects Coriolis/gravity terms and the control torque. Instead of
+deriving these by hand, we rely on Symbolics.jl to compute and **analytically invert**
+``M(q)`` for us.
+
+### State-Space Form
+
+Defining the state ``X = (x,\,\theta,\,\dot{x},\,\dot{\theta})^\top``, the equations of
+motion become the first-order system:
+
+```math
+\dot{X}(t) = f\!\left(X(t),\, u(t)\right) =
+\begin{pmatrix}
+  \dot{x} \\ \dot{\theta} \\ M^{-1}(q)\bigl(-C(q,\dot{q}) + \tau(q,u)\bigr)
+\end{pmatrix}.
+```
+
+## The Optimal Control Problem
+
+We look for a **limit cycle** of the nonlinear dynamics: a trajectory that returns exactly to
+its initial condition after a fixed period ``t_f``. The cost penalises the total control
+energy:
+
+```math
+\min_{u(\cdot)}\; \int_0^{t_f} u(t)^2\,\mathrm{d}t
+```
+
+```math
+\text{subject to} \quad \dot{X}(t) = f(X(t), u(t)), \quad t \in [0, t_f],
+```
+
+```math
+X(t_f) = X(0).
+```
+
+The periodicity constraint ``X(t_f) = X(0)`` — combined with a non-trivial initial condition
+— forces the solver to find an orbit rather than the trivial rest solution.
+
+## Implementation
+
+### Setup & Imports
+
 ```@example main
 # ==========================================
 # Setup & Imports
@@ -18,7 +136,14 @@ using StaticArrays
 import NLPModelsIpopt
 
 using Symbolics
+```
 
+### Physical Parameters and Symbolic Variables
+
+We declare all parameters both as numerical constants (for the final function evaluation) and
+as symbolic variables (for the Lagrangian computation).
+
+```@example main
 # ==========================================
 # 1. Variables and Parameters
 # ==========================================
@@ -37,7 +162,14 @@ D = Differential(t)
 @variables v ω dv dω # Static variables to solve for accelerations
 
 q = [x, θ]
+```
+
+### Automated Kinematics and Lagrangian
 
+We express the positions, kinetic energy, potential energy, and virtual work symbolically.
+The time derivatives ``\dot{p}_c``, ``\dot{p}_p`` are computed automatically by `D.(...)`.
+
+```@example main
 # ==========================================
 # 2. Automated Kinematics & Jacobians
 # ==========================================
@@ -48,7 +180,17 @@ F = [u, 0.0]
 T = 0.5 * m_c * sum(D.(p_c) .^ 2) + 0.5 * m_p * sum(D.(p_p) .^ 2)
 V = g * (m_p * p_p[2])
 P_non_conservative = transpose(D.(p_c)) * F
+```
+
+### Euler–Lagrange Equations and Mass-Matrix Inversion
 
+Starting from ``\mathcal{L} = T - V``, Symbolics.jl computes the three terms of the
+Euler–Lagrange equations and assembles the residual vector. Substituting static aliases
+``(v,\omega,\dot{v},\dot\omega)`` for the time derivatives makes it possible to identify
+the mass matrix ``M`` as the Jacobian of the residual with respect to the accelerations
+``(\dot{v}, \dot\omega)``. The system ``M\,a = -b`` is then solved analytically.
+
+```@example main
 # ==========================================
 # 3. Euler-Lagrange & Matrix Inversion
 # ==========================================
@@ -75,7 +217,15 @@ accel = Symbolics.simplify_fractions.(accel)
 
 # The fully explicit state derivatives: X_dot = [v, ω, v_dot, ω_dot]
 dx_dt = [v, ω, accel[1], accel[2]]
+```
+
+### Code Generation
+
+`build_function` compiles the symbolic expression `dx_dt` into a native Julia function.
+The `force_SA=true` flag generates a **StaticArrays** kernel, which avoids heap allocations
+inside the ODE right-hand side — important for solver performance.
 
+```@example main
 # ==========================================
 # 4. Extract Explicit Dynamics into a Julia Function
 # ==========================================
@@ -86,7 +236,15 @@ f_mtk = f_expr[1]
 const p_vals = [m_c_val, m_p_val, l_val, g_val]
 
 cartpole_dynamics(X, U) = f_mtk(X, U, p_vals)
+```
+
+### Optimal Control Problem Definition
+
+We now formulate the optimal control problem using the `@def` macro from OptimalControl.jl.
+The periodicity condition is encoded as ``X(t_f) - X(0) = 0``, and the dynamics plug
+directly into the `Ẋ(t) == ...` constraint.
 
+```@example main
 # ==========================================
 # 5. OptimalControl.jl dynamic optimization problem
 # ==========================================
@@ -103,14 +261,25 @@ cartpole_dynamics(X, U) = f_mtk(X, U, p_vals)
 
     ∫(F_ctrl(t)^2) → min
 end
+```
+
+### Solving the NLP
 
+The problem is transcribed into a nonlinear program using direct collocation on a uniform
+grid of 100 intervals, then handed to **Ipopt** via NLPModelsIpopt.
+
+```@example main
 # ==========================================
 # 6. Solvers and Collocation
 # ==========================================
 initial_guess = (state=[0.0, 0.0, 0.0, 0.1], control=0.0)
 
 @time sol = solve(cartpole_ocp; display=true, grid_size=100, init=initial_guess)
+```
+
+### Results
 
+```@example main
 # ==========================================
 # 7. Extracting Results
 # ==========================================
@@ -133,3 +302,8 @@ uplot = plot(tsol, Fsol, label="u", title="Control", linetype=:steppost)
 
 plot(qplot, dqplot, uplot, layout=3, size=(800, 600))
 ```
+
+The three panels show, from top to bottom: the cart position ``x`` and pendulum angle
+``\theta``, the corresponding velocities ``\dot{x}`` and ``\dot\theta``, and the optimal
+control force ``u``. The periodicity of the state trajectory confirms that a genuine limit
+cycle has been found.

From 10c5f1388a3070496556b9b22ad0ff58a85bb091 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 12:22:39 +0200
Subject: [PATCH 06/21] add v2

---
 docs/src/tutorial-symbolics-2.md | 290 +++++++++++++++++++++++++++++++
 1 file changed, 290 insertions(+)
 create mode 100644 docs/src/tutorial-symbolics-2.md

diff --git a/docs/src/tutorial-symbolics-2.md b/docs/src/tutorial-symbolics-2.md
new file mode 100644
index 0000000..74f4f3f
--- /dev/null
+++ b/docs/src/tutorial-symbolics-2.md
@@ -0,0 +1,290 @@
+# Cart-Pole Limit Cycle via Symbolic Lagrangian Mechanics
+
+```@meta
+Draft = false
+```
+
+This tutorial demonstrates how to combine **symbolic derivation of equations of motion** (via
+[Symbolics.jl](https://symbolics.juliasymbolics.org/)) with **direct optimal control**
+(via [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl/)) to find a
+**periodic orbit** (limit cycle) for a cart-pole system.
+
+The key idea is to let the computer do the hard mechanics: we write the Lagrangian in a few
+lines, and the Euler–Lagrange equations — including mass-matrix inversion — are derived
+automatically.
+
+## The Cart-Pole System
+
+The system consists of a cart of mass ``m_c`` sliding on a frictionless horizontal rail, with
+a rigid pendulum of mass ``m_p`` and length ``l`` attached to it. A horizontal force ``u``
+(the control input) acts on the cart.
+
+```@raw html
+<figure>
+  <img src="cartpole_sketch.svg" alt="Cart-pole diagram" width="350"/>
+  <figcaption>Fig. 1 — Cart-pole system. The angle θ is measured from the upright position.</figcaption>
+</figure>
+```
+
+The configuration vector is ``q = (x,\, \theta)^\top``, where ``x`` is the cart position and
+``\theta = 0`` corresponds to the **upright** (unstable) equilibrium of the pendulum.
+
+### Positions
+
+The Cartesian positions of the two bodies are:
+
+```math
+p_c = \begin{pmatrix} x \\ 0 \end{pmatrix}, \qquad
+p_p = \begin{pmatrix} x + l\sin\theta \\ l\cos\theta \end{pmatrix}.
+```
+
+### Lagrangian
+
+The kinetic and potential energies are:
+
+```math
+T = \tfrac{1}{2}m_c\,\|\dot{p}_c\|^2 + \tfrac{1}{2}m_p\,\|\dot{p}_p\|^2
+  = \tfrac{1}{2}(m_c+m_p)\dot{x}^2
+    + m_p l\,\dot{x}\dot{\theta}\cos\theta
+    + \tfrac{1}{2}m_p l^2\dot{\theta}^2,
+```
+
+```math
+V = m_p\,g\,l\cos\theta.
+```
+
+The Lagrangian is ``\mathcal{L} = T - V``. The virtual work of the control force ``u``
+acting on the cart gives the generalised force vector:
+
+```math
+W_\text{ctrl} = F \cdot \dot{p}_c = u\,\dot{x},
+\quad\Longrightarrow\quad
+Q = \frac{\partial W_\text{ctrl}}{\partial \dot{q}} = \begin{pmatrix} u \\ 0 \end{pmatrix}.
+```
+
+### Euler–Lagrange Equations
+
+The equations of motion follow from:
+
+```math
+\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{q}_i}
+- \frac{\partial \mathcal{L}}{\partial q_i} = Q_i, \qquad i = 1,2.
+```
+
+They can be written in the standard **manipulator form**:
+
+```math
+M(q)\,\ddot{q} = -C(q,\dot{q}) + \tau(q, u),
+```
+
+where the symmetric positive-definite **mass matrix** is:
+
+```math
+M(q) =
+\begin{pmatrix}
+  m_c + m_p & m_p l \cos\theta \\
+  m_p l \cos\theta & m_p l^2
+\end{pmatrix},
+```
+
+and the right-hand side collects Coriolis/gravity terms and the control torque. Instead of
+deriving these by hand, we rely on Symbolics.jl to compute and **analytically invert**
+``M(q)`` for us.
+
+### State-Space Form
+
+Defining the state ``X = (x,\,\theta,\,\dot{x},\,\dot{\theta})^\top``, the equations of
+motion become the first-order system:
+
+```math
+\dot{X}(t) = f\!\left(X(t),\, u(t)\right) =
+\begin{pmatrix}
+  \dot{x} \\ \dot{\theta} \\ M^{-1}(q)\bigl(-C(q,\dot{q}) + \tau(q,u)\bigr)
+\end{pmatrix}.
+```
+
+## The Optimal Control Problem
+
+We look for a **limit cycle** of the nonlinear dynamics: a trajectory that returns exactly to
+its initial condition after a fixed period ``t_f``. The cost penalises the total control
+energy:
+
+```math
+\min_{u(\cdot)}\; \int_0^{t_f} u(t)^2\,\mathrm{d}t
+```
+
+```math
+\text{subject to} \quad \dot{X}(t) = f(X(t), u(t)), \quad t \in [0, t_f],
+```
+
+```math
+X(0) = X(t_f) = (0,\, 0,\, 0,\, 0.1)^\top.
+```
+
+The common boundary condition ``X(0) = X(t_f)`` — combined with a non-trivial initial
+angular velocity — forces the solver to find an orbit rather than the trivial rest solution.
+
+## Implementation
+
+### Setup & Imports
+
+```@setup main
+# Project environment for this tutorial — edit the path for your local setup.
+import Pkg
+Pkg.activate(".")
+```
+
+```@example main
+using OptimalControl
+using Plots
+using StaticArrays
+using NLPModelsIpopt
+using Symbolics
+```
+
+### Physical Parameters and Symbolic Variables
+
+We declare all parameters both as numerical constants (for the final function evaluation) and
+as symbolic variables (for the Lagrangian computation). Note that the symbolic time variable
+`t` is used only inside the Symbolics.jl derivation and does not interfere with the time
+variable `t` introduced later by the `@def` macro.
+
+```@example main
+# Physical constants
+const m_c_val = 5.0
+const m_p_val = 1.0
+const l_val   = 2.0
+const g_val   = 9.81
+const tf_val  = 2.0
+
+# Symbolic variables
+@variables t
+D = Differential(t)
+@variables m_c m_p l g u
+@variables x(t) θ(t)
+@variables v ω dv dω   # static aliases for velocities and accelerations
+
+q = [x, θ]
+```
+
+### Automated Kinematics and Lagrangian
+
+We express the positions, kinetic energy, potential energy, and virtual power symbolically.
+The time derivatives ``\dot{p}_c``, ``\dot{p}_p`` are computed automatically by `D.(...)`.
+
+```@example main
+p_c = [x, 0.0]
+p_p = [x + l * sin(θ), l * cos(θ)]
+F   = [u, 0.0]
+
+T      = 0.5 * m_c * sum(D.(p_c) .^ 2) + 0.5 * m_p * sum(D.(p_p) .^ 2)
+V      = g * (m_p * p_p[2])
+W_ctrl = transpose(D.(p_c)) * F   # virtual power of the control force
+```
+
+### Euler–Lagrange Equations and Mass-Matrix Inversion
+
+Starting from ``\mathcal{L} = T - V``, Symbolics.jl computes the three terms of the
+Euler–Lagrange equations and assembles the residual vector. Substituting static aliases
+``(v,\omega,\dot{v},\dot\omega)`` for the time derivatives makes it possible to identify
+the mass matrix ``M`` as the Jacobian of the residual with respect to the accelerations
+``(\dot{v}, \dot\omega)``. The system ``M\,a = -b`` is then solved analytically.
+
+```@example main
+L = T - V
+
+A = D.(Symbolics.gradient(L, D.(q)))
+B = Symbolics.gradient(L, q)
+Q = Symbolics.gradient(W_ctrl, D.(q))   # generalised forces
+
+# Euler-Lagrange residual: d/dt(∂L/∂q̇) - ∂L/∂q - Q = 0
+el_eqs = expand_derivatives.(A - B - Q)
+
+# Freeze time derivatives into static algebraic variables
+sub_rules = Dict(D(x) => v, D(θ) => ω, D(D(x)) => dv, D(D(θ)) => dω)
+res = Symbolics.substitute.(el_eqs, (sub_rules,))
+
+# Identify mass matrix M and bias vector b such that M·[dv; dω] + b = 0
+Mass = Symbolics.jacobian(res, [dv, dω])
+bias = Symbolics.substitute.(res, (Dict(dv => 0.0, dω => 0.0),))
+
+# Analytically invert the 2×2 mass matrix
+accel = Symbolics.simplify_fractions.(Mass \ (-bias))
+
+# Fully explicit state derivative: Ẋ = [v, ω, v̇, ω̇]
+dx_dt = [v, ω, accel[1], accel[2]]
+```
+
+### Code Generation
+
+`build_function` compiles the symbolic expression `dx_dt` into a native Julia function.
+The `force_SA=true` flag generates a **StaticArrays** kernel, which avoids heap allocations
+inside the ODE right-hand side — important for solver performance. Parameters are stored as
+an `SVector` to match the generated kernel's expected input type.
+
+```@example main
+f_expr = build_function(dx_dt, [x, θ, v, ω], u, [m_c, m_p, l, g];
+    expression=Val{false}, force_SA=true)
+f_cartpole = f_expr[1]   # out-of-place variant: (state, u, params) → SVector
+
+const p_vals = SA[m_c_val, m_p_val, l_val, g_val]   # SVector, matches SA kernel
+
+cartpole_dynamics(X, U) = f_cartpole(X, U, p_vals)
+```
+
+### Optimal Control Problem Definition
+
+We now formulate the optimal control problem using the `@def` macro from OptimalControl.jl.
+The boundary conditions ``X(0) = X(t_f)`` encode the periodicity of the orbit directly.
+
+```@example main
+@def cartpole_ocp begin
+    t ∈ [0, tf_val], time
+    X = (x_ocp, θ_ocp, v_ocp, ω_ocp) ∈ R⁴, state
+    F_ctrl ∈ R, control
+
+    X(0)      == [0, 0, 0, 0.1]
+    X(tf_val) == [0, 0, 0, 0.1]   # periodicity: X(tf) = X(0)
+
+    Ẋ(t) == cartpole_dynamics(X(t), F_ctrl(t))
+
+    ∫(F_ctrl(t)^2) → min
+end
+```
+
+### Solving the NLP
+
+The problem is transcribed into a nonlinear program using direct collocation on a uniform
+grid of 100 intervals, then handed to **Ipopt** via NLPModelsIpopt. The `@time` macro
+reports wall-clock time, which on first call includes Julia's JIT compilation.
+
+```@example main
+initial_guess = (state=[0.0, 0.0, 0.0, 0.1], control=0.0)
+
+@time sol = solve(cartpole_ocp; display=true, grid_size=100, init=initial_guess)
+```
+
+### Results
+
+```@example main
+println("--- Optimal Limit Cycle Found ---")
+println("Total Control Energy (∫F² dt): ", sol.objective)
+
+tsol = time_grid(sol)
+Xsol = state(sol).(tsol)
+Fsol = control(sol).(tsol)
+
+# Stack state vectors into a 4×N matrix and unpack rows
+xsol, θsol, vsol, ωsol = eachrow(reduce(hcat, Xsol))
+
+qplot  = plot(tsol, [xsol θsol], label=["x" "θ"],   title="Configuration")
+dqplot = plot(tsol, [vsol ωsol], label=["v" "ω"],   title="Velocities")
+uplot  = plot(tsol, Fsol,        label="u",          title="Control", linetype=:steppost)
+
+plot(qplot, dqplot, uplot, layout=3, size=(800, 600))
+```
+
+The three panels show the cart position ``x`` and pendulum angle ``\theta``, the
+corresponding velocities ``\dot{x}`` and ``\dot\theta``, and the optimal control force
+``u``. The periodicity of the state trajectory confirms that a genuine limit cycle has been
+found.
\ No newline at end of file

From 534f8ef4a842f807cadb70854dc7eea88c94b49d Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 14:15:48 +0200
Subject: [PATCH 07/21] Add interactive JS canvas animation to
 tutorial-symbolics-2

---
 docs/src/tutorial-symbolics-2.md | 112 ++++++++++++++++++++++++++++++-
 1 file changed, 111 insertions(+), 1 deletion(-)

diff --git a/docs/src/tutorial-symbolics-2.md b/docs/src/tutorial-symbolics-2.md
index 74f4f3f..891f958 100644
--- a/docs/src/tutorial-symbolics-2.md
+++ b/docs/src/tutorial-symbolics-2.md
@@ -287,4 +287,114 @@ plot(qplot, dqplot, uplot, layout=3, size=(800, 600))
 The three panels show the cart position ``x`` and pendulum angle ``\theta``, the
 corresponding velocities ``\dot{x}`` and ``\dot\theta``, and the optimal control force
 ``u``. The periodicity of the state trajectory confirms that a genuine limit cycle has been
-found.
\ No newline at end of file
+found.
+
+### Animation
+
+The animation below shows the cart-pole evolving along the optimal limit-cycle trajectory.
+The blue cart slides on the horizontal rail while the pendulum swings around the upright
+equilibrium. A ghost trace follows the red bob to make the motion easier to read. Use the
+controls to play, pause, or reset the animation and to adjust the playback speed.
+
+```@example main
+# Serialise trajectory for the JS animation (manual, no extra dependencies)
+json_t  = "[" * join(string.(tsol),  ",") * "]"
+json_x  = "[" * join(string.(xsol),  ",") * "]"
+json_th = "[" * join(string.(θsol),  ",") * "]"
+
+uid    = string(rand(UInt32), base=16)   # unique suffix → safe for multi-example pages
+tf_str = string(round(tsol[end], digits=3))
+
+html_str = """
+<div style="text-align:center;font-family:sans-serif;margin:1.5em 0;">
+  <canvas id="cpC$(uid)" width="700" height="380"
+    style="border:1px solid #ccc;border-radius:6px;background:#fafafa;display:block;margin:0 auto;"></canvas>
+  <div style="margin:6px 0;font-size:0.9em;color:#444;">
+    <span id="cpT$(uid)">t = 0.000 / $(tf_str) s</span>
+  </div>
+  <div style="display:flex;justify-content:center;gap:10px;align-items:center;flex-wrap:wrap;margin-top:4px;">
+    <button id="cpPl$(uid)" style="padding:5px 14px;cursor:pointer;">&#9654; Play</button>
+    <button id="cpPa$(uid)" style="padding:5px 14px;cursor:pointer;">&#9646;&#9646; Pause</button>
+    <button id="cpRe$(uid)" style="padding:5px 14px;cursor:pointer;">&#x21BA; Reset</button>
+    <label style="font-size:0.9em;">Speed:
+      <input type="range" id="cpSp$(uid)" min="0.25" max="3" step="0.25" value="1"
+             style="vertical-align:middle;width:100px;">
+      <span id="cpSv$(uid)">1.00&times;</span>
+    </label>
+  </div>
+</div>
+<script>
+(function(){
+  var T=$(json_t), X=$(json_x), TH=$(json_th), L=$(l_val);
+  var cv=document.getElementById("cpC$(uid)");
+  var ctx=cv.getContext("2d");
+  var W=cv.width, H=cv.height;
+  var railY=H*0.62, sc=80, cW=60, cH=28;
+  var tSim=0, spd=1, rid=null, lts=null, ghost=[];
+  function lerp(a,b,t){return a+(b-a)*t;}
+  function interp(tq){
+    var n=T.length;
+    if(tq<=T[0])return{x:X[0],th:TH[0]};
+    if(tq>=T[n-1])return{x:X[n-1],th:TH[n-1]};
+    var lo=0,hi=n-2;
+    while(lo<hi){var m=(lo+hi)>>1;if(T[m+1]<tq)lo=m+1;else hi=m;}
+    var a=(tq-T[lo])/(T[lo+1]-T[lo]);
+    return{x:lerp(X[lo],X[lo+1],a),th:lerp(TH[lo],TH[lo+1],a)};
+  }
+  function draw(tq){
+    var d=interp(tq);
+    var cx=W/2+d.x*sc, cy=railY;
+    var px=cx, py=cy-cH/2;
+    var bx=px+L*sc*Math.sin(d.th), by=py-L*sc*Math.cos(d.th);
+    ghost.push({x:bx,y:by}); if(ghost.length>25)ghost.shift();
+    ctx.clearRect(0,0,W,H);
+    ctx.beginPath();ctx.moveTo(20,railY);ctx.lineTo(W-20,railY);
+    ctx.strokeStyle="#aaa";ctx.lineWidth=3;ctx.stroke();
+    var wr=8;
+    [[cx-18,cy+cH/2+wr],[cx+18,cy+cH/2+wr]].forEach(function(p){
+      ctx.beginPath();ctx.arc(p[0],p[1],wr,0,2*Math.PI);
+      ctx.fillStyle="#888";ctx.fill();
+    });
+    ctx.beginPath();ctx.rect(cx-cW/2,cy-cH/2,cW,cH);
+    ctx.fillStyle="#4a90d9";ctx.fill();
+    ctx.strokeStyle="#2c5f8a";ctx.lineWidth=1.5;ctx.stroke();
+    ghost.forEach(function(g,i){
+      ctx.beginPath();ctx.arc(g.x,g.y,5,0,2*Math.PI);
+      ctx.fillStyle="rgba(220,60,60,"+((i+1)/ghost.length*0.35)+")";ctx.fill();
+    });
+    ctx.beginPath();ctx.moveTo(px,py);ctx.lineTo(bx,by);
+    ctx.strokeStyle="#222";ctx.lineWidth=3;ctx.stroke();
+    ctx.beginPath();ctx.arc(px,py,5,0,2*Math.PI);ctx.fillStyle="#555";ctx.fill();
+    ctx.beginPath();ctx.arc(bx,by,10,0,2*Math.PI);
+    ctx.fillStyle="#dc3c3c";ctx.fill();ctx.strokeStyle="#8a1a1a";ctx.lineWidth=1.5;ctx.stroke();
+    document.getElementById("cpT$(uid)").textContent=
+      "t = "+tq.toFixed(3)+" / "+T[T.length-1].toFixed(2)+" s";
+  }
+  function step(ts){
+    if(lts!==null){
+      var dt=(ts-lts)*0.001*spd; tSim+=dt;
+      var tf=T[T.length-1]; if(tSim>tf){tSim-=tf;ghost=[];}
+    }
+    lts=ts; draw(tSim); rid=requestAnimationFrame(step);
+  }
+  document.getElementById("cpPl$(uid)").onclick=function(){
+    if(rid===null){lts=null;rid=requestAnimationFrame(step);}
+  };
+  document.getElementById("cpPa$(uid)").onclick=function(){
+    if(rid!==null){cancelAnimationFrame(rid);rid=null;}
+  };
+  document.getElementById("cpRe$(uid)").onclick=function(){
+    if(rid!==null){cancelAnimationFrame(rid);rid=null;}
+    tSim=0;ghost=[];lts=null;draw(0);
+  };
+  document.getElementById("cpSp$(uid)").oninput=function(){
+    spd=parseFloat(this.value);
+    document.getElementById("cpSv$(uid)").textContent=spd.toFixed(2)+"\u00d7";
+  };
+  draw(0);
+})();
+</script>
+"""
+
+HTML(html_str)
+```
\ No newline at end of file

From 9af509b3542880ec1bff39f15888131917318dfd Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 14:31:06 +0200
Subject: [PATCH 08/21] review cart pole

---
 docs/src/assets/Manifest.toml       |   6 +-
 docs/src/assets/cartpole_sketch.svg |  39 +++
 docs/src/tutorial-symbolics-2.md    | 400 ----------------------------
 docs/src/tutorial-symbolics-save.md | 313 ++++++++++++++++++++++
 docs/src/tutorial-symbolics.md      | 247 +++++++++++------
 5 files changed, 526 insertions(+), 479 deletions(-)
 create mode 100644 docs/src/assets/cartpole_sketch.svg
 delete mode 100644 docs/src/tutorial-symbolics-2.md
 create mode 100644 docs/src/tutorial-symbolics-save.md

diff --git a/docs/src/assets/Manifest.toml b/docs/src/assets/Manifest.toml
index 65c19ad..ced7cf3 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.6"
+julia_version = "1.12.5"
 manifest_format = "2.0"
 project_hash = "1d73fc8b7306dbb42eb2a1fdd2ae0c8be50a6d19"
 
@@ -227,9 +227,9 @@ version = "0.18.7"
 
 [[deps.CTDirect]]
 deps = ["ADNLPModels", "CTModels", "CTSolvers", "DocStringExtensions", "ExaModels", "SolverCore", "SparseArrays", "SparseConnectivityTracer"]
-git-tree-sha1 = "db1699ffab2f31aeda5dbb61dfa9015773cb70fc"
+git-tree-sha1 = "a4d812f60412f47bf29e685a921d0d675f581d55"
 uuid = "790bbbee-bee9-49ee-8912-a9de031322d5"
-version = "1.0.6"
+version = "1.0.7-beta"
 
 [[deps.CTFlows]]
 deps = ["CTBase", "CTModels", "DocStringExtensions", "ForwardDiff", "LinearAlgebra", "MLStyle", "MacroTools"]
diff --git a/docs/src/assets/cartpole_sketch.svg b/docs/src/assets/cartpole_sketch.svg
new file mode 100644
index 0000000..ea4addb
--- /dev/null
+++ b/docs/src/assets/cartpole_sketch.svg
@@ -0,0 +1,39 @@
+<svg width="350" height="250" viewBox="0 0 350 250" xmlns="http://www.w3.org/2000/svg">
+  <!-- Background -->
+  <rect width="100%" height="100%" fill="white"/>
+  
+  <!-- Rail -->
+  <line x1="30" y1="180" x2="320" y2="180" stroke="#666" stroke-width="3"/>
+  
+  <!-- Cart wheels -->
+  <circle cx="130" cy="188" r="8" fill="#888"/>
+  <circle cx="170" cy="188" r="8" fill="#888"/>
+  
+  <!-- Cart body -->
+  <rect x="120" y="160" width="60" height="25" rx="3" fill="#4a90d9" stroke="#2c5f8a" stroke-width="2"/>
+  
+  <!-- Pendulum rod -->
+  <line x1="150" y1="160" x2="150" y2="80" stroke="#333" stroke-width="3"/>
+  
+  <!-- Pendulum bob -->
+  <circle cx="150" cy="80" r="12" fill="#dc3c3c" stroke="#8a1a1a" stroke-width="2"/>
+  
+  <!-- Angle θ annotation -->
+  <path d="M 150 160 L 150 140" stroke="#666" stroke-width="1" stroke-dasharray="3,3"/>
+  <text x="155" y="145" font-family="Arial, sans-serif" font-size="14" fill="#666">θ</text>
+  
+  <!-- Force u arrow -->
+  <line x1="180" y1="172" x2="220" y2="172" stroke="#e74c3c" stroke-width="2" marker-end="url(#arrowhead)"/>
+  <text x="200" y="165" font-family="Arial, sans-serif" font-size="12" fill="#e74c3c" text-anchor="middle">u</text>
+  
+  <!-- Arrow marker definition -->
+  <defs>
+    <marker id="arrowhead" markerWidth="10" markerHeight="7" refX="9" refY="3.5" orient="auto">
+      <polygon points="0 0, 10 3.5, 0 7" fill="#e74c3c"/>
+    </marker>
+  </defs>
+  
+  <!-- Labels -->
+  <text x="150" y="220" font-family="Arial, sans-serif" font-size="12" fill="#333" text-anchor="middle">x</text>
+  <text x="150" y="240" font-family="Arial, sans-serif" font-size="12" fill="#666" text-anchor="middle">rail</text>
+</svg>
diff --git a/docs/src/tutorial-symbolics-2.md b/docs/src/tutorial-symbolics-2.md
deleted file mode 100644
index 891f958..0000000
--- a/docs/src/tutorial-symbolics-2.md
+++ /dev/null
@@ -1,400 +0,0 @@
-# Cart-Pole Limit Cycle via Symbolic Lagrangian Mechanics
-
-```@meta
-Draft = false
-```
-
-This tutorial demonstrates how to combine **symbolic derivation of equations of motion** (via
-[Symbolics.jl](https://symbolics.juliasymbolics.org/)) with **direct optimal control**
-(via [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl/)) to find a
-**periodic orbit** (limit cycle) for a cart-pole system.
-
-The key idea is to let the computer do the hard mechanics: we write the Lagrangian in a few
-lines, and the Euler–Lagrange equations — including mass-matrix inversion — are derived
-automatically.
-
-## The Cart-Pole System
-
-The system consists of a cart of mass ``m_c`` sliding on a frictionless horizontal rail, with
-a rigid pendulum of mass ``m_p`` and length ``l`` attached to it. A horizontal force ``u``
-(the control input) acts on the cart.
-
-```@raw html
-<figure>
-  <img src="cartpole_sketch.svg" alt="Cart-pole diagram" width="350"/>
-  <figcaption>Fig. 1 — Cart-pole system. The angle θ is measured from the upright position.</figcaption>
-</figure>
-```
-
-The configuration vector is ``q = (x,\, \theta)^\top``, where ``x`` is the cart position and
-``\theta = 0`` corresponds to the **upright** (unstable) equilibrium of the pendulum.
-
-### Positions
-
-The Cartesian positions of the two bodies are:
-
-```math
-p_c = \begin{pmatrix} x \\ 0 \end{pmatrix}, \qquad
-p_p = \begin{pmatrix} x + l\sin\theta \\ l\cos\theta \end{pmatrix}.
-```
-
-### Lagrangian
-
-The kinetic and potential energies are:
-
-```math
-T = \tfrac{1}{2}m_c\,\|\dot{p}_c\|^2 + \tfrac{1}{2}m_p\,\|\dot{p}_p\|^2
-  = \tfrac{1}{2}(m_c+m_p)\dot{x}^2
-    + m_p l\,\dot{x}\dot{\theta}\cos\theta
-    + \tfrac{1}{2}m_p l^2\dot{\theta}^2,
-```
-
-```math
-V = m_p\,g\,l\cos\theta.
-```
-
-The Lagrangian is ``\mathcal{L} = T - V``. The virtual work of the control force ``u``
-acting on the cart gives the generalised force vector:
-
-```math
-W_\text{ctrl} = F \cdot \dot{p}_c = u\,\dot{x},
-\quad\Longrightarrow\quad
-Q = \frac{\partial W_\text{ctrl}}{\partial \dot{q}} = \begin{pmatrix} u \\ 0 \end{pmatrix}.
-```
-
-### Euler–Lagrange Equations
-
-The equations of motion follow from:
-
-```math
-\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{q}_i}
-- \frac{\partial \mathcal{L}}{\partial q_i} = Q_i, \qquad i = 1,2.
-```
-
-They can be written in the standard **manipulator form**:
-
-```math
-M(q)\,\ddot{q} = -C(q,\dot{q}) + \tau(q, u),
-```
-
-where the symmetric positive-definite **mass matrix** is:
-
-```math
-M(q) =
-\begin{pmatrix}
-  m_c + m_p & m_p l \cos\theta \\
-  m_p l \cos\theta & m_p l^2
-\end{pmatrix},
-```
-
-and the right-hand side collects Coriolis/gravity terms and the control torque. Instead of
-deriving these by hand, we rely on Symbolics.jl to compute and **analytically invert**
-``M(q)`` for us.
-
-### State-Space Form
-
-Defining the state ``X = (x,\,\theta,\,\dot{x},\,\dot{\theta})^\top``, the equations of
-motion become the first-order system:
-
-```math
-\dot{X}(t) = f\!\left(X(t),\, u(t)\right) =
-\begin{pmatrix}
-  \dot{x} \\ \dot{\theta} \\ M^{-1}(q)\bigl(-C(q,\dot{q}) + \tau(q,u)\bigr)
-\end{pmatrix}.
-```
-
-## The Optimal Control Problem
-
-We look for a **limit cycle** of the nonlinear dynamics: a trajectory that returns exactly to
-its initial condition after a fixed period ``t_f``. The cost penalises the total control
-energy:
-
-```math
-\min_{u(\cdot)}\; \int_0^{t_f} u(t)^2\,\mathrm{d}t
-```
-
-```math
-\text{subject to} \quad \dot{X}(t) = f(X(t), u(t)), \quad t \in [0, t_f],
-```
-
-```math
-X(0) = X(t_f) = (0,\, 0,\, 0,\, 0.1)^\top.
-```
-
-The common boundary condition ``X(0) = X(t_f)`` — combined with a non-trivial initial
-angular velocity — forces the solver to find an orbit rather than the trivial rest solution.
-
-## Implementation
-
-### Setup & Imports
-
-```@setup main
-# Project environment for this tutorial — edit the path for your local setup.
-import Pkg
-Pkg.activate(".")
-```
-
-```@example main
-using OptimalControl
-using Plots
-using StaticArrays
-using NLPModelsIpopt
-using Symbolics
-```
-
-### Physical Parameters and Symbolic Variables
-
-We declare all parameters both as numerical constants (for the final function evaluation) and
-as symbolic variables (for the Lagrangian computation). Note that the symbolic time variable
-`t` is used only inside the Symbolics.jl derivation and does not interfere with the time
-variable `t` introduced later by the `@def` macro.
-
-```@example main
-# Physical constants
-const m_c_val = 5.0
-const m_p_val = 1.0
-const l_val   = 2.0
-const g_val   = 9.81
-const tf_val  = 2.0
-
-# Symbolic variables
-@variables t
-D = Differential(t)
-@variables m_c m_p l g u
-@variables x(t) θ(t)
-@variables v ω dv dω   # static aliases for velocities and accelerations
-
-q = [x, θ]
-```
-
-### Automated Kinematics and Lagrangian
-
-We express the positions, kinetic energy, potential energy, and virtual power symbolically.
-The time derivatives ``\dot{p}_c``, ``\dot{p}_p`` are computed automatically by `D.(...)`.
-
-```@example main
-p_c = [x, 0.0]
-p_p = [x + l * sin(θ), l * cos(θ)]
-F   = [u, 0.0]
-
-T      = 0.5 * m_c * sum(D.(p_c) .^ 2) + 0.5 * m_p * sum(D.(p_p) .^ 2)
-V      = g * (m_p * p_p[2])
-W_ctrl = transpose(D.(p_c)) * F   # virtual power of the control force
-```
-
-### Euler–Lagrange Equations and Mass-Matrix Inversion
-
-Starting from ``\mathcal{L} = T - V``, Symbolics.jl computes the three terms of the
-Euler–Lagrange equations and assembles the residual vector. Substituting static aliases
-``(v,\omega,\dot{v},\dot\omega)`` for the time derivatives makes it possible to identify
-the mass matrix ``M`` as the Jacobian of the residual with respect to the accelerations
-``(\dot{v}, \dot\omega)``. The system ``M\,a = -b`` is then solved analytically.
-
-```@example main
-L = T - V
-
-A = D.(Symbolics.gradient(L, D.(q)))
-B = Symbolics.gradient(L, q)
-Q = Symbolics.gradient(W_ctrl, D.(q))   # generalised forces
-
-# Euler-Lagrange residual: d/dt(∂L/∂q̇) - ∂L/∂q - Q = 0
-el_eqs = expand_derivatives.(A - B - Q)
-
-# Freeze time derivatives into static algebraic variables
-sub_rules = Dict(D(x) => v, D(θ) => ω, D(D(x)) => dv, D(D(θ)) => dω)
-res = Symbolics.substitute.(el_eqs, (sub_rules,))
-
-# Identify mass matrix M and bias vector b such that M·[dv; dω] + b = 0
-Mass = Symbolics.jacobian(res, [dv, dω])
-bias = Symbolics.substitute.(res, (Dict(dv => 0.0, dω => 0.0),))
-
-# Analytically invert the 2×2 mass matrix
-accel = Symbolics.simplify_fractions.(Mass \ (-bias))
-
-# Fully explicit state derivative: Ẋ = [v, ω, v̇, ω̇]
-dx_dt = [v, ω, accel[1], accel[2]]
-```
-
-### Code Generation
-
-`build_function` compiles the symbolic expression `dx_dt` into a native Julia function.
-The `force_SA=true` flag generates a **StaticArrays** kernel, which avoids heap allocations
-inside the ODE right-hand side — important for solver performance. Parameters are stored as
-an `SVector` to match the generated kernel's expected input type.
-
-```@example main
-f_expr = build_function(dx_dt, [x, θ, v, ω], u, [m_c, m_p, l, g];
-    expression=Val{false}, force_SA=true)
-f_cartpole = f_expr[1]   # out-of-place variant: (state, u, params) → SVector
-
-const p_vals = SA[m_c_val, m_p_val, l_val, g_val]   # SVector, matches SA kernel
-
-cartpole_dynamics(X, U) = f_cartpole(X, U, p_vals)
-```
-
-### Optimal Control Problem Definition
-
-We now formulate the optimal control problem using the `@def` macro from OptimalControl.jl.
-The boundary conditions ``X(0) = X(t_f)`` encode the periodicity of the orbit directly.
-
-```@example main
-@def cartpole_ocp begin
-    t ∈ [0, tf_val], time
-    X = (x_ocp, θ_ocp, v_ocp, ω_ocp) ∈ R⁴, state
-    F_ctrl ∈ R, control
-
-    X(0)      == [0, 0, 0, 0.1]
-    X(tf_val) == [0, 0, 0, 0.1]   # periodicity: X(tf) = X(0)
-
-    Ẋ(t) == cartpole_dynamics(X(t), F_ctrl(t))
-
-    ∫(F_ctrl(t)^2) → min
-end
-```
-
-### Solving the NLP
-
-The problem is transcribed into a nonlinear program using direct collocation on a uniform
-grid of 100 intervals, then handed to **Ipopt** via NLPModelsIpopt. The `@time` macro
-reports wall-clock time, which on first call includes Julia's JIT compilation.
-
-```@example main
-initial_guess = (state=[0.0, 0.0, 0.0, 0.1], control=0.0)
-
-@time sol = solve(cartpole_ocp; display=true, grid_size=100, init=initial_guess)
-```
-
-### Results
-
-```@example main
-println("--- Optimal Limit Cycle Found ---")
-println("Total Control Energy (∫F² dt): ", sol.objective)
-
-tsol = time_grid(sol)
-Xsol = state(sol).(tsol)
-Fsol = control(sol).(tsol)
-
-# Stack state vectors into a 4×N matrix and unpack rows
-xsol, θsol, vsol, ωsol = eachrow(reduce(hcat, Xsol))
-
-qplot  = plot(tsol, [xsol θsol], label=["x" "θ"],   title="Configuration")
-dqplot = plot(tsol, [vsol ωsol], label=["v" "ω"],   title="Velocities")
-uplot  = plot(tsol, Fsol,        label="u",          title="Control", linetype=:steppost)
-
-plot(qplot, dqplot, uplot, layout=3, size=(800, 600))
-```
-
-The three panels show the cart position ``x`` and pendulum angle ``\theta``, the
-corresponding velocities ``\dot{x}`` and ``\dot\theta``, and the optimal control force
-``u``. The periodicity of the state trajectory confirms that a genuine limit cycle has been
-found.
-
-### Animation
-
-The animation below shows the cart-pole evolving along the optimal limit-cycle trajectory.
-The blue cart slides on the horizontal rail while the pendulum swings around the upright
-equilibrium. A ghost trace follows the red bob to make the motion easier to read. Use the
-controls to play, pause, or reset the animation and to adjust the playback speed.
-
-```@example main
-# Serialise trajectory for the JS animation (manual, no extra dependencies)
-json_t  = "[" * join(string.(tsol),  ",") * "]"
-json_x  = "[" * join(string.(xsol),  ",") * "]"
-json_th = "[" * join(string.(θsol),  ",") * "]"
-
-uid    = string(rand(UInt32), base=16)   # unique suffix → safe for multi-example pages
-tf_str = string(round(tsol[end], digits=3))
-
-html_str = """
-<div style="text-align:center;font-family:sans-serif;margin:1.5em 0;">
-  <canvas id="cpC$(uid)" width="700" height="380"
-    style="border:1px solid #ccc;border-radius:6px;background:#fafafa;display:block;margin:0 auto;"></canvas>
-  <div style="margin:6px 0;font-size:0.9em;color:#444;">
-    <span id="cpT$(uid)">t = 0.000 / $(tf_str) s</span>
-  </div>
-  <div style="display:flex;justify-content:center;gap:10px;align-items:center;flex-wrap:wrap;margin-top:4px;">
-    <button id="cpPl$(uid)" style="padding:5px 14px;cursor:pointer;">&#9654; Play</button>
-    <button id="cpPa$(uid)" style="padding:5px 14px;cursor:pointer;">&#9646;&#9646; Pause</button>
-    <button id="cpRe$(uid)" style="padding:5px 14px;cursor:pointer;">&#x21BA; Reset</button>
-    <label style="font-size:0.9em;">Speed:
-      <input type="range" id="cpSp$(uid)" min="0.25" max="3" step="0.25" value="1"
-             style="vertical-align:middle;width:100px;">
-      <span id="cpSv$(uid)">1.00&times;</span>
-    </label>
-  </div>
-</div>
-<script>
-(function(){
-  var T=$(json_t), X=$(json_x), TH=$(json_th), L=$(l_val);
-  var cv=document.getElementById("cpC$(uid)");
-  var ctx=cv.getContext("2d");
-  var W=cv.width, H=cv.height;
-  var railY=H*0.62, sc=80, cW=60, cH=28;
-  var tSim=0, spd=1, rid=null, lts=null, ghost=[];
-  function lerp(a,b,t){return a+(b-a)*t;}
-  function interp(tq){
-    var n=T.length;
-    if(tq<=T[0])return{x:X[0],th:TH[0]};
-    if(tq>=T[n-1])return{x:X[n-1],th:TH[n-1]};
-    var lo=0,hi=n-2;
-    while(lo<hi){var m=(lo+hi)>>1;if(T[m+1]<tq)lo=m+1;else hi=m;}
-    var a=(tq-T[lo])/(T[lo+1]-T[lo]);
-    return{x:lerp(X[lo],X[lo+1],a),th:lerp(TH[lo],TH[lo+1],a)};
-  }
-  function draw(tq){
-    var d=interp(tq);
-    var cx=W/2+d.x*sc, cy=railY;
-    var px=cx, py=cy-cH/2;
-    var bx=px+L*sc*Math.sin(d.th), by=py-L*sc*Math.cos(d.th);
-    ghost.push({x:bx,y:by}); if(ghost.length>25)ghost.shift();
-    ctx.clearRect(0,0,W,H);
-    ctx.beginPath();ctx.moveTo(20,railY);ctx.lineTo(W-20,railY);
-    ctx.strokeStyle="#aaa";ctx.lineWidth=3;ctx.stroke();
-    var wr=8;
-    [[cx-18,cy+cH/2+wr],[cx+18,cy+cH/2+wr]].forEach(function(p){
-      ctx.beginPath();ctx.arc(p[0],p[1],wr,0,2*Math.PI);
-      ctx.fillStyle="#888";ctx.fill();
-    });
-    ctx.beginPath();ctx.rect(cx-cW/2,cy-cH/2,cW,cH);
-    ctx.fillStyle="#4a90d9";ctx.fill();
-    ctx.strokeStyle="#2c5f8a";ctx.lineWidth=1.5;ctx.stroke();
-    ghost.forEach(function(g,i){
-      ctx.beginPath();ctx.arc(g.x,g.y,5,0,2*Math.PI);
-      ctx.fillStyle="rgba(220,60,60,"+((i+1)/ghost.length*0.35)+")";ctx.fill();
-    });
-    ctx.beginPath();ctx.moveTo(px,py);ctx.lineTo(bx,by);
-    ctx.strokeStyle="#222";ctx.lineWidth=3;ctx.stroke();
-    ctx.beginPath();ctx.arc(px,py,5,0,2*Math.PI);ctx.fillStyle="#555";ctx.fill();
-    ctx.beginPath();ctx.arc(bx,by,10,0,2*Math.PI);
-    ctx.fillStyle="#dc3c3c";ctx.fill();ctx.strokeStyle="#8a1a1a";ctx.lineWidth=1.5;ctx.stroke();
-    document.getElementById("cpT$(uid)").textContent=
-      "t = "+tq.toFixed(3)+" / "+T[T.length-1].toFixed(2)+" s";
-  }
-  function step(ts){
-    if(lts!==null){
-      var dt=(ts-lts)*0.001*spd; tSim+=dt;
-      var tf=T[T.length-1]; if(tSim>tf){tSim-=tf;ghost=[];}
-    }
-    lts=ts; draw(tSim); rid=requestAnimationFrame(step);
-  }
-  document.getElementById("cpPl$(uid)").onclick=function(){
-    if(rid===null){lts=null;rid=requestAnimationFrame(step);}
-  };
-  document.getElementById("cpPa$(uid)").onclick=function(){
-    if(rid!==null){cancelAnimationFrame(rid);rid=null;}
-  };
-  document.getElementById("cpRe$(uid)").onclick=function(){
-    if(rid!==null){cancelAnimationFrame(rid);rid=null;}
-    tSim=0;ghost=[];lts=null;draw(0);
-  };
-  document.getElementById("cpSp$(uid)").oninput=function(){
-    spd=parseFloat(this.value);
-    document.getElementById("cpSv$(uid)").textContent=spd.toFixed(2)+"\u00d7";
-  };
-  draw(0);
-})();
-</script>
-"""
-
-HTML(html_str)
-```
\ No newline at end of file
diff --git a/docs/src/tutorial-symbolics-save.md b/docs/src/tutorial-symbolics-save.md
new file mode 100644
index 0000000..6c7d2dd
--- /dev/null
+++ b/docs/src/tutorial-symbolics-save.md
@@ -0,0 +1,313 @@
+# Cart-Pole Limit Cycle via Symbolic Lagrangian Mechanics
+
+```@meta
+Draft = true
+```
+
+This tutorial demonstrates how to combine **symbolic derivation of equations of motion** (via
+[Symbolics.jl](https://symbolics.juliasymbolics.org/)) with **direct optimal control**
+(via [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl/)) to find a
+**periodic orbit** (limit cycle) for a cart-pole system.
+
+The key idea is to let the computer do the hard mechanics: we write the Lagrangian in a few
+lines, and the Euler–Lagrange equations — including mass-matrix inversion — are derived
+automatically.
+
+## The Cart-Pole System
+
+The system consists of a cart of mass ``m_c`` sliding on a frictionless horizontal rail, with
+a rigid pendulum of mass ``m_p`` and length ``l`` attached to it. A horizontal force ``u``
+(the control input) acts on the cart.
+
+```@raw html
+<figure>
+  <img src="cartpole_sketch.svg" alt="Cart-pole diagram" width="350"/>
+  <figcaption>Fig. 1 — Cart-pole system. The angle θ is measured from the upright position.</figcaption>
+</figure>
+```
+
+The configuration vector is ``q = (x,\, \theta)^\top``, where ``x`` is the cart position and
+``\theta = 0`` corresponds to the **upright** (unstable) equilibrium of the pendulum.
+
+### Positions
+
+The Cartesian positions of the two bodies are:
+
+```math
+p_c = \begin{pmatrix} x \\ 0 \end{pmatrix}, \qquad
+p_p = \begin{pmatrix} x + l\sin\theta \\ l\cos\theta \end{pmatrix}.
+```
+
+### Lagrangian
+
+The kinetic and potential energies are:
+
+```math
+T = \tfrac{1}{2}m_c\,\|\dot{p}_c\|^2 + \tfrac{1}{2}m_p\,\|\dot{p}_p\|^2
+  = \tfrac{1}{2}(m_c+m_p)\dot{x}^2
+    + m_p l\,\dot{x}\dot{\theta}\cos\theta
+    + \tfrac{1}{2}m_p l^2\dot{\theta}^2,
+```
+
+```math
+V = m_p\,g\,l\cos\theta.
+```
+
+The Lagrangian is ``\mathcal{L} = T - V``, and the virtual work of the control force gives the
+generalised force vector ``Q = (u,\, 0)^\top``.
+
+### Euler–Lagrange Equations
+
+The equations of motion follow from:
+
+```math
+\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{q}_i}
+- \frac{\partial \mathcal{L}}{\partial q_i} = Q_i, \qquad i = 1,2.
+```
+
+They can be written in the standard **manipulator form**:
+
+```math
+M(q)\,\ddot{q} = -C(q,\dot{q}) + \tau(q, u),
+```
+
+where the symmetric positive-definite **mass matrix** is:
+
+```math
+M(q) =
+\begin{pmatrix}
+  m_c + m_p & m_p l \cos\theta \\
+  m_p l \cos\theta & m_p l^2
+\end{pmatrix},
+```
+
+and the right-hand side collects Coriolis/gravity terms and the control torque. Instead of
+deriving these by hand, we rely on Symbolics.jl to compute and **analytically invert**
+``M(q)`` for us.
+
+### State-Space Form
+
+Defining the state ``X = (x,\,\theta,\,\dot{x},\,\dot{\theta})^\top``, the equations of
+motion become the first-order system:
+
+```math
+\dot{X}(t) = f\!\left(X(t),\, u(t)\right) =
+\begin{pmatrix}
+  \dot{x} \\ \dot{\theta} \\ M^{-1}(q)\bigl(-C(q,\dot{q}) + \tau(q,u)\bigr)
+\end{pmatrix}.
+```
+
+## The Optimal Control Problem
+
+We look for a **limit cycle** of the nonlinear dynamics: a trajectory that returns exactly to
+its initial condition after a fixed period ``t_f``. The cost penalises the total control
+energy:
+
+```math
+\min_{u(\cdot)}\; \int_0^{t_f} u(t)^2\,\mathrm{d}t
+```
+
+```math
+\text{subject to} \quad \dot{X}(t) = f(X(t), u(t)), \quad t \in [0, t_f],
+```
+
+```math
+X(t_f) = X(0).
+```
+
+The periodicity constraint ``X(t_f) = X(0)`` — combined with a non-trivial initial condition
+— forces the solver to find an orbit rather than the trivial rest solution.
+
+## Implementation
+
+### Setup & Imports
+
+```@example main
+# ==========================================
+# Setup & Imports
+# ==========================================
+import Pkg
+Pkg.activate("Control.jl/ModelingToolkit.jl")
+
+using OptimalControl
+using Plots
+using StaticArrays
+
+import NLPModelsIpopt
+
+using Symbolics
+```
+
+### Physical Parameters and Symbolic Variables
+
+We declare all parameters both as numerical constants (for the final function evaluation) and
+as symbolic variables (for the Lagrangian computation).
+
+```@example main
+# ==========================================
+# 1. Variables and Parameters
+# ==========================================
+# Physical Constants
+const m_c_val = 5.0
+const m_p_val = 1.0
+const l_val = 2.0
+const g_val = 9.81
+const tf_val = 2.0
+
+# Symbolic Variables
+@variables t
+D = Differential(t)
+@variables m_c m_p l g u
+@variables x(t) θ(t)
+@variables v ω dv dω # Static variables to solve for accelerations
+
+q = [x, θ]
+nothing # hide
+```
+
+### Automated Kinematics and Lagrangian
+
+We express the positions, kinetic energy, potential energy, and virtual work symbolically.
+The time derivatives ``\dot{p}_c``, ``\dot{p}_p`` are computed automatically by `D.(...)`.
+
+```@example main
+# ==========================================
+# 2. Automated Kinematics & Jacobians
+# ==========================================
+p_c = [x, 0.0]
+p_p = [x + l * sin(θ), l * cos(θ)]
+F = [u, 0.0]
+
+T = 0.5 * m_c * sum(D.(p_c) .^ 2) + 0.5 * m_p * sum(D.(p_p) .^ 2)
+V = g * (m_p * p_p[2])
+P_non_conservative = transpose(D.(p_c)) * F
+nothing # hide
+```
+
+### Euler–Lagrange Equations and Mass-Matrix Inversion
+
+Starting from ``\mathcal{L} = T - V``, Symbolics.jl computes the three terms of the
+Euler–Lagrange equations and assembles the residual vector. Substituting static aliases
+``(v,\omega,\dot{v},\dot\omega)`` for the time derivatives makes it possible to identify
+the mass matrix ``M`` as the Jacobian of the residual with respect to the accelerations
+``(\dot{v}, \dot\omega)``. The system ``M\,a = -b`` is then solved analytically.
+
+```@example main
+# ==========================================
+# 3. Euler-Lagrange & Matrix Inversion
+# ==========================================
+L = T - V
+
+A = D.(Symbolics.gradient(L, D.(q)))
+B = Symbolics.gradient(L, q)
+Q = Symbolics.gradient(P_non_conservative, D.(q))
+
+# Euler-Lagrange: d/dt( dL/dq_dot ) - dL/dq = Forces
+el_eqs = expand_derivatives.(A - B - Q)
+
+# Freeze time derivatives into static algebraic variables
+sub_rules = Dict(D(x) => v, D(θ) => ω, D(D(x)) => dv, D(D(θ)) => dω)
+res = Symbolics.substitute.(el_eqs, (sub_rules,))
+
+# Extract Mass Matrix and Bias vector
+Mass = Symbolics.jacobian(res, [dv, dω])
+bias = Symbolics.substitute.(res, (Dict(dv => 0.0, dω => 0.0),))
+
+# Analytically invert the Mass Matrix (Symbolics handles 2x2 natively)
+accel = Mass \ (-bias)
+accel = Symbolics.simplify_fractions.(accel)
+
+# The fully explicit state derivatives: X_dot = [v, ω, v_dot, ω_dot]
+dx_dt = [v, ω, accel[1], accel[2]]
+nothing # hide
+```
+
+### Code Generation
+
+`build_function` compiles the symbolic expression `dx_dt` into a native Julia function.
+The `force_SA=true` flag generates a **StaticArrays** kernel, which avoids heap allocations
+inside the ODE right-hand side — important for solver performance.
+
+```@example main
+# ==========================================
+# 4. Extract Explicit Dynamics into a Julia Function
+# ==========================================
+f_expr = build_function(dx_dt, [x, θ, v, ω], u, [m_c, m_p, l, g];
+    expression=Val{false}, force_SA=true)
+f_mtk = f_expr[1]
+
+const p_vals = [m_c_val, m_p_val, l_val, g_val]
+
+cartpole_dynamics(X, U) = f_mtk(X, U, p_vals)
+nothing # hide
+```
+
+### Optimal Control Problem Definition
+
+We now formulate the optimal control problem using the `@def` macro from OptimalControl.jl.
+The periodicity condition is encoded as ``X(t_f) - X(0) = 0``, and the dynamics plug
+directly into the `Ẋ(t) == ...` constraint.
+
+```@example main
+# ==========================================
+# 5. OptimalControl.jl dynamic optimization problem
+# ==========================================
+@def cartpole_ocp begin
+    t ∈ [0, tf_val], time
+    X = (x_ocp, θ_ocp, v_ocp, ω_ocp) ∈ R⁴, state
+    F_ctrl ∈ R, control
+
+    X(0) == [0, 0, 0, 0.1]
+    X(tf_val) - X(0) == [0, 0, 0, 0]
+
+    # Inject our explicitly inverted Symbolic dynamics!
+    Ẋ(t) == cartpole_dynamics(X(t), F_ctrl(t))
+
+    ∫(F_ctrl(t)^2) → min
+end
+```
+
+### Solving the NLP
+
+The problem is transcribed into a nonlinear program using direct collocation on a uniform
+grid of 100 intervals, then handed to **Ipopt** via NLPModelsIpopt.
+
+```@example main
+# ==========================================
+# 6. Solvers and Collocation
+# ==========================================
+initial_guess = (state=[0.0, 0.0, 0.0, 0.1], control=0.0)
+
+sol = solve(cartpole_ocp; display=true, grid_size=100, init=initial_guess)
+```
+
+### Results
+
+```@example main
+# ==========================================
+# 7. Extracting Results
+# ==========================================
+println("--- Optimal Limit Cycle Found ---")
+println("Total Control Energy (∫F² dt): ", sol.objective)
+
+tsol = time_grid(sol)
+# tsol = range(sol.model.times.initial.time, sol.model.times.final.time, 1001)
+Xsol = state(sol).(tsol)
+Fsol = control(sol).(tsol)
+
+xsol = [X[1] for X in Xsol]
+θsol = [X[2] for X in Xsol]
+vsol = [X[3] for X in Xsol]
+ωsol = [X[4] for X in Xsol]
+
+qplot = plot(tsol, [xsol θsol], label=["x" "θ"], title="Configuration")
+dqplot = plot(tsol, [vsol ωsol], label=["v" "ω"], title="Velocities")
+uplot = plot(tsol, Fsol, label="u", title="Control", linetype=:steppost)
+
+plot(qplot, dqplot, uplot, layout=3, size=(800, 600))
+```
+
+The three panels show, from top to bottom: the cart position ``x`` and pendulum angle
+``\theta``, the corresponding velocities ``\dot{x}`` and ``\dot\theta``, and the optimal
+control force ``u``. The periodicity of the state trajectory confirms that a genuine limit
+cycle has been found.
diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index 515278d..ec7c303 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -53,8 +53,14 @@ T = \tfrac{1}{2}m_c\,\|\dot{p}_c\|^2 + \tfrac{1}{2}m_p\,\|\dot{p}_p\|^2
 V = m_p\,g\,l\cos\theta.
 ```
 
-The Lagrangian is ``\mathcal{L} = T - V``, and the virtual work of the control force gives the
-generalised force vector ``Q = (u,\, 0)^\top``.
+The Lagrangian is ``\mathcal{L} = T - V``. The virtual work of the control force ``u``
+acting on the cart gives the generalised force vector:
+
+```math
+W_\text{ctrl} = F \cdot \dot{p}_c = u\,\dot{x},
+\quad\Longrightarrow\quad
+Q = \frac{\partial W_\text{ctrl}}{\partial \dot{q}} = \begin{pmatrix} u \\ 0 \end{pmatrix}.
+```
 
 ### Euler–Lagrange Equations
 
@@ -112,74 +118,70 @@ energy:
 ```
 
 ```math
-X(t_f) = X(0).
+X(0) = X(t_f) = (0,\, 0,\, 0,\, 0.1)^\top.
 ```
 
-The periodicity constraint ``X(t_f) = X(0)`` — combined with a non-trivial initial condition
-— forces the solver to find an orbit rather than the trivial rest solution.
+The common boundary condition ``X(0) = X(t_f)`` — combined with a non-trivial initial
+angular velocity — forces the solver to find an orbit rather than the trivial rest solution.
 
 ## Implementation
 
 ### Setup & Imports
 
-```@example main
-# ==========================================
-# Setup & Imports
-# ==========================================
+```@setup main
+# Project environment for this tutorial — edit the path for your local setup.
 import Pkg
-Pkg.activate("Control.jl/ModelingToolkit.jl")
+Pkg.activate(".")
+```
 
+```@example main
 using OptimalControl
 using Plots
 using StaticArrays
-
-import NLPModelsIpopt
-
+using NLPModelsIpopt
 using Symbolics
 ```
 
 ### Physical Parameters and Symbolic Variables
 
 We declare all parameters both as numerical constants (for the final function evaluation) and
-as symbolic variables (for the Lagrangian computation).
+as symbolic variables (for the Lagrangian computation). Note that the symbolic time variable
+`t` is used only inside the Symbolics.jl derivation and does not interfere with the time
+variable `t` introduced later by the `@def` macro.
 
 ```@example main
-# ==========================================
-# 1. Variables and Parameters
-# ==========================================
-# Physical Constants
+# Physical constants
 const m_c_val = 5.0
 const m_p_val = 1.0
-const l_val = 2.0
-const g_val = 9.81
-const tf_val = 2.0
+const l_val   = 2.0
+const g_val   = 9.81
+const tf_val  = 2.0
 
-# Symbolic Variables
+# Symbolic variables
 @variables t
 D = Differential(t)
 @variables m_c m_p l g u
 @variables x(t) θ(t)
-@variables v ω dv dω # Static variables to solve for accelerations
+@variables v ω dv dω   # static aliases for velocities and accelerations
 
 q = [x, θ]
+nothing # hide
 ```
 
 ### Automated Kinematics and Lagrangian
 
-We express the positions, kinetic energy, potential energy, and virtual work symbolically.
+We express the positions, kinetic energy, potential energy, and virtual power symbolically.
 The time derivatives ``\dot{p}_c``, ``\dot{p}_p`` are computed automatically by `D.(...)`.
 
 ```@example main
-# ==========================================
-# 2. Automated Kinematics & Jacobians
-# ==========================================
 p_c = [x, 0.0]
 p_p = [x + l * sin(θ), l * cos(θ)]
-F = [u, 0.0]
+F   = [u, 0.0]
 
-T = 0.5 * m_c * sum(D.(p_c) .^ 2) + 0.5 * m_p * sum(D.(p_p) .^ 2)
-V = g * (m_p * p_p[2])
-P_non_conservative = transpose(D.(p_c)) * F
+T      = 0.5 * m_c * sum(D.(p_c) .^ 2) + 0.5 * m_p * sum(D.(p_p) .^ 2)
+V      = g * (m_p * p_p[2])
+W_ctrl = transpose(D.(p_c)) * F   # virtual power of the control force
+nothing # hide
 ```
 
 ### Euler–Lagrange Equations and Mass-Matrix Inversion
@@ -191,72 +193,63 @@ the mass matrix ``M`` as the Jacobian of the residual with respect to the accele
 ``(\dot{v}, \dot\omega)``. The system ``M\,a = -b`` is then solved analytically.
 
 ```@example main
-# ==========================================
-# 3. Euler-Lagrange & Matrix Inversion
-# ==========================================
 L = T - V
 
 A = D.(Symbolics.gradient(L, D.(q)))
 B = Symbolics.gradient(L, q)
-Q = Symbolics.gradient(P_non_conservative, D.(q))
+Q = Symbolics.gradient(W_ctrl, D.(q))   # generalised forces
 
-# Euler-Lagrange: d/dt( dL/dq_dot ) - dL/dq = Forces
+# Euler-Lagrange residual: d/dt(∂L/∂q̇) - ∂L/∂q - Q = 0
 el_eqs = expand_derivatives.(A - B - Q)
 
 # Freeze time derivatives into static algebraic variables
 sub_rules = Dict(D(x) => v, D(θ) => ω, D(D(x)) => dv, D(D(θ)) => dω)
 res = Symbolics.substitute.(el_eqs, (sub_rules,))
 
-# Extract Mass Matrix and Bias vector
+# Identify mass matrix M and bias vector b such that M·[dv; dω] + b = 0
 Mass = Symbolics.jacobian(res, [dv, dω])
 bias = Symbolics.substitute.(res, (Dict(dv => 0.0, dω => 0.0),))
 
-# Analytically invert the Mass Matrix (Symbolics handles 2x2 natively)
-accel = Mass \ (-bias)
-accel = Symbolics.simplify_fractions.(accel)
+# Analytically invert the 2×2 mass matrix
+accel = Symbolics.simplify_fractions.(Mass \ (-bias))
 
-# The fully explicit state derivatives: X_dot = [v, ω, v_dot, ω_dot]
+# Fully explicit state derivative: Ẋ = [v, ω, v̇, ω̇]
 dx_dt = [v, ω, accel[1], accel[2]]
+nothing # hide
 ```
 
 ### Code Generation
 
 `build_function` compiles the symbolic expression `dx_dt` into a native Julia function.
 The `force_SA=true` flag generates a **StaticArrays** kernel, which avoids heap allocations
-inside the ODE right-hand side — important for solver performance.
+inside the ODE right-hand side — important for solver performance. Parameters are stored as
+an `SVector` to match the generated kernel's expected input type.
 
 ```@example main
-# ==========================================
-# 4. Extract Explicit Dynamics into a Julia Function
-# ==========================================
 f_expr = build_function(dx_dt, [x, θ, v, ω], u, [m_c, m_p, l, g];
     expression=Val{false}, force_SA=true)
-f_mtk = f_expr[1]
+f_cartpole = f_expr[1]   # out-of-place variant: (state, u, params) → SVector
 
-const p_vals = [m_c_val, m_p_val, l_val, g_val]
+const p_vals = SA[m_c_val, m_p_val, l_val, g_val]   # SVector, matches SA kernel
 
-cartpole_dynamics(X, U) = f_mtk(X, U, p_vals)
+cartpole_dynamics(X, U) = f_cartpole(X, U, p_vals)
+nothing # hide
 ```
 
 ### Optimal Control Problem Definition
 
 We now formulate the optimal control problem using the `@def` macro from OptimalControl.jl.
-The periodicity condition is encoded as ``X(t_f) - X(0) = 0``, and the dynamics plug
-directly into the `Ẋ(t) == ...` constraint.
+The boundary conditions ``X(0) = X(t_f)`` encode the periodicity of the orbit directly.
 
 ```@example main
-# ==========================================
-# 5. OptimalControl.jl dynamic optimization problem
-# ==========================================
 @def cartpole_ocp begin
     t ∈ [0, tf_val], time
     X = (x_ocp, θ_ocp, v_ocp, ω_ocp) ∈ R⁴, state
     F_ctrl ∈ R, control
 
-    X(0) == [0, 0, 0, 0.1]
-    X(tf_val) - X(0) == [0, 0, 0, 0]
+    X(0)      == [0, 0, 0, 0.1]
+    X(tf_val) == [0, 0, 0, 0.1]   # periodicity: X(tf) = X(0)
 
-    # Inject our explicitly inverted Symbolic dynamics!
     Ẋ(t) == cartpole_dynamics(X(t), F_ctrl(t))
 
     ∫(F_ctrl(t)^2) → min
@@ -266,44 +259,146 @@ end
 ### Solving the NLP
 
 The problem is transcribed into a nonlinear program using direct collocation on a uniform
-grid of 100 intervals, then handed to **Ipopt** via NLPModelsIpopt.
+grid of 100 intervals, then handed to **Ipopt** via NLPModelsIpopt. The `@time` macro
+reports wall-clock time, which on first call includes Julia's JIT compilation.
 
 ```@example main
-# ==========================================
-# 6. Solvers and Collocation
-# ==========================================
 initial_guess = (state=[0.0, 0.0, 0.0, 0.1], control=0.0)
 
-@time sol = solve(cartpole_ocp; display=true, grid_size=100, init=initial_guess)
+sol = solve(cartpole_ocp; display=true, grid_size=100, init=initial_guess)
 ```
 
 ### Results
 
 ```@example main
-# ==========================================
-# 7. Extracting Results
-# ==========================================
 println("--- Optimal Limit Cycle Found ---")
 println("Total Control Energy (∫F² dt): ", sol.objective)
 
 tsol = time_grid(sol)
-# tsol = range(sol.model.times.initial.time, sol.model.times.final.time, 1001)
 Xsol = state(sol).(tsol)
 Fsol = control(sol).(tsol)
 
-xsol = [X[1] for X in Xsol]
-θsol = [X[2] for X in Xsol]
-vsol = [X[3] for X in Xsol]
-ωsol = [X[4] for X in Xsol]
+# Stack state vectors into a 4×N matrix and unpack rows
+xsol, θsol, vsol, ωsol = eachrow(reduce(hcat, Xsol))
 
-qplot = plot(tsol, [xsol θsol], label=["x" "θ"], title="Configuration")
-dqplot = plot(tsol, [vsol ωsol], label=["v" "ω"], title="Velocities")
-uplot = plot(tsol, Fsol, label="u", title="Control", linetype=:steppost)
+qplot  = plot(tsol, [xsol θsol], label=["x" "θ"],   title="Configuration")
+dqplot = plot(tsol, [vsol ωsol], label=["v" "ω"],   title="Velocities")
+uplot  = plot(tsol, Fsol,        label="u",          title="Control", linetype=:steppost)
 
 plot(qplot, dqplot, uplot, layout=3, size=(800, 600))
 ```
 
-The three panels show, from top to bottom: the cart position ``x`` and pendulum angle
-``\theta``, the corresponding velocities ``\dot{x}`` and ``\dot\theta``, and the optimal
-control force ``u``. The periodicity of the state trajectory confirms that a genuine limit
-cycle has been found.
+The three panels show the cart position ``x`` and pendulum angle ``\theta``, the
+corresponding velocities ``\dot{x}`` and ``\dot\theta``, and the optimal control force
+``u``. The periodicity of the state trajectory confirms that a genuine limit cycle has been
+found.
+
+### Animation
+
+The animation below shows the cart-pole evolving along the optimal limit-cycle trajectory.
+The blue cart slides on the horizontal rail while the pendulum swings around the upright
+equilibrium. A ghost trace follows the red bob to make the motion easier to read. Use the
+controls to play, pause, or reset the animation and to adjust the playback speed.
+
+```@setup main
+# Serialise trajectory for the JS animation (manual, no extra dependencies)
+json_t  = "[" * join(string.(tsol),  ",") * "]"
+json_x  = "[" * join(string.(xsol),  ",") * "]"
+json_th = "[" * join(string.(θsol),  ",") * "]"
+
+uid    = string(rand(UInt32), base=16)   # unique suffix → safe for multi-example pages
+tf_str = string(round(tsol[end], digits=3))
+
+html_str = """
+<div style="text-align:center;font-family:sans-serif;margin:1.5em 0;">
+  <canvas id="cpC$(uid)" width="700" height="380"
+    style="border:1px solid #ccc;border-radius:6px;background:#fafafa;display:block;margin:0 auto;"></canvas>
+  <div style="margin:6px 0;font-size:0.9em;color:#444;">
+    <span id="cpT$(uid)">t = 0.000 / $(tf_str) s</span>
+  </div>
+  <div style="display:flex;justify-content:center;gap:10px;align-items:center;flex-wrap:wrap;margin-top:4px;">
+    <button id="cpPl$(uid)" style="padding:5px 14px;cursor:pointer;">&#9654; Play</button>
+    <button id="cpPa$(uid)" style="padding:5px 14px;cursor:pointer;">&#9646;&#9646; Pause</button>
+    <button id="cpRe$(uid)" style="padding:5px 14px;cursor:pointer;">&#x21BA; Reset</button>
+    <label style="font-size:0.9em;">Speed:
+      <input type="range" id="cpSp$(uid)" min="0.25" max="3" step="0.25" value="1"
+             style="vertical-align:middle;width:100px;">
+      <span id="cpSv$(uid)">1.00&times;</span>
+    </label>
+  </div>
+</div>
+<script>
+(function(){
+  var T=$(json_t), X=$(json_x), TH=$(json_th), L=$(l_val);
+  var cv=document.getElementById("cpC$(uid)");
+  var ctx=cv.getContext("2d");
+  var W=cv.width, H=cv.height;
+  var railY=H*0.62, sc=80, cW=60, cH=28;
+  var tSim=0, spd=1, rid=null, lts=null, ghost=[];
+  function lerp(a,b,t){return a+(b-a)*t;}
+  function interp(tq){
+    var n=T.length;
+    if(tq<=T[0])return{x:X[0],th:TH[0]};
+    if(tq>=T[n-1])return{x:X[n-1],th:TH[n-1]};
+    var lo=0,hi=n-2;
+    while(lo<hi){var m=(lo+hi)>>1;if(T[m+1]<tq)lo=m+1;else hi=m;}
+    var a=(tq-T[lo])/(T[lo+1]-T[lo]);
+    return{x:lerp(X[lo],X[lo+1],a),th:lerp(TH[lo],TH[lo+1],a)};
+  }
+  function draw(tq){
+    var d=interp(tq);
+    var cx=W/2+d.x*sc, cy=railY;
+    var px=cx, py=cy-cH/2;
+    var bx=px+L*sc*Math.sin(d.th), by=py-L*sc*Math.cos(d.th);
+    ghost.push({x:bx,y:by}); if(ghost.length>25)ghost.shift();
+    ctx.clearRect(0,0,W,H);
+    ctx.beginPath();ctx.moveTo(20,railY);ctx.lineTo(W-20,railY);
+    ctx.strokeStyle="#aaa";ctx.lineWidth=3;ctx.stroke();
+    var wr=8;
+    [[cx-18,cy+cH/2+wr],[cx+18,cy+cH/2+wr]].forEach(function(p){
+      ctx.beginPath();ctx.arc(p[0],p[1],wr,0,2*Math.PI);
+      ctx.fillStyle="#888";ctx.fill();
+    });
+    ctx.beginPath();ctx.rect(cx-cW/2,cy-cH/2,cW,cH);
+    ctx.fillStyle="#4a90d9";ctx.fill();
+    ctx.strokeStyle="#2c5f8a";ctx.lineWidth=1.5;ctx.stroke();
+    ghost.forEach(function(g,i){
+      ctx.beginPath();ctx.arc(g.x,g.y,5,0,2*Math.PI);
+      ctx.fillStyle="rgba(220,60,60,"+((i+1)/ghost.length*0.35)+")";ctx.fill();
+    });
+    ctx.beginPath();ctx.moveTo(px,py);ctx.lineTo(bx,by);
+    ctx.strokeStyle="#222";ctx.lineWidth=3;ctx.stroke();
+    ctx.beginPath();ctx.arc(px,py,5,0,2*Math.PI);ctx.fillStyle="#555";ctx.fill();
+    ctx.beginPath();ctx.arc(bx,by,10,0,2*Math.PI);
+    ctx.fillStyle="#dc3c3c";ctx.fill();ctx.strokeStyle="#8a1a1a";ctx.lineWidth=1.5;ctx.stroke();
+    document.getElementById("cpT$(uid)").textContent=
+      "t = "+tq.toFixed(3)+" / "+T[T.length-1].toFixed(2)+" s";
+  }
+  function step(ts){
+    if(lts!==null){
+      var dt=(ts-lts)*0.001*spd; tSim+=dt;
+      var tf=T[T.length-1]; if(tSim>tf){tSim-=tf;ghost=[];}
+    }
+    lts=ts; draw(tSim); rid=requestAnimationFrame(step);
+  }
+  document.getElementById("cpPl$(uid)").onclick=function(){
+    if(rid===null){lts=null;rid=requestAnimationFrame(step);}
+  };
+  document.getElementById("cpPa$(uid)").onclick=function(){
+    if(rid!==null){cancelAnimationFrame(rid);rid=null;}
+  };
+  document.getElementById("cpRe$(uid)").onclick=function(){
+    if(rid!==null){cancelAnimationFrame(rid);rid=null;}
+    tSim=0;ghost=[];lts=null;draw(0);
+  };
+  document.getElementById("cpSp$(uid)").oninput=function(){
+    spd=parseFloat(this.value);
+    document.getElementById("cpSv$(uid)").textContent=spd.toFixed(2)+"\u00d7";
+  };
+  draw(0);
+})();
+</script>
+"""
+
+HTML(html_str)
+```

From 61cc7ec6f693c8363e75263bf894336969da581e Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 14:33:00 +0200
Subject: [PATCH 09/21] typos

---
 .github/workflows/SpellCheck.yml |  2 ++
 _typos.toml                      | 12 ++++++++++++
 2 files changed, 14 insertions(+)
 create mode 100644 _typos.toml

diff --git a/.github/workflows/SpellCheck.yml b/.github/workflows/SpellCheck.yml
index fe1c7c4..5513363 100644
--- a/.github/workflows/SpellCheck.yml
+++ b/.github/workflows/SpellCheck.yml
@@ -7,3 +7,5 @@ on:
 jobs:
   call:
     uses: control-toolbox/CTActions/.github/workflows/spell-check.yml@main
+    with:
+      config-path: '_typos.toml'
diff --git a/_typos.toml b/_typos.toml
new file mode 100644
index 0000000..6263e4b
--- /dev/null
+++ b/_typos.toml
@@ -0,0 +1,12 @@
+[default]
+locale = "en"
+extend-ignore-re = [
+    "ded",
+]
+
+[files]
+extend-exclude = [
+    "*.json",
+    "*.toml",
+    "*.svg",
+]
\ No newline at end of file

From 65732525f9748a9e9c8ee22ca704eee3a3a8aa74 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 14:34:12 +0200
Subject: [PATCH 10/21] add svg

---
 docs/src/tutorial-symbolics.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index ec7c303..ed4947f 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -21,7 +21,7 @@ a rigid pendulum of mass ``m_p`` and length ``l`` attached to it. A horizontal f
 
 ```@raw html
 <figure>
-  <img src="cartpole_sketch.svg" alt="Cart-pole diagram" width="350"/>
+  <img src="./assets/cartpole_sketch.svg" alt="Cart-pole diagram" width="350"/>
   <figcaption>Fig. 1 — Cart-pole system. The angle θ is measured from the upright position.</figcaption>
 </figure>
 ```

From 68226c286a883462b50bc3ed42fb9fc3a3293544 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 14:35:38 +0200
Subject: [PATCH 11/21] review anim

---
 docs/src/tutorial-symbolics.md | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index ed4947f..c9a52d6 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -399,6 +399,8 @@ html_str = """
 })();
 </script>
 """
+```
 
-HTML(html_str)
+```@example main
+HTML(html_str) # hide
 ```

From 723f6735ce8610fea7bf72551adb1793c380e5c6 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 14:36:38 +0200
Subject: [PATCH 12/21] Add typos config and update spell-check workflow

---
 _typos.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/_typos.toml b/_typos.toml
index 6263e4b..d004eee 100644
--- a/_typos.toml
+++ b/_typos.toml
@@ -1,7 +1,7 @@
 [default]
 locale = "en"
 extend-ignore-re = [
-    "ded",
+    "BA",
 ]
 
 [files]

From e269e36e422caa53763a01c09d90e3dd65047723 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 14:40:13 +0200
Subject: [PATCH 13/21] Use @init macro for initial guess in tutorial-symbolics

---
 docs/src/tutorial-symbolics.md | 9 ++++++---
 1 file changed, 6 insertions(+), 3 deletions(-)

diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index c9a52d6..6f2d6f2 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -242,7 +242,7 @@ We now formulate the optimal control problem using the `@def` macro from Optimal
 The boundary conditions ``X(0) = X(t_f)`` encode the periodicity of the orbit directly.
 
 ```@example main
-@def cartpole_ocp begin
+@def cartpole begin
     t ∈ [0, tf_val], time
     X = (x_ocp, θ_ocp, v_ocp, ω_ocp) ∈ R⁴, state
     F_ctrl ∈ R, control
@@ -263,9 +263,12 @@ grid of 100 intervals, then handed to **Ipopt** via NLPModelsIpopt. The `@time`
 reports wall-clock time, which on first call includes Julia's JIT compilation.
 
 ```@example main
-initial_guess = (state=[0.0, 0.0, 0.0, 0.1], control=0.0)
+initial_guess = @init cartpole begin
+    X(t) := [0.0, 0.0, 0.0, 0.1]
+    F_ctrl(t) := 0.0
+end
 
-sol = solve(cartpole_ocp; display=true, grid_size=100, init=initial_guess)
+sol = solve(cartpole; grid_size=100, init=initial_guess)
 ```
 
 ### Results

From 36cb010d4a7d76c651b39b6108d4f3ae3d0f3690 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 15:19:00 +0200
Subject: [PATCH 14/21] Increase example_size_threshold to 12000 to accommodate
 animation HTML

---
 docs/make.jl                   | 10 ++++++----
 docs/src/tutorial-symbolics.md |  4 ++--
 2 files changed, 8 insertions(+), 6 deletions(-)

diff --git a/docs/make.jl b/docs/make.jl
index 041ed0c..5301a9e 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -60,10 +60,12 @@ makedocs(;
     format=Documenter.HTML(;
         repolink="https://" * repo_url,
         prettyurls=false,
-        size_threshold_ignore=[
-            "tutorial-discretisation.md",
-            "tutorial-nlp.md",
-        ],
+        example_size_threshold=1_000_000,
+        size_threshold=1_000_000,
+        # size_threshold_ignore=[
+        #     "tutorial-discretisation.md",
+        #     "tutorial-nlp.md",
+        # ],
         assets=[
             asset("https://control-toolbox.org/assets/css/documentation.css"),
             asset("https://control-toolbox.org/assets/js/documentation.js"),
diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index 6f2d6f2..3e4004b 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -268,7 +268,7 @@ initial_guess = @init cartpole begin
     F_ctrl(t) := 0.0
 end
 
-sol = solve(cartpole; grid_size=100, init=initial_guess)
+sol = solve(cartpole; display=false, grid_size=100, init=initial_guess)
 ```
 
 ### Results
@@ -396,7 +396,7 @@ html_str = """
   };
   document.getElementById("cpSp$(uid)").oninput=function(){
     spd=parseFloat(this.value);
-    document.getElementById("cpSv$(uid)").textContent=spd.toFixed(2)+"\u00d7";
+    document.getElementById("cpSv$(uid)").textContent=spd.toFixed(2)+"×";
   };
   draw(0);
 })();

From 0be0465acc66c42a1bb940942aeb5bd2cdfdbb60 Mon Sep 17 00:00:00 2001
From: abavoil <17646791+abavoil@users.noreply.github.com>
Date: Wed, 15 Apr 2026 20:49:05 +0200
Subject: [PATCH 15/21] human revised version

---
 Project.toml                   |   4 -
 docs/src/assets/Manifest.toml  |   6 +-
 docs/src/tutorial-symbolics.md | 376 +++++++++++++++------------------
 3 files changed, 178 insertions(+), 208 deletions(-)

diff --git a/Project.toml b/Project.toml
index 9139201..26db0c8 100644
--- a/Project.toml
+++ b/Project.toml
@@ -3,9 +3,5 @@ uuid = "cb10daa6-a5e5-4c25-a171-ae181b8ea3c9"
 version = "0.4.1"
 authors = ["Olivier Cots <olivier.cots@toulouse-inp.fr>"]
 
-[deps]
-StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
-
 [compat]
-StaticArrays = "1.9.18"
 julia = "1.10"
diff --git a/docs/src/assets/Manifest.toml b/docs/src/assets/Manifest.toml
index ced7cf3..65c19ad 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.5"
+julia_version = "1.12.6"
 manifest_format = "2.0"
 project_hash = "1d73fc8b7306dbb42eb2a1fdd2ae0c8be50a6d19"
 
@@ -227,9 +227,9 @@ version = "0.18.7"
 
 [[deps.CTDirect]]
 deps = ["ADNLPModels", "CTModels", "CTSolvers", "DocStringExtensions", "ExaModels", "SolverCore", "SparseArrays", "SparseConnectivityTracer"]
-git-tree-sha1 = "a4d812f60412f47bf29e685a921d0d675f581d55"
+git-tree-sha1 = "db1699ffab2f31aeda5dbb61dfa9015773cb70fc"
 uuid = "790bbbee-bee9-49ee-8912-a9de031322d5"
-version = "1.0.7-beta"
+version = "1.0.6"
 
 [[deps.CTFlows]]
 deps = ["CTBase", "CTModels", "DocStringExtensions", "ForwardDiff", "LinearAlgebra", "MLStyle", "MacroTools"]
diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index 3e4004b..1842184 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -5,19 +5,17 @@ Draft = false
 ```
 
 This tutorial demonstrates how to combine **symbolic derivation of equations of motion** (via
-[Symbolics.jl](https://symbolics.juliasymbolics.org/)) with **direct optimal control**
+[`Symbolics.jl`](https://symbolics.juliasymbolics.org/)) with **direct optimal control**
 (via [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl/)) to find a
 **periodic orbit** (limit cycle) for a cart-pole system.
 
-The key idea is to let the computer do the hard mechanics: we write the Lagrangian in a few
-lines, and the Euler–Lagrange equations — including mass-matrix inversion — are derived
-automatically.
+The key is to define the kinetic and potential energies and the power of the non-conservative forces, and to let `Symbolics.jl` handle the derivations of the equations of motion using the Lagrange-Euler equation.
+
 
 ## The Cart-Pole System
 
 The system consists of a cart of mass ``m_c`` sliding on a frictionless horizontal rail, with
-a rigid pendulum of mass ``m_p`` and length ``l`` attached to it. A horizontal force ``u``
-(the control input) acts on the cart.
+a rigid pendulum of mass ``m_p`` and length ``l`` attached to it. A horizontal force ``u`` (the control input) acts on the cart.
 
 ```@raw html
 <figure>
@@ -26,7 +24,7 @@ a rigid pendulum of mass ``m_p`` and length ``l`` attached to it. A horizontal f
 </figure>
 ```
 
-The configuration vector is ``q = (x,\, \theta)^\top``, where ``x`` is the cart position and
+The configuration vector is ``q = (x,\, \theta)``, where ``x`` is the cart position and
 ``\theta = 0`` corresponds to the **upright** (unstable) equilibrium of the pendulum.
 
 ### Positions
@@ -53,25 +51,24 @@ T = \tfrac{1}{2}m_c\,\|\dot{p}_c\|^2 + \tfrac{1}{2}m_p\,\|\dot{p}_p\|^2
 V = m_p\,g\,l\cos\theta.
 ```
 
-The Lagrangian is ``\mathcal{L} = T - V``. The virtual work of the control force ``u``
-acting on the cart gives the generalised force vector:
+The Lagrangian is ``\mathcal{L} = T - V``. The virtual power ``P_{nc}`` of the non-conservative force (control) ``F = (u, 0)`` acting on the cart gives the generalised force vector:
 
 ```math
-W_\text{ctrl} = F \cdot \dot{p}_c = u\,\dot{x},
+P_{nc} = F \cdot \dot{p}_c = u\,\dot{x},
 \quad\Longrightarrow\quad
-Q = \frac{\partial W_\text{ctrl}}{\partial \dot{q}} = \begin{pmatrix} u \\ 0 \end{pmatrix}.
+Q = \frac{\partial P_{nc}}{\partial \dot{q}} = \begin{pmatrix} u \\ 0 \end{pmatrix}.
 ```
 
-### Euler–Lagrange Equations
+### Euler–Lagrange Equation
 
 The equations of motion follow from:
 
 ```math
-\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{q}_i}
-- \frac{\partial \mathcal{L}}{\partial q_i} = Q_i, \qquad i = 1,2.
+\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{q}}
+- \frac{\partial \mathcal{L}}{\partial q} = Q.
 ```
 
-They can be written in the standard **manipulator form**:
+They can be written in the standard **manipulator form** (also known as **robotic equation of motion**):
 
 ```math
 M(q)\,\ddot{q} = -C(q,\dot{q}) + \tau(q, u),
@@ -88,12 +85,12 @@ M(q) =
 ```
 
 and the right-hand side collects Coriolis/gravity terms and the control torque. Instead of
-deriving these by hand, we rely on Symbolics.jl to compute and **analytically invert**
+deriving these by hand, we rely on `Symbolics.jl` to compute and **analytically invert**
 ``M(q)`` for us.
 
 ### State-Space Form
 
-Defining the state ``X = (x,\,\theta,\,\dot{x},\,\dot{\theta})^\top``, the equations of
+Defining the state $X = (x,\,\theta,\,\dot{x},\,\dot{\theta})$, the equations of
 motion become the first-order system:
 
 ```math
@@ -102,52 +99,25 @@ motion become the first-order system:
   \dot{x} \\ \dot{\theta} \\ M^{-1}(q)\bigl(-C(q,\dot{q}) + \tau(q,u)\bigr)
 \end{pmatrix}.
 ```
+# Optimal Control of a Cart-Pole System using Symbolics.jl
 
-## The Optimal Control Problem
-
-We look for a **limit cycle** of the nonlinear dynamics: a trajectory that returns exactly to
-its initial condition after a fixed period ``t_f``. The cost penalises the total control
-energy:
-
-```math
-\min_{u(\cdot)}\; \int_0^{t_f} u(t)^2\,\mathrm{d}t
-```
-
-```math
-\text{subject to} \quad \dot{X}(t) = f(X(t), u(t)), \quad t \in [0, t_f],
-```
-
-```math
-X(0) = X(t_f) = (0,\, 0,\, 0,\, 0.1)^\top.
-```
-
-The common boundary condition ``X(0) = X(t_f)`` — combined with a non-trivial initial
-angular velocity — forces the solver to find an orbit rather than the trivial rest solution.
+This tutorial demonstrates how to use `Symbolics.jl` to automate the derivation of equations of motion (EOM) for a mechanical system and subsequently solve an optimal control problem using `OptimalControl.jl`.
 
 ## Implementation
 
 ### Setup & Imports
 
-```@setup main
-# Project environment for this tutorial — edit the path for your local setup.
-import Pkg
-Pkg.activate(".")
-```
-
 ```@example main
 using OptimalControl
-using Plots
-using StaticArrays
 using NLPModelsIpopt
 using Symbolics
+using LinearAlgebra: dot
+using Plots
 ```
 
 ### Physical Parameters and Symbolic Variables
 
-We declare all parameters both as numerical constants (for the final function evaluation) and
-as symbolic variables (for the Lagrangian computation). Note that the symbolic time variable
-`t` is used only inside the Symbolics.jl derivation and does not interfere with the time
-variable `t` introduced later by the `@def` macro.
+We declare all parameters both as numerical constants (for the final function evaluation) and as symbolic variables (for the Lagrangian computation). We define the configuration vector ``q = (x, \theta)``.
 
 ```@example main
 # Physical constants
@@ -162,110 +132,110 @@ const tf_val  = 2.0
 D = Differential(t)
 @variables m_c m_p l g u
 @variables x(t) θ(t)
-@variables v ω dv dω   # static aliases for velocities and accelerations
 
 q = [x, θ]
+
 nothing # hide
 ```
 
 ### Automated Kinematics and Lagrangian
 
-We express the positions, kinetic energy, potential energy, and virtual power symbolically.
-The time derivatives ``\dot{p}_c``, ``\dot{p}_p`` are computed automatically by `D.(...)`.
+We express the positions, kinetic energy, potential energy, and (non-conservative) power symbolically. The time derivatives ``\dot{p}_c`` and ``\dot{p}_p`` are computed automatically by `D.(...)`.
 
 ```@example main
 p_c = [x, 0.0]
 p_p = [x + l * sin(θ), l * cos(θ)]
-F   = [u, 0.0]
+F_ext = [u, 0.0]
+
+T = 0.5 * m_c * sum(abs2, D.(p_c)) + 0.5 * m_p * sum(abs2, D.(p_p))
+V = g * (m_p * p_p[2])
+L = T - V
 
-T      = 0.5 * m_c * sum(D.(p_c) .^ 2) + 0.5 * m_p * sum(D.(p_p) .^ 2)
-V      = g * (m_p * p_p[2])
-W_ctrl = transpose(D.(p_c)) * F   # virtual power of the control force
+P_non_conservative = dot(D.(p_c), F_ext)
 nothing # hide
 ```
 
 ### Euler–Lagrange Equations and Mass-Matrix Inversion
 
-Starting from ``\mathcal{L} = T - V``, Symbolics.jl computes the three terms of the
-Euler–Lagrange equations and assembles the residual vector. Substituting static aliases
-``(v,\omega,\dot{v},\dot\omega)`` for the time derivatives makes it possible to identify
-the mass matrix ``M`` as the Jacobian of the residual with respect to the accelerations
-``(\dot{v}, \dot\omega)``. The system ``M\,a = -b`` is then solved analytically.
+Starting from ``\mathcal{L} = T - V``, `Symbolics.jl` computes the terms of the Euler–Lagrange equations. To isolate the accelerations, we substitute the symbolic time derivatives with algebraic variables. This allows us to identify the **standard manipulator form** components: the mass matrix is the Jacobian of the residual with respect to the accelerations ``\ddot{q}``, and the bias vector contains the remaining terms.
 
 ```@example main
-L = T - V
-
 A = D.(Symbolics.gradient(L, D.(q)))
 B = Symbolics.gradient(L, q)
-Q = Symbolics.gradient(W_ctrl, D.(q))   # generalised forces
+Q = Symbolics.gradient(P_non_conservative, D.(q))
 
 # Euler-Lagrange residual: d/dt(∂L/∂q̇) - ∂L/∂q - Q = 0
 el_eqs = expand_derivatives.(A - B - Q)
 
+# Helper to create algebraic variables for derivatives
+diff_var(qi, suffix) = Symbolics.variable(Symbol(Symbolics.operation(Symbolics.value(qi)), suffix))
+
+# Static aliases for velocities and accelerations
+v = diff_var.(q, :_t)
+dv = diff_var.(q, :_tt)
+
 # Freeze time derivatives into static algebraic variables
-sub_rules = Dict(D(x) => v, D(θ) => ω, D(D(x)) => dv, D(D(θ)) => dω)
-res = Symbolics.substitute.(el_eqs, (sub_rules,))
+sub_rules = Dict([D.(q) .=> v; D.(D.(q)) .=> dv])
+residual = Symbolics.substitute.(el_eqs, (sub_rules,))
 
-# Identify mass matrix M and bias vector b such that M·[dv; dω] + b = 0
-Mass = Symbolics.jacobian(res, [dv, dω])
-bias = Symbolics.substitute.(res, (Dict(dv => 0.0, dω => 0.0),))
+# Identify mass matrix M and bias vector b such that M·dv + b = 0
+mass = Symbolics.jacobian(residual, dv)
+bias = Symbolics.substitute.(residual, (Dict(dv .=> 0.0),))
 
-# Analytically invert the 2×2 mass matrix
-accel = Symbolics.simplify_fractions.(Mass \ (-bias))
+# Solve for accelerations analytically: dv = M⁻¹(-b)
+accel = Symbolics.simplify_fractions.(mass \ (-bias))
 
-# Fully explicit state derivative: Ẋ = [v, ω, v̇, ω̇]
-dx_dt = [v, ω, accel[1], accel[2]]
+# Fully explicit state derivative: Ẋ = [v, accel]
+X = [q; v]
+dX = [v; accel]
 nothing # hide
 ```
 
 ### Code Generation
 
-`build_function` compiles the symbolic expression `dx_dt` into a native Julia function.
-The `force_SA=true` flag generates a **StaticArrays** kernel, which avoids heap allocations
-inside the ODE right-hand side — important for solver performance. Parameters are stored as
-an `SVector` to match the generated kernel's expected input type.
+`build_function` compiles the symbolic expression `dX` into a native Julia function with arguments `X`, `u`, and parameter values. The `force_SA=true` flag generates a **StaticArrays** kernel, which avoids heap allocations inside the ODE right-hand side — crucial for solver performance because dimension is small. For larger problems (``X \in \mathrm{R}^n``, ``n > 100``), we would use a mutating dynamics function instead, cf. [Julia Perfomance Tips](https://docs.julialang.org/en/v1/manual/performance-tips/#Consider-StaticArrays.jl-for-small-fixed-size-vector/matrix-operations).
 
 ```@example main
-f_expr = build_function(dx_dt, [x, θ, v, ω], u, [m_c, m_p, l, g];
+f_expr = build_function(dX, X, u, [m_c, m_p, l, g];
     expression=Val{false}, force_SA=true)
-f_cartpole = f_expr[1]   # out-of-place variant: (state, u, params) → SVector
+f_sym = f_expr[1]   # out-of-place variant: (state, u, params) → SVector
 
-const p_vals = SA[m_c_val, m_p_val, l_val, g_val]   # SVector, matches SA kernel
+const p_vals  = [m_c_val, m_p_val, l_val, g_val]
 
-cartpole_dynamics(X, U) = f_cartpole(X, U, p_vals)
+cartpole_dynamics(X, u) = f_sym(X, u, p_vals)
 nothing # hide
 ```
 
 ### Optimal Control Problem Definition
 
-We now formulate the optimal control problem using the `@def` macro from OptimalControl.jl.
-The boundary conditions ``X(0) = X(t_f)`` encode the periodicity of the orbit directly.
+We now formulate the optimal control problem using the `@def` macro from `OptimalControl.jl`. The initial state is not at rest because ``ω(0) = 0.2``, while the boundary condition ``X(0) - X(tf) = 0`` encodes the periodicity of the orbit.
 
 ```@example main
 @def cartpole begin
     t ∈ [0, tf_val], time
-    X = (x_ocp, θ_ocp, v_ocp, ω_ocp) ∈ R⁴, state
-    F_ctrl ∈ R, control
+    X = (x, θ, v, ω) ∈ R⁴, state
+    u ∈ R, control
 
-    X(0)      == [0, 0, 0, 0.1]
-    X(tf_val) == [0, 0, 0, 0.1]   # periodicity: X(tf) = X(0)
+    x(0) == 0
+    θ(0) == 0
+    v(0) == 0
+    ω(0) == 0.2
+    X(tf_val) - X(0) == [0, 0, 0, 0]  # Periodic orbit
 
-    Ẋ(t) == cartpole_dynamics(X(t), F_ctrl(t))
+    Ẋ(t) == cartpole_dynamics(X(t), u(t))
 
-    ∫(F_ctrl(t)^2) → min
+    ∫(u(t)^2) → min
 end
 ```
 
 ### Solving the NLP
 
-The problem is transcribed into a nonlinear program using direct collocation on a uniform
-grid of 100 intervals, then handed to **Ipopt** via NLPModelsIpopt. The `@time` macro
-reports wall-clock time, which on first call includes Julia's JIT compilation.
+The problem is transcribed into a nonlinear program using direct collocation on a uniform grid of 100 intervals, then handed to `Ipopt` via ``NLPModelsIpopt.jl``. We provide a simple initial guess for the state and control trajectories. See [the documentation](https://control-toolbox.org/OptimalControl.jl/stable/manual-solve.html) for more informations.
 
 ```@example main
 initial_guess = @init cartpole begin
-    X(t) := [0.0, 0.0, 0.0, 0.1]
-    F_ctrl(t) := 0.0
+    X(t) := [0.0, 0.0, 0.0, 0.2]
+    u(t) := 0.0
 end
 
 sol = solve(cartpole; display=false, grid_size=100, init=initial_guess)
@@ -274,136 +244,140 @@ sol = solve(cartpole; display=false, grid_size=100, init=initial_guess)
 ### Results
 
 ```@example main
-println("--- Optimal Limit Cycle Found ---")
-println("Total Control Energy (∫F² dt): ", sol.objective)
-
 tsol = time_grid(sol)
 Xsol = state(sol).(tsol)
-Fsol = control(sol).(tsol)
+usol = control(sol).(tsol)
 
-# Stack state vectors into a 4×N matrix and unpack rows
-xsol, θsol, vsol, ωsol = eachrow(reduce(hcat, Xsol))
+X_mat = reduce(hcat, Xsol)
+q_sol = X_mat[1:2, :]'
+dq_sol = X_mat[3:4, :]'
 
-qplot  = plot(tsol, [xsol θsol], label=["x" "θ"],   title="Configuration")
-dqplot = plot(tsol, [vsol ωsol], label=["v" "ω"],   title="Velocities")
-uplot  = plot(tsol, Fsol,        label="u",          title="Control", linetype=:steppost)
+p1 = plot(tsol, q_sol, label=["x" "θ"], title="Configuration")
+p2 = plot(tsol, dq_sol, label=["v" "ω"], title="Velocities")
+p3 = plot(tsol, usol, label="u", title="Control", linetype=:steppost)
 
-plot(qplot, dqplot, uplot, layout=3, size=(800, 600))
+plot(p1, p2, p3, layout=(3, 1), size=(800, 700))
 ```
 
-The three panels show the cart position ``x`` and pendulum angle ``\theta``, the
-corresponding velocities ``\dot{x}`` and ``\dot\theta``, and the optimal control force
-``u``. The periodicity of the state trajectory confirms that a genuine limit cycle has been
-found.
+The plots show the cart position ``x`` and pendulum angle ``\theta``, the corresponding velocities, and the optimal control force ``u`` required to stabilize the system back to its initial state within 2 seconds.
+
 
 ### Animation
 
-The animation below shows the cart-pole evolving along the optimal limit-cycle trajectory.
-The blue cart slides on the horizontal rail while the pendulum swings around the upright
-equilibrium. A ghost trace follows the red bob to make the motion easier to read. Use the
-controls to play, pause, or reset the animation and to adjust the playback speed.
+The animation below shows the cart-pole evolving along the optimal limit-cycle trajectory. The blue cart slides on the horizontal rail while the pendulum swings around the upright equilibrium.
 
 ```@setup main
-# Serialise trajectory for the JS animation (manual, no extra dependencies)
+# Extract states for the animation
+xsol = X_mat[1, :]
+θsol = X_mat[2, :]
+
+# Serialise trajectory for JS (manual, no extra dependencies)
 json_t  = "[" * join(string.(tsol),  ",") * "]"
 json_x  = "[" * join(string.(xsol),  ",") * "]"
 json_th = "[" * join(string.(θsol),  ",") * "]"
 
-uid    = string(rand(UInt32), base=16)   # unique suffix → safe for multi-example pages
-tf_str = string(round(tsol[end], digits=3))
-
-html_str = """
-<div style="text-align:center;font-family:sans-serif;margin:1.5em 0;">
-  <canvas id="cpC$(uid)" width="700" height="380"
-    style="border:1px solid #ccc;border-radius:6px;background:#fafafa;display:block;margin:0 auto;"></canvas>
-  <div style="margin:6px 0;font-size:0.9em;color:#444;">
-    <span id="cpT$(uid)">t = 0.000 / $(tf_str) s</span>
-  </div>
-  <div style="display:flex;justify-content:center;gap:10px;align-items:center;flex-wrap:wrap;margin-top:4px;">
-    <button id="cpPl$(uid)" style="padding:5px 14px;cursor:pointer;">&#9654; Play</button>
-    <button id="cpPa$(uid)" style="padding:5px 14px;cursor:pointer;">&#9646;&#9646; Pause</button>
-    <button id="cpRe$(uid)" style="padding:5px 14px;cursor:pointer;">&#x21BA; Reset</button>
-    <label style="font-size:0.9em;">Speed:
-      <input type="range" id="cpSp$(uid)" min="0.25" max="3" step="0.25" value="1"
-             style="vertical-align:middle;width:100px;">
-      <span id="cpSv$(uid)">1.00&times;</span>
-    </label>
-  </div>
+# Define a wrapper to instruct Documenter to render the output as HTML
+struct RawHTML
+    raw::String
+end
+Base.show(io::IO, ::MIME"text/html", h::RawHTML) = print(io, h.raw)
+
+html_anim = """
+<div style="display: flex; justify-content: center; margin: 20px 0;">
+    <canvas id="cartpoleCanvas" width="600" height="300" 
+            style="border:1px solid #ddd; background:#fafafa; border-radius: 8px; overflow: hidden;">
+    </canvas>
 </div>
+
 <script>
-(function(){
-  var T=$(json_t), X=$(json_x), TH=$(json_th), L=$(l_val);
-  var cv=document.getElementById("cpC$(uid)");
-  var ctx=cv.getContext("2d");
-  var W=cv.width, H=cv.height;
-  var railY=H*0.62, sc=80, cW=60, cH=28;
-  var tSim=0, spd=1, rid=null, lts=null, ghost=[];
-  function lerp(a,b,t){return a+(b-a)*t;}
-  function interp(tq){
-    var n=T.length;
-    if(tq<=T[0])return{x:X[0],th:TH[0]};
-    if(tq>=T[n-1])return{x:X[n-1],th:TH[n-1]};
-    var lo=0,hi=n-2;
-    while(lo<hi){var m=(lo+hi)>>1;if(T[m+1]<tq)lo=m+1;else hi=m;}
-    var a=(tq-T[lo])/(T[lo+1]-T[lo]);
-    return{x:lerp(X[lo],X[lo+1],a),th:lerp(TH[lo],TH[lo+1],a)};
-  }
-  function draw(tq){
-    var d=interp(tq);
-    var cx=W/2+d.x*sc, cy=railY;
-    var px=cx, py=cy-cH/2;
-    var bx=px+L*sc*Math.sin(d.th), by=py-L*sc*Math.cos(d.th);
-    ghost.push({x:bx,y:by}); if(ghost.length>25)ghost.shift();
-    ctx.clearRect(0,0,W,H);
-    ctx.beginPath();ctx.moveTo(20,railY);ctx.lineTo(W-20,railY);
-    ctx.strokeStyle="#aaa";ctx.lineWidth=3;ctx.stroke();
-    var wr=8;
-    [[cx-18,cy+cH/2+wr],[cx+18,cy+cH/2+wr]].forEach(function(p){
-      ctx.beginPath();ctx.arc(p[0],p[1],wr,0,2*Math.PI);
-      ctx.fillStyle="#888";ctx.fill();
-    });
-    ctx.beginPath();ctx.rect(cx-cW/2,cy-cH/2,cW,cH);
-    ctx.fillStyle="#4a90d9";ctx.fill();
-    ctx.strokeStyle="#2c5f8a";ctx.lineWidth=1.5;ctx.stroke();
-    ghost.forEach(function(g,i){
-      ctx.beginPath();ctx.arc(g.x,g.y,5,0,2*Math.PI);
-      ctx.fillStyle="rgba(220,60,60,"+((i+1)/ghost.length*0.35)+")";ctx.fill();
-    });
-    ctx.beginPath();ctx.moveTo(px,py);ctx.lineTo(bx,by);
-    ctx.strokeStyle="#222";ctx.lineWidth=3;ctx.stroke();
-    ctx.beginPath();ctx.arc(px,py,5,0,2*Math.PI);ctx.fillStyle="#555";ctx.fill();
-    ctx.beginPath();ctx.arc(bx,by,10,0,2*Math.PI);
-    ctx.fillStyle="#dc3c3c";ctx.fill();ctx.strokeStyle="#8a1a1a";ctx.lineWidth=1.5;ctx.stroke();
-    document.getElementById("cpT$(uid)").textContent=
-      "t = "+tq.toFixed(3)+" / "+T[T.length-1].toFixed(2)+" s";
-  }
-  function step(ts){
-    if(lts!==null){
-      var dt=(ts-lts)*0.001*spd; tSim+=dt;
-      var tf=T[T.length-1]; if(tSim>tf){tSim-=tf;ghost=[];}
+(function() {
+    // Julia directly interpolates the arrays into the JS here
+    const t = $json_t;
+    const x = $json_x;
+    const th = $json_th;
+    
+    const canvas = document.getElementById('cartpoleCanvas');
+    const ctx = canvas.getContext('2d');
+    
+    const scale = 50; 
+    const l_px = $(l_val) * scale; 
+    const duration = t[t.length - 1];
+
+    let start_time = null;
+
+    function draw(time) {
+        if (!start_time) start_time = time;
+        
+        // Seamless periodic loop
+        let elapsed = time - start_time;
+        let sim_t = (elapsed / 1000.0) % duration;
+        
+        // Find the current time interval
+        let i = 0;
+        while(i < t.length - 1 && t[i + 1] < sim_t) { i++; }
+
+        // Linear interpolation for buttery smooth frames
+        let dt = t[i+1] - t[i];
+        let alpha = (dt > 0) ? (sim_t - t[i]) / dt : 0;
+        let cur_x = x[i] + alpha * (x[i+1] - x[i]);
+        let cur_th = th[i] + alpha * (th[i+1] - th[i]);
+
+        ctx.clearRect(0, 0, canvas.width, canvas.height);
+        
+        const cx = canvas.width / 2 + cur_x * scale;
+        const cy = canvas.height / 2 + 40; 
+
+        // Draw Track
+        ctx.beginPath();
+        ctx.moveTo(0, cy + 15);
+        ctx.lineTo(canvas.width, cy + 15);
+        ctx.lineWidth = 2;
+        ctx.strokeStyle = '#ccc';
+        ctx.stroke();
+
+        // Draw Cart
+        ctx.fillStyle = '#3498db';
+        ctx.fillRect(cx - 30, cy - 15, 60, 30);
+        
+        // Draw Pole
+        const px = cx + l_px * Math.sin(cur_th);
+        const py = cy - l_px * Math.cos(cur_th);
+        
+        ctx.beginPath();
+        ctx.moveTo(cx, cy);
+        ctx.lineTo(px, py);
+        ctx.lineWidth = 6;
+        ctx.lineCap = 'round';
+        ctx.strokeStyle = '#e74c3c';
+        ctx.stroke();
+        
+        // Draw Hinge
+        ctx.beginPath();
+        ctx.arc(cx, cy, 6, 0, 2 * Math.PI);
+        ctx.fillStyle = '#2c3e50';
+        ctx.fill();
+
+        // Draw Bob
+        ctx.beginPath();
+        ctx.arc(px, py, 8, 0, 2 * Math.PI);
+        ctx.fillStyle = '#e74c3c';
+        ctx.fill();
+
+        // Draw Progress Bar at the very bottom
+        const bar_height = 4;
+        ctx.fillStyle = '#ecf0f1'; // Background
+        ctx.fillRect(0, canvas.height - bar_height, canvas.width, bar_height);
+        ctx.fillStyle = '#3498db'; // Active progress
+        ctx.fillRect(0, canvas.height - bar_height, canvas.width * (sim_t / duration), bar_height);
+
+        requestAnimationFrame(draw);
     }
-    lts=ts; draw(tSim); rid=requestAnimationFrame(step);
-  }
-  document.getElementById("cpPl$(uid)").onclick=function(){
-    if(rid===null){lts=null;rid=requestAnimationFrame(step);}
-  };
-  document.getElementById("cpPa$(uid)").onclick=function(){
-    if(rid!==null){cancelAnimationFrame(rid);rid=null;}
-  };
-  document.getElementById("cpRe$(uid)").onclick=function(){
-    if(rid!==null){cancelAnimationFrame(rid);rid=null;}
-    tSim=0;ghost=[];lts=null;draw(0);
-  };
-  document.getElementById("cpSp$(uid)").oninput=function(){
-    spd=parseFloat(this.value);
-    document.getElementById("cpSv$(uid)").textContent=spd.toFixed(2)+"×";
-  };
-  draw(0);
+    requestAnimationFrame(draw);
 })();
 </script>
 """
 ```
 
 ```@example main
-HTML(html_str) # hide
+RawHTML(html_anim)  # hide
 ```

From b6b8ab2ac7a711ff4b48796247e2520a9f863e0e Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 21:08:36 +0200
Subject: [PATCH 16/21] Increase size thresholds and fix typo in
 tutorial-symbolics

---
 docs/make.jl                   | 1 +
 docs/src/assets/Manifest.toml  | 6 +++---
 docs/src/tutorial-symbolics.md | 3 +--
 3 files changed, 5 insertions(+), 5 deletions(-)

diff --git a/docs/make.jl b/docs/make.jl
index 5301a9e..ff40741 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -61,6 +61,7 @@ makedocs(;
         repolink="https://" * repo_url,
         prettyurls=false,
         example_size_threshold=1_000_000,
+        size_threshold_warn=1_000_000,
         size_threshold=1_000_000,
         # size_threshold_ignore=[
         #     "tutorial-discretisation.md",
diff --git a/docs/src/assets/Manifest.toml b/docs/src/assets/Manifest.toml
index 65c19ad..ced7cf3 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.6"
+julia_version = "1.12.5"
 manifest_format = "2.0"
 project_hash = "1d73fc8b7306dbb42eb2a1fdd2ae0c8be50a6d19"
 
@@ -227,9 +227,9 @@ version = "0.18.7"
 
 [[deps.CTDirect]]
 deps = ["ADNLPModels", "CTModels", "CTSolvers", "DocStringExtensions", "ExaModels", "SolverCore", "SparseArrays", "SparseConnectivityTracer"]
-git-tree-sha1 = "db1699ffab2f31aeda5dbb61dfa9015773cb70fc"
+git-tree-sha1 = "a4d812f60412f47bf29e685a921d0d675f581d55"
 uuid = "790bbbee-bee9-49ee-8912-a9de031322d5"
-version = "1.0.6"
+version = "1.0.7-beta"
 
 [[deps.CTFlows]]
 deps = ["CTBase", "CTModels", "DocStringExtensions", "ForwardDiff", "LinearAlgebra", "MLStyle", "MacroTools"]
diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index 1842184..27ac627 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -230,7 +230,7 @@ end
 
 ### Solving the NLP
 
-The problem is transcribed into a nonlinear program using direct collocation on a uniform grid of 100 intervals, then handed to `Ipopt` via ``NLPModelsIpopt.jl``. We provide a simple initial guess for the state and control trajectories. See [the documentation](https://control-toolbox.org/OptimalControl.jl/stable/manual-solve.html) for more informations.
+The problem is transcribed into a nonlinear program using direct collocation on a uniform grid of 100 intervals, then handed to `Ipopt` via `NLPModelsIpopt.jl`. We provide a simple initial guess for the state and control trajectories. See [the documentation](@extref OptimalControl manual-solve) for more information.
 
 ```@example main
 initial_guess = @init cartpole begin
@@ -261,7 +261,6 @@ plot(p1, p2, p3, layout=(3, 1), size=(800, 700))
 
 The plots show the cart position ``x`` and pendulum angle ``\theta``, the corresponding velocities, and the optimal control force ``u`` required to stabilize the system back to its initial state within 2 seconds.
 
-
 ### Animation
 
 The animation below shows the cart-pole evolving along the optimal limit-cycle trajectory. The blue cart slides on the horizontal rail while the pendulum swings around the upright equilibrium.

From 0d0d449b13ddd90dc17139085e27b2b20c5c6ebb Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 22:10:46 +0200
Subject: [PATCH 17/21] foo

---
 docs/src/tutorial-symbolics.md | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index 27ac627..2c457cc 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -193,7 +193,7 @@ nothing # hide
 
 ### Code Generation
 
-`build_function` compiles the symbolic expression `dX` into a native Julia function with arguments `X`, `u`, and parameter values. The `force_SA=true` flag generates a **StaticArrays** kernel, which avoids heap allocations inside the ODE right-hand side — crucial for solver performance because dimension is small. For larger problems (``X \in \mathrm{R}^n``, ``n > 100``), we would use a mutating dynamics function instead, cf. [Julia Perfomance Tips](https://docs.julialang.org/en/v1/manual/performance-tips/#Consider-StaticArrays.jl-for-small-fixed-size-vector/matrix-operations).
+`build_function` compiles the symbolic expression `dX` into a native Julia function with arguments `X`, `u`, and parameter values. The `force_SA=true` flag generates a **StaticArrays** kernel, which avoids heap allocations inside the ODE right-hand side — crucial for solver performance because dimension is small. For larger problems (``X \in \mathrm{R}^n``, ``n > 100``), we would use a mutating dynamics function instead, cf. [Julia Performance Tips](https://docs.julialang.org/en/v1/manual/performance-tips/#Consider-StaticArrays.jl-for-small-fixed-size-vector/matrix-operations).
 
 ```@example main
 f_expr = build_function(dX, X, u, [m_c, m_p, l, g];
@@ -284,7 +284,7 @@ Base.show(io::IO, ::MIME"text/html", h::RawHTML) = print(io, h.raw)
 html_anim = """
 <div style="display: flex; justify-content: center; margin: 20px 0;">
     <canvas id="cartpoleCanvas" width="600" height="300" 
-            style="border:1px solid #ddd; background:#fafafa; border-radius: 8px; overflow: hidden;">
+            style="border:1px solid #ddd; background:#fafafa; border-radius: 8px; max-width: 100%;">
     </canvas>
 </div>
 

From aff54a860f94d747220ca827cbffc8b806ea68fd Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 22:42:17 +0200
Subject: [PATCH 18/21] Add dark/light mode to cartpole SVG and animation, use
 Julia colors

---
 docs/src/assets/cartpole_sketch.svg | 237 ++++++++++++++++++++++++----
 docs/src/tutorial-symbolics.md      |  45 ++++--
 2 files changed, 237 insertions(+), 45 deletions(-)

diff --git a/docs/src/assets/cartpole_sketch.svg b/docs/src/assets/cartpole_sketch.svg
index ea4addb..d9e080f 100644
--- a/docs/src/assets/cartpole_sketch.svg
+++ b/docs/src/assets/cartpole_sketch.svg
@@ -1,39 +1,206 @@
-<svg width="350" height="250" viewBox="0 0 350 250" xmlns="http://www.w3.org/2000/svg">
-  <!-- Background -->
-  <rect width="100%" height="100%" fill="white"/>
-  
-  <!-- Rail -->
-  <line x1="30" y1="180" x2="320" y2="180" stroke="#666" stroke-width="3"/>
-  
-  <!-- Cart wheels -->
-  <circle cx="130" cy="188" r="8" fill="#888"/>
-  <circle cx="170" cy="188" r="8" fill="#888"/>
-  
-  <!-- Cart body -->
-  <rect x="120" y="160" width="60" height="25" rx="3" fill="#4a90d9" stroke="#2c5f8a" stroke-width="2"/>
-  
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<svg
+   width="313.82666"
+   height="158.66362"
+   viewBox="0 0 313.82666 158.66362"
+   version="1.1"
+   id="svg216"
+   sodipodi:docname="cartpole_sketch.svg"
+   inkscape:version="1.2 (dc2aeda, 2022-05-15)"
+   xmlns:inkscape="http://www.inkscape.org/namespaces/inkscape"
+   xmlns:sodipodi="http://sodipodi.sourceforge.net/DTD/sodipodi-0.dtd"
+   xmlns="http://www.w3.org/2000/svg"
+   xmlns:svg="http://www.w3.org/2000/svg">
+  <sodipodi:namedview
+     id="namedview218"
+     pagecolor="#ffffff"
+     bordercolor="#000000"
+     borderopacity="0.25"
+     inkscape:showpageshadow="2"
+     inkscape:pageopacity="0.0"
+     inkscape:pagecheckerboard="0"
+     inkscape:deskcolor="#d1d1d1"
+     showgrid="false"
+     inkscape:zoom="1.6047635"
+     inkscape:cx="123.69424"
+     inkscape:cy="61.068188"
+     inkscape:window-width="1512"
+     inkscape:window-height="916"
+     inkscape:window-x="0"
+     inkscape:window-y="38"
+     inkscape:window-maximized="1"
+     inkscape:current-layer="svg216" />
+  <defs
+     id="defs190">
+    <marker
+       style="overflow:visible"
+       id="TriangleStart"
+       refX="0"
+       refY="0"
+       orient="auto-start-reverse"
+       inkscape:stockid="TriangleStart"
+       markerWidth="6"
+       markerHeight="5"
+       viewBox="0 0 5.3244081 6.1553851"
+       inkscape:isstock="true"
+       inkscape:collect="always"
+       preserveAspectRatio="none">
+      <path
+         transform="scale(0.5)"
+         style="fill:context-stroke;fill-rule:evenodd;stroke:context-stroke;stroke-width:1pt"
+         d="M 5.77,0 -2.88,5 V -5 Z"
+         id="path135" />
+    </marker>
+    <style
+       id="style182">
+      .track   { stroke: #999; }
+      .cart    { fill: #4063D8; stroke: #2c4fa0; }
+      .hinge   { fill: #2c3e50; }
+      .rod     { stroke: #555; }
+      .bob     { fill: #CB3C33; stroke: #8a1a1a; }
+      .annot   { stroke: #888; fill: none; }
+      .label   { fill: #555; font-family: Arial, sans-serif; }
+      .force   { stroke: #CB3C33; }
+      .arrow   { fill: #CB3C33; }
+      .xarrow  { stroke: #555; }
+      .xlabel  { fill: #555; font-family: Arial, sans-serif; font-size: 13px; }
+      @media (prefers-color-scheme: dark) {
+        .track  { stroke: #555; }
+        .cart   { fill: #4063D8; stroke: #6080e8; }
+        .hinge  { fill: #ecf0f1; }
+        .rod    { stroke: #bdc3c7; }
+        .bob    { fill: #CB3C33; stroke: #e05050; }
+        .annot  { stroke: #aaa; fill: none; }
+        .label  { fill: #bdc3c7; }
+        .force  { stroke: #e05050; }
+        .arrow  { fill: #e05050; }
+        .xarrow { stroke: #bdc3c7; }
+        .xlabel { fill: #bdc3c7; }
+      }
+    </style>
+    <marker
+       id="arrowhead"
+       markerWidth="10"
+       markerHeight="7"
+       refX="9"
+       refY="3.5"
+       orient="auto"
+       viewBox="0 0 10 7"
+       preserveAspectRatio="xMidYMid">
+      <polygon
+         points="10,3.5 0,7 0,0 "
+         class="arrow"
+         id="polygon184" />
+    </marker>
+    <marker
+       id="xarrowhead"
+       markerWidth="8"
+       markerHeight="6"
+       refX="7"
+       refY="3"
+       orient="auto"
+       viewBox="0 0 8 6"
+       preserveAspectRatio="xMidYMid">
+      <polygon
+         points="0,0 8,3 0,6 "
+         class="arrow"
+         style="fill:#555555"
+         id="polygon187" />
+    </marker>
+  </defs>
+  <path
+     id="line192-0"
+     style="fill:none;stroke:#999999;stroke-width:0.749806"
+     class="track"
+     d="m 109.01617,66.335535 c 0,0 2.66987,-7.02161 10.83691,-9.55147 8.16704,-2.52986 13.65788,1.23699 13.65788,1.23699"
+     sodipodi:nodetypes="czc" />
+  <line
+     x1="133.40085"
+     y1="99.381401"
+     x2="82.760437"
+     y2="31.113207"
+     class="rod"
+     stroke-width="3"
+     stroke-linecap="round"
+     id="line198" />
+  <!-- Track -->
+  <line
+     x1="13.590039"
+     y1="125.55579"
+     x2="303.59003"
+     y2="125.55579"
+     class="track"
+     stroke-width="2.5"
+     id="line192" />
+  <!-- Cart body (no wheels) -->
+  <rect
+     x="101.59004"
+     y="100.55579"
+     width="64"
+     height="26"
+     rx="4"
+     class="cart"
+     stroke-width="1.5"
+     id="rect194" />
+  <!-- Hinge -->
+  <circle
+     cx="133.59004"
+     cy="100.55579"
+     r="5"
+     class="hinge"
+     id="circle196" />
   <!-- Pendulum rod -->
-  <line x1="150" y1="160" x2="150" y2="80" stroke="#333" stroke-width="3"/>
-  
   <!-- Pendulum bob -->
-  <circle cx="150" cy="80" r="12" fill="#dc3c3c" stroke="#8a1a1a" stroke-width="2"/>
-  
-  <!-- Angle θ annotation -->
-  <path d="M 150 160 L 150 140" stroke="#666" stroke-width="1" stroke-dasharray="3,3"/>
-  <text x="155" y="145" font-family="Arial, sans-serif" font-size="14" fill="#666">θ</text>
-  
+  <circle
+     cx="47.933208"
+     cy="74.295029"
+     r="12"
+     class="bob"
+     stroke-width="1.5"
+     id="circle200"
+     transform="rotate(-36.567523)" />
+  <!-- Angle θ dashed reference + arc -->
+  <line
+     x1="133.59004"
+     y1="100.55579"
+     x2="133.59004"
+     y2="12.343751"
+     class="annot"
+     stroke-width="1.66031"
+     stroke-dasharray="4, 3"
+     id="line202" />
+  <text
+     x="114.16084"
+     y="52.973534"
+     class="label"
+     font-size="14px"
+     font-style="italic"
+     id="text206">θ</text>
   <!-- Force u arrow -->
-  <line x1="180" y1="172" x2="220" y2="172" stroke="#e74c3c" stroke-width="2" marker-end="url(#arrowhead)"/>
-  <text x="200" y="165" font-family="Arial, sans-serif" font-size="12" fill="#e74c3c" text-anchor="middle">u</text>
-  
-  <!-- Arrow marker definition -->
-  <defs>
-    <marker id="arrowhead" markerWidth="10" markerHeight="7" refX="9" refY="3.5" orient="auto">
-      <polygon points="0 0, 10 3.5, 0 7" fill="#e74c3c"/>
-    </marker>
-  </defs>
-  
-  <!-- Labels -->
-  <text x="150" y="220" font-family="Arial, sans-serif" font-size="12" fill="#333" text-anchor="middle">x</text>
-  <text x="150" y="240" font-family="Arial, sans-serif" font-size="12" fill="#666" text-anchor="middle">rail</text>
+  <line
+     x1="165.59004"
+     y1="113.55579"
+     x2="225.27538"
+     y2="113.55579"
+     class="force"
+     stroke-width="2.44306"
+     marker-end="url(#arrowhead)"
+     id="line208"
+     style="stroke-width:2.443;stroke-dasharray:none;marker-end:url(#TriangleStart)" />
+  <text
+     x="193.18916"
+     y="107.16372"
+     class="label"
+     font-size="12px"
+     text-anchor="middle"
+     font-style="italic"
+     id="text210">u</text>
+  <!-- x position arrow + label -->
+  <text
+     x="133.78026"
+     y="145.27246"
+     class="xlabel"
+     text-anchor="middle"
+     font-style="italic"
+     id="text214">x</text>
 </svg>
diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index 2c457cc..2fd3ba4 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -11,7 +11,6 @@ This tutorial demonstrates how to combine **symbolic derivation of equations of
 
 The key is to define the kinetic and potential energies and the power of the non-conservative forces, and to let `Symbolics.jl` handle the derivations of the equations of motion using the Lagrange-Euler equation.
 
-
 ## The Cart-Pole System
 
 The system consists of a cart of mass ``m_c`` sliding on a frictionless horizontal rail, with
@@ -99,6 +98,7 @@ motion become the first-order system:
   \dot{x} \\ \dot{\theta} \\ M^{-1}(q)\bigl(-C(q,\dot{q}) + \tau(q,u)\bigr)
 \end{pmatrix}.
 ```
+
 # Optimal Control of a Cart-Pole System using Symbolics.jl
 
 This tutorial demonstrates how to use `Symbolics.jl` to automate the derivation of equations of motion (EOM) for a mechanical system and subsequently solve an optimal control problem using `OptimalControl.jl`.
@@ -284,7 +284,7 @@ Base.show(io::IO, ::MIME"text/html", h::RawHTML) = print(io, h.raw)
 html_anim = """
 <div style="display: flex; justify-content: center; margin: 20px 0;">
     <canvas id="cartpoleCanvas" width="600" height="300" 
-            style="border:1px solid #ddd; background:#fafafa; border-radius: 8px; max-width: 100%;">
+            style="border:1px solid #ddd; border-radius: 8px; max-width: 100%;">
     </canvas>
 </div>
 
@@ -298,6 +298,30 @@ html_anim = """
     const canvas = document.getElementById('cartpoleCanvas');
     const ctx = canvas.getContext('2d');
     
+    // Detect dark/light mode
+    const isDark = window.matchMedia && window.matchMedia('(prefers-color-scheme: dark)').matches;
+    
+    // Color palettes
+    const colors = isDark ? {
+        canvas_bg: '#1e1e2e',
+        track: '#4a4a4a',
+        cart: '#4063D8',
+        pole: '#CB3C33',
+        hinge: '#ecf0f1',
+        bob: '#CB3C33',
+        progress_bg: '#3a3a3a',
+        progress: '#4063D8'
+    } : {
+        canvas_bg: '#fafafa',
+        track: '#999',
+        cart: '#4063D8',
+        pole: '#555',
+        hinge: '#2c3e50',
+        bob: '#CB3C33',
+        progress_bg: '#ecf0f1',
+        progress: '#4063D8'
+    };
+    
     const scale = 50; 
     const l_px = $(l_val) * scale; 
     const duration = t[t.length - 1];
@@ -321,7 +345,8 @@ html_anim = """
         let cur_x = x[i] + alpha * (x[i+1] - x[i]);
         let cur_th = th[i] + alpha * (th[i+1] - th[i]);
 
-        ctx.clearRect(0, 0, canvas.width, canvas.height);
+        ctx.fillStyle = colors.canvas_bg;
+        ctx.fillRect(0, 0, canvas.width, canvas.height);
         
         const cx = canvas.width / 2 + cur_x * scale;
         const cy = canvas.height / 2 + 40; 
@@ -331,11 +356,11 @@ html_anim = """
         ctx.moveTo(0, cy + 15);
         ctx.lineTo(canvas.width, cy + 15);
         ctx.lineWidth = 2;
-        ctx.strokeStyle = '#ccc';
+        ctx.strokeStyle = colors.track;
         ctx.stroke();
 
         // Draw Cart
-        ctx.fillStyle = '#3498db';
+        ctx.fillStyle = colors.cart;
         ctx.fillRect(cx - 30, cy - 15, 60, 30);
         
         // Draw Pole
@@ -347,26 +372,26 @@ html_anim = """
         ctx.lineTo(px, py);
         ctx.lineWidth = 6;
         ctx.lineCap = 'round';
-        ctx.strokeStyle = '#e74c3c';
+        ctx.strokeStyle = colors.pole;
         ctx.stroke();
         
         // Draw Hinge
         ctx.beginPath();
         ctx.arc(cx, cy, 6, 0, 2 * Math.PI);
-        ctx.fillStyle = '#2c3e50';
+        ctx.fillStyle = colors.hinge;
         ctx.fill();
 
         // Draw Bob
         ctx.beginPath();
         ctx.arc(px, py, 8, 0, 2 * Math.PI);
-        ctx.fillStyle = '#e74c3c';
+        ctx.fillStyle = colors.bob;
         ctx.fill();
 
         // Draw Progress Bar at the very bottom
         const bar_height = 4;
-        ctx.fillStyle = '#ecf0f1'; // Background
+        ctx.fillStyle = colors.progress_bg;
         ctx.fillRect(0, canvas.height - bar_height, canvas.width, bar_height);
-        ctx.fillStyle = '#3498db'; // Active progress
+        ctx.fillStyle = colors.progress; // Active progress
         ctx.fillRect(0, canvas.height - bar_height, canvas.width * (sim_t / duration), bar_height);
 
         requestAnimationFrame(draw);

From d150b47faaaf4c1f7ed85ed2651e03f799096d7b Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Wed, 15 Apr 2026 22:44:18 +0200
Subject: [PATCH 19/21] Increase cartpole canvas width to 900

---
 docs/src/tutorial-symbolics.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index 2fd3ba4..12471c4 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -283,7 +283,7 @@ Base.show(io::IO, ::MIME"text/html", h::RawHTML) = print(io, h.raw)
 
 html_anim = """
 <div style="display: flex; justify-content: center; margin: 20px 0;">
-    <canvas id="cartpoleCanvas" width="600" height="300" 
+    <canvas id="cartpoleCanvas" width="900" height="300" 
             style="border:1px solid #ddd; border-radius: 8px; max-width: 100%;">
     </canvas>
 </div>

From 368d085b8b2ebe4b4d29a3eea2f894bb81404723 Mon Sep 17 00:00:00 2001
From: Olivier Cots <olivier.cots@irit.fr>
Date: Thu, 16 Apr 2026 09:34:45 +0200
Subject: [PATCH 20/21] Merge main changes into symbolics branch

---
 docs/make.jl                        |   2 +-
 docs/src/tutorial-continuation.md   |   4 -
 docs/src/tutorial-symbolics-save.md | 313 ----------------------------
 docs/src/tutorial-symbolics.md      |   4 -
 4 files changed, 1 insertion(+), 322 deletions(-)
 delete mode 100644 docs/src/tutorial-symbolics-save.md

diff --git a/docs/make.jl b/docs/make.jl
index 989c5bb..aad4b85 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -48,7 +48,7 @@ repo_url = "github.com/control-toolbox/Tutorials.jl"
 Draft = false
 ```
 =#
-draft = true  # Draft mode: if true, @example blocks in markdown are not executed
+draft = false  # Draft mode: if true, @example blocks in markdown are not executed
 
 # ═══════════════════════════════════════════════════════════════════════════════
 # Build documentation
diff --git a/docs/src/tutorial-continuation.md b/docs/src/tutorial-continuation.md
index acabcb7..b58d1aa 100644
--- a/docs/src/tutorial-continuation.md
+++ b/docs/src/tutorial-continuation.md
@@ -1,9 +1,5 @@
 # [Discrete continuation](@id tutorial-continuation)
 
-```@meta
-Draft = false
-```
-
 By using the warm start option, it is easy to implement a basic discrete continuation method, in which a sequence of problems is solved by using each solution as the initial guess for the next problem. This approach typically leads to faster and more reliable convergence than solving each problem with the same initial guess and is particularly useful for problems that require a good initial guess to converge.
 
 ## Continuation on parametric OCP
diff --git a/docs/src/tutorial-symbolics-save.md b/docs/src/tutorial-symbolics-save.md
deleted file mode 100644
index 6c7d2dd..0000000
--- a/docs/src/tutorial-symbolics-save.md
+++ /dev/null
@@ -1,313 +0,0 @@
-# Cart-Pole Limit Cycle via Symbolic Lagrangian Mechanics
-
-```@meta
-Draft = true
-```
-
-This tutorial demonstrates how to combine **symbolic derivation of equations of motion** (via
-[Symbolics.jl](https://symbolics.juliasymbolics.org/)) with **direct optimal control**
-(via [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl/)) to find a
-**periodic orbit** (limit cycle) for a cart-pole system.
-
-The key idea is to let the computer do the hard mechanics: we write the Lagrangian in a few
-lines, and the Euler–Lagrange equations — including mass-matrix inversion — are derived
-automatically.
-
-## The Cart-Pole System
-
-The system consists of a cart of mass ``m_c`` sliding on a frictionless horizontal rail, with
-a rigid pendulum of mass ``m_p`` and length ``l`` attached to it. A horizontal force ``u``
-(the control input) acts on the cart.
-
-```@raw html
-<figure>
-  <img src="cartpole_sketch.svg" alt="Cart-pole diagram" width="350"/>
-  <figcaption>Fig. 1 — Cart-pole system. The angle θ is measured from the upright position.</figcaption>
-</figure>
-```
-
-The configuration vector is ``q = (x,\, \theta)^\top``, where ``x`` is the cart position and
-``\theta = 0`` corresponds to the **upright** (unstable) equilibrium of the pendulum.
-
-### Positions
-
-The Cartesian positions of the two bodies are:
-
-```math
-p_c = \begin{pmatrix} x \\ 0 \end{pmatrix}, \qquad
-p_p = \begin{pmatrix} x + l\sin\theta \\ l\cos\theta \end{pmatrix}.
-```
-
-### Lagrangian
-
-The kinetic and potential energies are:
-
-```math
-T = \tfrac{1}{2}m_c\,\|\dot{p}_c\|^2 + \tfrac{1}{2}m_p\,\|\dot{p}_p\|^2
-  = \tfrac{1}{2}(m_c+m_p)\dot{x}^2
-    + m_p l\,\dot{x}\dot{\theta}\cos\theta
-    + \tfrac{1}{2}m_p l^2\dot{\theta}^2,
-```
-
-```math
-V = m_p\,g\,l\cos\theta.
-```
-
-The Lagrangian is ``\mathcal{L} = T - V``, and the virtual work of the control force gives the
-generalised force vector ``Q = (u,\, 0)^\top``.
-
-### Euler–Lagrange Equations
-
-The equations of motion follow from:
-
-```math
-\frac{d}{dt}\frac{\partial \mathcal{L}}{\partial \dot{q}_i}
-- \frac{\partial \mathcal{L}}{\partial q_i} = Q_i, \qquad i = 1,2.
-```
-
-They can be written in the standard **manipulator form**:
-
-```math
-M(q)\,\ddot{q} = -C(q,\dot{q}) + \tau(q, u),
-```
-
-where the symmetric positive-definite **mass matrix** is:
-
-```math
-M(q) =
-\begin{pmatrix}
-  m_c + m_p & m_p l \cos\theta \\
-  m_p l \cos\theta & m_p l^2
-\end{pmatrix},
-```
-
-and the right-hand side collects Coriolis/gravity terms and the control torque. Instead of
-deriving these by hand, we rely on Symbolics.jl to compute and **analytically invert**
-``M(q)`` for us.
-
-### State-Space Form
-
-Defining the state ``X = (x,\,\theta,\,\dot{x},\,\dot{\theta})^\top``, the equations of
-motion become the first-order system:
-
-```math
-\dot{X}(t) = f\!\left(X(t),\, u(t)\right) =
-\begin{pmatrix}
-  \dot{x} \\ \dot{\theta} \\ M^{-1}(q)\bigl(-C(q,\dot{q}) + \tau(q,u)\bigr)
-\end{pmatrix}.
-```
-
-## The Optimal Control Problem
-
-We look for a **limit cycle** of the nonlinear dynamics: a trajectory that returns exactly to
-its initial condition after a fixed period ``t_f``. The cost penalises the total control
-energy:
-
-```math
-\min_{u(\cdot)}\; \int_0^{t_f} u(t)^2\,\mathrm{d}t
-```
-
-```math
-\text{subject to} \quad \dot{X}(t) = f(X(t), u(t)), \quad t \in [0, t_f],
-```
-
-```math
-X(t_f) = X(0).
-```
-
-The periodicity constraint ``X(t_f) = X(0)`` — combined with a non-trivial initial condition
-— forces the solver to find an orbit rather than the trivial rest solution.
-
-## Implementation
-
-### Setup & Imports
-
-```@example main
-# ==========================================
-# Setup & Imports
-# ==========================================
-import Pkg
-Pkg.activate("Control.jl/ModelingToolkit.jl")
-
-using OptimalControl
-using Plots
-using StaticArrays
-
-import NLPModelsIpopt
-
-using Symbolics
-```
-
-### Physical Parameters and Symbolic Variables
-
-We declare all parameters both as numerical constants (for the final function evaluation) and
-as symbolic variables (for the Lagrangian computation).
-
-```@example main
-# ==========================================
-# 1. Variables and Parameters
-# ==========================================
-# Physical Constants
-const m_c_val = 5.0
-const m_p_val = 1.0
-const l_val = 2.0
-const g_val = 9.81
-const tf_val = 2.0
-
-# Symbolic Variables
-@variables t
-D = Differential(t)
-@variables m_c m_p l g u
-@variables x(t) θ(t)
-@variables v ω dv dω # Static variables to solve for accelerations
-
-q = [x, θ]
-nothing # hide
-```
-
-### Automated Kinematics and Lagrangian
-
-We express the positions, kinetic energy, potential energy, and virtual work symbolically.
-The time derivatives ``\dot{p}_c``, ``\dot{p}_p`` are computed automatically by `D.(...)`.
-
-```@example main
-# ==========================================
-# 2. Automated Kinematics & Jacobians
-# ==========================================
-p_c = [x, 0.0]
-p_p = [x + l * sin(θ), l * cos(θ)]
-F = [u, 0.0]
-
-T = 0.5 * m_c * sum(D.(p_c) .^ 2) + 0.5 * m_p * sum(D.(p_p) .^ 2)
-V = g * (m_p * p_p[2])
-P_non_conservative = transpose(D.(p_c)) * F
-nothing # hide
-```
-
-### Euler–Lagrange Equations and Mass-Matrix Inversion
-
-Starting from ``\mathcal{L} = T - V``, Symbolics.jl computes the three terms of the
-Euler–Lagrange equations and assembles the residual vector. Substituting static aliases
-``(v,\omega,\dot{v},\dot\omega)`` for the time derivatives makes it possible to identify
-the mass matrix ``M`` as the Jacobian of the residual with respect to the accelerations
-``(\dot{v}, \dot\omega)``. The system ``M\,a = -b`` is then solved analytically.
-
-```@example main
-# ==========================================
-# 3. Euler-Lagrange & Matrix Inversion
-# ==========================================
-L = T - V
-
-A = D.(Symbolics.gradient(L, D.(q)))
-B = Symbolics.gradient(L, q)
-Q = Symbolics.gradient(P_non_conservative, D.(q))
-
-# Euler-Lagrange: d/dt( dL/dq_dot ) - dL/dq = Forces
-el_eqs = expand_derivatives.(A - B - Q)
-
-# Freeze time derivatives into static algebraic variables
-sub_rules = Dict(D(x) => v, D(θ) => ω, D(D(x)) => dv, D(D(θ)) => dω)
-res = Symbolics.substitute.(el_eqs, (sub_rules,))
-
-# Extract Mass Matrix and Bias vector
-Mass = Symbolics.jacobian(res, [dv, dω])
-bias = Symbolics.substitute.(res, (Dict(dv => 0.0, dω => 0.0),))
-
-# Analytically invert the Mass Matrix (Symbolics handles 2x2 natively)
-accel = Mass \ (-bias)
-accel = Symbolics.simplify_fractions.(accel)
-
-# The fully explicit state derivatives: X_dot = [v, ω, v_dot, ω_dot]
-dx_dt = [v, ω, accel[1], accel[2]]
-nothing # hide
-```
-
-### Code Generation
-
-`build_function` compiles the symbolic expression `dx_dt` into a native Julia function.
-The `force_SA=true` flag generates a **StaticArrays** kernel, which avoids heap allocations
-inside the ODE right-hand side — important for solver performance.
-
-```@example main
-# ==========================================
-# 4. Extract Explicit Dynamics into a Julia Function
-# ==========================================
-f_expr = build_function(dx_dt, [x, θ, v, ω], u, [m_c, m_p, l, g];
-    expression=Val{false}, force_SA=true)
-f_mtk = f_expr[1]
-
-const p_vals = [m_c_val, m_p_val, l_val, g_val]
-
-cartpole_dynamics(X, U) = f_mtk(X, U, p_vals)
-nothing # hide
-```
-
-### Optimal Control Problem Definition
-
-We now formulate the optimal control problem using the `@def` macro from OptimalControl.jl.
-The periodicity condition is encoded as ``X(t_f) - X(0) = 0``, and the dynamics plug
-directly into the `Ẋ(t) == ...` constraint.
-
-```@example main
-# ==========================================
-# 5. OptimalControl.jl dynamic optimization problem
-# ==========================================
-@def cartpole_ocp begin
-    t ∈ [0, tf_val], time
-    X = (x_ocp, θ_ocp, v_ocp, ω_ocp) ∈ R⁴, state
-    F_ctrl ∈ R, control
-
-    X(0) == [0, 0, 0, 0.1]
-    X(tf_val) - X(0) == [0, 0, 0, 0]
-
-    # Inject our explicitly inverted Symbolic dynamics!
-    Ẋ(t) == cartpole_dynamics(X(t), F_ctrl(t))
-
-    ∫(F_ctrl(t)^2) → min
-end
-```
-
-### Solving the NLP
-
-The problem is transcribed into a nonlinear program using direct collocation on a uniform
-grid of 100 intervals, then handed to **Ipopt** via NLPModelsIpopt.
-
-```@example main
-# ==========================================
-# 6. Solvers and Collocation
-# ==========================================
-initial_guess = (state=[0.0, 0.0, 0.0, 0.1], control=0.0)
-
-sol = solve(cartpole_ocp; display=true, grid_size=100, init=initial_guess)
-```
-
-### Results
-
-```@example main
-# ==========================================
-# 7. Extracting Results
-# ==========================================
-println("--- Optimal Limit Cycle Found ---")
-println("Total Control Energy (∫F² dt): ", sol.objective)
-
-tsol = time_grid(sol)
-# tsol = range(sol.model.times.initial.time, sol.model.times.final.time, 1001)
-Xsol = state(sol).(tsol)
-Fsol = control(sol).(tsol)
-
-xsol = [X[1] for X in Xsol]
-θsol = [X[2] for X in Xsol]
-vsol = [X[3] for X in Xsol]
-ωsol = [X[4] for X in Xsol]
-
-qplot = plot(tsol, [xsol θsol], label=["x" "θ"], title="Configuration")
-dqplot = plot(tsol, [vsol ωsol], label=["v" "ω"], title="Velocities")
-uplot = plot(tsol, Fsol, label="u", title="Control", linetype=:steppost)
-
-plot(qplot, dqplot, uplot, layout=3, size=(800, 600))
-```
-
-The three panels show, from top to bottom: the cart position ``x`` and pendulum angle
-``\theta``, the corresponding velocities ``\dot{x}`` and ``\dot\theta``, and the optimal
-control force ``u``. The periodicity of the state trajectory confirms that a genuine limit
-cycle has been found.
diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index 12471c4..98f6705 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -1,9 +1,5 @@
 # Cart-Pole Limit Cycle via Symbolic Lagrangian Mechanics
 
-```@meta
-Draft = false
-```
-
 This tutorial demonstrates how to combine **symbolic derivation of equations of motion** (via
 [`Symbolics.jl`](https://symbolics.juliasymbolics.org/)) with **direct optimal control**
 (via [OptimalControl.jl](https://control-toolbox.org/OptimalControl.jl/)) to find a

From 7202a60e0f0e96b077c4b2bc2aeeb3d4985cbe91 Mon Sep 17 00:00:00 2001
From: abavoil <17646791+abavoil@users.noreply.github.com>
Date: Thu, 16 Apr 2026 14:08:44 +0200
Subject: [PATCH 21/21] add comment

---
 docs/src/tutorial-symbolics.md | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/docs/src/tutorial-symbolics.md b/docs/src/tutorial-symbolics.md
index 98f6705..990041b 100644
--- a/docs/src/tutorial-symbolics.md
+++ b/docs/src/tutorial-symbolics.md
@@ -130,7 +130,6 @@ D = Differential(t)
 @variables x(t) θ(t)
 
 q = [x, θ]
-
 nothing # hide
 ```
 
@@ -156,8 +155,8 @@ nothing # hide
 Starting from ``\mathcal{L} = T - V``, `Symbolics.jl` computes the terms of the Euler–Lagrange equations. To isolate the accelerations, we substitute the symbolic time derivatives with algebraic variables. This allows us to identify the **standard manipulator form** components: the mass matrix is the Jacobian of the residual with respect to the accelerations ``\ddot{q}``, and the bias vector contains the remaining terms.
 
 ```@example main
-A = D.(Symbolics.gradient(L, D.(q)))
-B = Symbolics.gradient(L, q)
+A = D.(Symbolics.gradient(L, D.(q)))  # d/dt(∂L/∂q̇)
+B = Symbolics.gradient(L, q)          # ∂L/∂q
 Q = Symbolics.gradient(P_non_conservative, D.(q))
 
 # Euler-Lagrange residual: d/dt(∂L/∂q̇) - ∂L/∂q - Q = 0