diff --git a/docs/make.jl b/docs/make.jl
index 7f9d8d2d..1fac0fec 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -2,8 +2,8 @@ using Documenter
 using CTDirect
 using CTBase
 
-#DocMeta.setdocmeta!(CTBase, :DocTestSetup, :(using CTBase); recursive = true)
-#DocMeta.setdocmeta!(CTDirect, :DocTestSetup, :(using CTDirect); recursive = true)
+DocMeta.setdocmeta!(CTBase, :DocTestSetup, :(using CTBase); recursive = true)
+DocMeta.setdocmeta!(CTDirect, :DocTestSetup, :(using CTDirect); recursive = true)
 
 makedocs(
     warnonly = [:cross_references, :autodocs_block],
@@ -11,12 +11,9 @@ makedocs(
     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"],
-        #"API"           => "api.md",
+        "Introduction"  => "index.md",
         "Tutorial"      => "tutorial.md",
+        "API"           => "api.md",
         #"Developers"    => "dev-api.md",
     ]
 )
diff --git a/docs/src/assets/goddard.png b/docs/src/assets/goddard.png
index 3130fb2e..76e5054e 100644
Binary files a/docs/src/assets/goddard.png 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
index cd50010d..62c0abb7 100644
--- a/docs/src/tutorial.md
+++ b/docs/src/tutorial.md
@@ -1,18 +1,15 @@
-# [First example: Goddard problem](@id goddard)
+# 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.
 
-```@raw html
-<img src="./assets/goddard.png" style="display: block; margin: 0 auto 20px auto>
-```
-
-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.
-
-## Basic solve
+## Problem definition
 
 First import the CTDirect and CTBase modules
 ```@example main
@@ -21,7 +18,6 @@ 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.
-**+++ can we get here the nice table with the green checkmarks for the ocp ? it seems to appear in the terminal instead**
 ```@example main
 Cd = 310
 β = 500
@@ -37,7 +33,7 @@ function F1(x)
     r, v, m = x
     return [ 0, Tmax/m, -b*Tmax ]
 end
-@def ocp1 begin
+@def ocp begin
     tf ∈ R, variable
     t ∈ [ 0, tf ], time
     x ∈ R^3, state
@@ -57,9 +53,11 @@ end
 nothing # hide
 ```
 
-We can solve the problem directly using the default options, i.e 100 time steps and a constant initial guess set at 0.1. 
+## Basic solve
+
+We can solve the problem directly using the default options.
 ```@example main
-sol1 = solve(ocp1, print_level=5)
+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. 
@@ -67,24 +65,65 @@ Then plot the solution with the state and control variables, as well as the cost
 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
 
-+++ variantes initial guess: constant
-sol2
-, fonctionnel
-sol3
-, mixed
-sol4
+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")
+``` 
 
-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).
+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
-sol4 = solve(ocp1, grid_size=200, print_level=5, init=OptimalControlInit(sol1))
+init4 = OptimalControlInit(sol1)
+sol4 = solve(ocp, grid_size=200, print_level=5, init=init4)
 nothing # hide
 ```
 ```@example main
 plot(sol4)
 ```
 
-## Working explicitely on the discretized problem
-+++ variantes appel: transcription puis solve docp, avec init changee dans le docp (ex: mixed) et passee au solve (ex: warmstart)
+## 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/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 e50c46d4..790ece4a 100644
--- a/test/test_initial_guess.jl
+++ b/test/test_initial_guess.jl
@@ -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)