diff --git a/.github/workflows/Documentation.yml.disabled b/.github/workflows/Documentation.yml
similarity index 100%
rename from .github/workflows/Documentation.yml.disabled
rename to .github/workflows/Documentation.yml
diff --git a/docs/Project.toml b/docs/Project.toml
index 82a5f094..dec9ddc4 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -1,3 +1,4 @@
 [deps]
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80"
diff --git a/docs/make.jl b/docs/make.jl
index 4b04a7de..1fac0fec 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -6,13 +6,13 @@ DocMeta.setdocmeta!(CTBase, :DocTestSetup, :(using CTBase); recursive = true)
 DocMeta.setdocmeta!(CTDirect, :DocTestSetup, :(using CTDirect); recursive = true)
 
 makedocs(
+    warnonly = [:cross_references, :autodocs_block],
     sitename = "CTDirect.jl",
-    format = Documenter.HTML(prettyurls = false),
+    format = Documenter.HTML(prettyurls = false,
+            size_threshold_ignore = ["api-ctbase.md"]),
     pages = [
         "Introduction"  => "index.md",
-        #"Methods and Options" => ["rk.md",
-        #                          "optimization.md",
-        #                          "solve-options.md"],
+        "Tutorial"      => "tutorial.md",
         "API"           => "api.md",
         #"Developers"    => "dev-api.md",
     ]
diff --git a/docs/src/api-ctbase.md b/docs/src/api-ctbase.md
index 1c543b60..6d8698f1 100644
--- a/docs/src/api-ctbase.md
+++ b/docs/src/api-ctbase.md
@@ -9,7 +9,6 @@ For more details about `CTBase.jl` package, see the [documentation](https://cont
 Pages   = ["api-ctbase.md"]
 Modules = [CTBase]
 Order = [:module, :constant, :type, :function, :macro]
-Private = false
 ```
 
 ## Documentation
diff --git a/docs/src/api.md b/docs/src/api.md
index fc55b2b9..d120f0ee 100644
--- a/docs/src/api.md
+++ b/docs/src/api.md
@@ -6,7 +6,6 @@
 Pages   = ["api.md"]
 Modules = [CTDirect]
 Order = [:module, :constant, :type, :function, :macro]
-Private = false
 ```
 
 ## Available methods
diff --git a/docs/src/assets/goddard.png b/docs/src/assets/goddard.png
new file mode 100644
index 00000000..76e5054e
Binary files /dev/null and b/docs/src/assets/goddard.png differ
diff --git a/docs/src/optimization.md b/docs/src/optimization.md
deleted file mode 100644
index 147f49d1..00000000
--- a/docs/src/optimization.md
+++ /dev/null
@@ -1,2 +0,0 @@
-# Optimization
-    For the moment we have implemented only Ipopt for the optimization soltware
\ No newline at end of file
diff --git a/docs/src/rk.md b/docs/src/rk.md
deleted file mode 100644
index 1a79836b..00000000
--- a/docs/src/rk.md
+++ /dev/null
@@ -1,2 +0,0 @@
-# Runge Kutta method
-For the moment we have implemented only the trapeze method for the discretization
\ No newline at end of file
diff --git a/docs/src/solve-options.md b/docs/src/solve-options.md
deleted file mode 100644
index 6803cb34..00000000
--- a/docs/src/solve-options.md
+++ /dev/null
@@ -1,4 +0,0 @@
-# solve options
-
-
-See [`solve`](@ref) method.
\ No newline at end of file
diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md
new file mode 100644
index 00000000..62c0abb7
--- /dev/null
+++ b/docs/src/tutorial.md
@@ -0,0 +1,129 @@
+# First example: Goddard problem
+
+![R.H. Goddard](./assets/goddard.png)
+
+We start with the well-known so-called [Goddard](http://en.wikipedia.org/wiki/Robert_H._Goddard) problem, which models the ascent of a rocket through the atmosphere (see for instance [^1],[^2]).
+We restrict here ourselves to  vertical (monodimensional) trajectories, and the state variables are the altitude, speed and mass of the rocket during the flight, for a total dimension of 3. The rocket is subject to gravity, thrust and drag forces. The final time is free, and the objective is here to reach a maximal altitude with a given fuel consumption, i.e a fixed final mass. We also impose a speed limit. All units are renormalized.
+
+[^1]: R.H. Goddard. A Method of Reaching Extreme Altitudes, volume 71(2) of Smithsonian Miscellaneous Collections. Smithsonian institution, City of Washington, 1919.
+
+[^2]: H. Seywald and E.M. Cliff. Goddard problem in presence of a dynamic pressure limit. Journal of Guidance, Control, and Dynamics, 16(4):776–781, 1993.
+
+## Problem definition
+
+First import the CTDirect and CTBase modules
+```@example main
+using CTDirect
+using CTBase
+```
+
+Then define the OCP for the Goddard problem. Note that the free final time is modeled as an *optimization variable*, and has both a lower bound to prevent the optimization getting stuck at tf=0. In this particular case an upper bound is not needed for the final time since the final mass is prescribed.
+```@example main
+Cd = 310
+β = 500
+Tmax = 3.5
+b = 2
+vmax = 0.1
+function F0(x)
+    r, v, m = x
+    D = Cd * v^2 * exp(-β*(r - 1))
+    return [ v, -D/m - 1/r^2, 0 ]
+end
+function F1(x)
+    r, v, m = x
+    return [ 0, Tmax/m, -b*Tmax ]
+end
+@def ocp begin
+    tf ∈ R, variable
+    t ∈ [ 0, tf ], time
+    x ∈ R^3, state
+    u ∈ R, control
+    0.1 ≤ tf ≤ Inf
+    r = x[1]
+    v = x[2]
+    m = x[3]
+    x(0) == [1, 0, 1]
+    m(tf) == 0.6
+    1 ≤ r(t) ≤ 1.1
+    0 ≤ v(t) ≤ vmax
+    0 ≤ u(t) ≤ 1
+    ẋ(t) == F0(x(t)) + u(t)*F1(x(t))
+    r(tf) → max
+end
+nothing # hide
+```
+
+## Basic solve
+
+We can solve the problem directly using the default options.
+```@example main
+sol1 = solve(ocp, print_level=5)
+nothing # hide
+```
+Then plot the solution with the state and control variables, as well as the costate recovered from the lagrange multipliers of the discretized problem. 
+```@example main
+plot(sol1)
+```
+
+The most common option for **solve** is the number of time steps for the discretized problem (default 100), that can be set with the argument *grid_size*. A larger grid size will increase the computational cost, while a smaller value may lead to a very coarse solution.
+
+## Initial guess options
+
+The function **solve** uses a default constant initialisation of 0.1 for all variables. More advanced options include constant and/or functional initialisation for each individual state or control component, as well as reusing an existing solution, also known as *warm start*[^3].
+
+[^3]: Currently only the primal variables are reused for the warm start, not the lagrange multipliers. It should be noted that previous experiments with the Bocop software seemed to indicate that initializing also the multipliers gave little benefit.
+
+Let us start with the simplest case, constant initialisation.
+```@example main
+x_const = [1.05, 0.2, 0.8]
+u_const = 0.5
+v_const = 0.15
+init1 = OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const)
+sol2 = solve(ocp, print_level=0, init=init1)
+println("Objective ", sol2.objective, " after ", sol2.iterations, " iterations")
+``` 
+
+Now we illustrate the functional initialisation, with some random functions. Note that we only consider the state and control variables, since the optimization variables are scalar and therefore a functional initialisation is not relevant. In the example notice that the call to **OptimalControlInit** does not provide an argument for the optimization variables, therefore the default initial guess will be used.  
+```@example main
+x_func = t->[1+t^2, sqrt(t), 1-t]
+u_func = t->(cos(t)+1)*0.5
+init2 = OptimalControlInit(x_init=x_func, u_init=u_func)
+sol3 = solve(ocp, print_level=0, init=init2)
+println("Objective ", sol3.objective, " after ", sol3.iterations, " iterations")
+``` 
+More generally, the default, constant and functional initialisations can be mixed, as shown in the example below that uses a functional initial guess for the state, a constant initial guess for the control, and the default initial guess for the optimization variables. 
+```@example main
+init3 = OptimalControlInit(x_init=x_func, u_init=u_const)
+sol4 = solve(ocp, print_level=0, init=init3)
+println("Objective ", sol4.objective, " after ", sol4.iterations, " iterations")
+``` 
+
+Finally, we can also use a so-called *warmstart* strategy and use an existing solution as initial guess (note that the OCP solution returned by the **solve** call is functional, thus it is not necessary to use the same time grid). Notice that the objective and constraint violation values start much closer to the solution than with the previous initialisations.
+```@example main
+init4 = OptimalControlInit(sol1)
+sol4 = solve(ocp, grid_size=200, print_level=5, init=init4)
+nothing # hide
+```
+```@example main
+plot(sol4)
+```
+
+## The discretized problem
+
+Instead of calling **solve** directly on the OCP problem, you can first obtain the discretized problem (DOCP) by calling **directTranscription**, then call **solve** on the DOCP.
+```@example main
+docp = directTranscription(ocp, grid_size=100)
+sol5 = solve(docp, print_level=5)
+nothing # hide
+```
+The initial guess can be passed to **solve** same as before.
+```@example main
+sol6 = solve(docp, print_level=0, init=OptimalControlInit(sol1))
+println("Objective ", sol6.objective, " after ", sol6.iterations, " iterations")
+```
+Another possibility is to set the initial guess associated to the DOCP, using the function **setDOCPInit**.
+```@example main
+setDOCPInit(docp, OptimalControlInit(sol1))
+sol7 = solve(docp, print_level=5)
+nothing # hide
+```
diff --git a/test/suite/abstract_ocp.jl b/test/suite/abstract_ocp.jl
index 83bffb07..b009efc3 100644
--- a/test/suite/abstract_ocp.jl
+++ b/test/suite/abstract_ocp.jl
@@ -45,4 +45,42 @@ end
 @testset verbose = true showtiming = true ":double_integrator :min_tf :abstract :constr" begin
     sol2 = solve(ocp2, grid_size=100, print_level=0, tol=1e-12)
     @test sol2.objective ≈ 5.46 rtol=1e-2
+end
+
+Cd = 310
+β = 500
+Tmax = 3.5
+b = 2
+vmax = 0.1
+function F0(x)
+    r, v, m = x
+    D = Cd * v^2 * exp(-β*(r - 1))
+    return [ v, -D/m - 1/r^2, 0 ]
+end
+function F1(x)
+    r, v, m = x
+    return [ 0, Tmax/m, -b*Tmax ]
+end
+# NB. the ≤ is not the same as <= (parse error for <=)
+@def ocp3 begin
+    tf ∈ R, variable
+    t ∈ [ 0, tf ], time
+    x ∈ R^3, state
+    u ∈ R, control
+    0.1 ≤ tf ≤ Inf
+    r = x[1]
+    v = x[2]
+    m = x[3]
+    x(0) == [1, 0, 1]
+    m(tf) == 0.6
+    1 ≤ r(t) ≤ 1.1
+    0 ≤ v(t) ≤ vmax
+    0 ≤ u(t) ≤ 1
+    ẋ(t) == F0(x(t)) + u(t)*F1(x(t))
+    r(tf) → max
+end
+
+@testset verbose = true showtiming = true ":double_integrator :min_tf :abstract :constr" begin
+    sol3 = solve(ocp3, grid_size=100, print_level=0, tol=1e-12)    
+    @test sol3.objective ≈ 1.0125 rtol=1e-2
 end
\ No newline at end of file
diff --git a/test/suite/initial_guess.jl b/test/suite/initial_guess.jl
index 1a117bdd..621f6d57 100644
--- a/test/suite/initial_guess.jl
+++ b/test/suite/initial_guess.jl
@@ -58,7 +58,7 @@ end
 # constant initial guess
 x_const = [1.05, 0.2, 0.8]
 u_const = 0.5
-v_init = 0.15
+v_const = 0.15
 
 @testset verbose = true showtiming = true ":constant_init_x" begin
     sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const))
@@ -69,7 +69,7 @@ end
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":constant_init_v" begin
-    sol = solve(ocp, print_level=0, init=OptimalControlInit(v_init=v_init))
+    sol = solve(ocp, print_level=0, init=OptimalControlInit(v_init=v_const))
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":constant_init_xu" begin
@@ -77,15 +77,15 @@ end
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":constant_init_xv" begin
-    sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, v_init=v_init))
+    sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, v_init=v_const))
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":constant_init_uv" begin
-    sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const, v_init=v_init))
+    sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const, v_init=v_const))
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":constant_init_xuv" begin
-    sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init))
+    sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const))
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 
@@ -119,7 +119,7 @@ end
 # set initial guess in DOCP
 docp = directTranscription(ocp)
 @testset verbose = true showtiming = true ":DOCPInit_constant" begin
-    setDOCPInit(docp, OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init))
+    setDOCPInit(docp, OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const))
     sol = solve(docp, print_level=0)
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
@@ -137,7 +137,7 @@ end
 # pass initial guess to solve
 setDOCPInit(docp, OptimalControlInit()) # reset init in docp
 @testset verbose = true showtiming = true ":solve_constant_init" begin
-    sol = solve(docp, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init), print_level=0)
+    sol = solve(docp, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const), print_level=0)
     @test sol.objective ≈ 1.0125 rtol=1e-2
 end
 @testset verbose = true showtiming = true ":solve_mixed_init" begin
diff --git a/test/test_initial_guess.jl b/test/test_initial_guess.jl
index f1f40648..790ece4a 100644
--- a/test/test_initial_guess.jl
+++ b/test/test_initial_guess.jl
@@ -52,7 +52,7 @@ end
 dynamics!(ocp, (x, u, v) -> F0(x) + u*F1(x) )
 
 # reference solution
-sol0 = solve(ocp, print_level=0, tol=1e-8)
+sol0 = solve(ocp, print_level=0)
 #= check ok
 tf0 = sol0.variable
 println("tf ", tf0, " obj ", sol0.objective)
@@ -82,7 +82,7 @@ sol = solve(ocp, print_level=0, max_iter=maxiter)
 # constant initial guess
 x_const = [1.05, 0.2, 0.8]
 u_const = 0.5
-v_init = 0.15
+v_const = 0.15
 
 # Constant initial guess (vector for x; default for u,v)
 sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const), max_iter=maxiter)
@@ -93,7 +93,7 @@ sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const), max_ite
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for u; default for x,v):", sol.objective, sol.iterations)
 
 # Constant initial guess (vector for v; default for x,u)
-sol = solve(ocp, print_level=0, init=OptimalControlInit(v_init=v_init), max_iter=maxiter)
+sol = solve(ocp, print_level=0, init=OptimalControlInit(v_init=v_const), max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for v; default for x,u):", sol.objective, sol.iterations)
 
 # Constant initial guess (vector for x,u; default for v)
@@ -101,15 +101,15 @@ sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for x,u; default for v):", sol.objective, sol.iterations)
 
 # Constant initial guess (vector for x,v; default for u)
-sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, v_init=v_init), max_iter=maxiter)
+sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, v_init=v_const), max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for x,v; default for u):", sol.objective, sol.iterations)
 
 # Constant initial guess (vector for u,v; default for x)
-sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const, v_init=v_init), max_iter=maxiter)
+sol = solve(ocp, print_level=0, init=OptimalControlInit(u_init=u_const, v_init=v_const), max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for u,v; default for x):", sol.objective, sol.iterations)
 
 # Constant initial guess (vector for x,u,v)
-sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init), max_iter=maxiter)
+sol = solve(ocp, print_level=0, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const), max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess (vector for x,u,v):", sol.objective, sol.iterations)
 
 # functional initial guess
@@ -141,7 +141,7 @@ sol = solve(ocp, print_level=0, init=init_function_u = OptimalControlInit(sol0),
 println("\nSetting the initial guess at the DOCP level")
 docp = directTranscription(ocp)
 # constant vector init
-setDOCPInit(docp, OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init))
+setDOCPInit(docp, OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const))
 sol = solve(docp, print_level=0, max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "Constant initial guess set in DOCP", sol.objective, sol.iterations)
 # mixed init
@@ -158,7 +158,7 @@ sol = solve(docp, print_level=0, max_iter=maxiter)
 println("\nPassing the initial guess to solve call")
 setDOCPInit(docp, OptimalControlInit()) # reset init in docp
 # constant vector init
-sol = solve(docp, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_init), print_level=0, max_iter=maxiter)
+sol = solve(docp, init=OptimalControlInit(x_init=x_const, u_init=u_const, v_init=v_const), print_level=0, max_iter=maxiter)
 @printf("%-56s %.3f at %d iterations\n", "constant initial guess passed to solve", sol.objective, sol.iterations)
 # mixed init
 sol = solve(docp, init=OptimalControlInit(x_init=x_func, u_init=u_const), print_level=0, max_iter=maxiter)