diff --git a/test/test_model.jl b/test/test_model.jl
index b0cd051..29a3444 100644
--- a/test/test_model.jl
+++ b/test/test_model.jl
@@ -781,8 +781,8 @@ function test_model() # 30 55 185
         )
         constraint!(ocp, :initial, lb = 0, ub = 0, label = :c0)
         constraint!(ocp, :final, lb = 1, ub = 1, label = :cf)
-        @test __constraint(ocp, :c0)(12, ∅) == 12
-        @test __constraint(ocp, :cf)(∅, 12) == 12
+        @test constraint(ocp, :c0)(12, ∅) == 12
+        @test constraint(ocp, :cf)(∅, 12) == 12
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
@@ -790,8 +790,8 @@ function test_model() # 30 55 185
         control!(ocp, 1)
         constraint!(ocp, :initial, lb = [0, 1], ub = [0, 1], label = :c0)
         constraint!(ocp, :final, lb = [1, 2], ub = [1, 2], label = :cf)
-        @test __constraint(ocp, :c0)([12, 13], ∅) == [12, 13]
-        @test __constraint(ocp, :cf)(∅, [12, 13]) == [12, 13]
+        @test constraint(ocp, :c0)([12, 13], ∅) == [12, 13]
+        @test constraint(ocp, :cf)(∅, [12, 13]) == [12, 13]
 
         # constraint already exists
         ocp = Model()
@@ -812,13 +812,13 @@ function test_model() # 30 55 185
         xf = 1
         constraint!(ocp, :initial, rg = 1, lb = x0, ub = x0, label = :c0)
         constraint!(ocp, :final, rg = 1, lb = xf, ub = xf, label = :cf)
-        @test __constraint(ocp, :c0)(x, ∅) == x
-        @test __constraint(ocp, :cf)(∅, x) == x
+        @test constraint(ocp, :c0)(x, ∅) == x
+        @test constraint(ocp, :cf)(∅, x) == x
 
         constraint!(ocp, :initial, rg = 1, lb = x0, ub = x0, label = :c00)
         constraint!(ocp, :final, rg = 1, lb = xf, ub = xf, label = :cff)
-        @test __constraint(ocp, :c00)(x, ∅) == x
-        @test __constraint(ocp, :cff)(∅, x) == x
+        @test constraint(ocp, :c00)(x, ∅) == x
+        @test constraint(ocp, :cff)(∅, x) == x
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
@@ -853,8 +853,8 @@ function test_model() # 30 55 185
         xf = [1, 2]
         constraint!(ocp, :initial, rg = 1:2, lb = x0, ub = x0, label = :c0)
         constraint!(ocp, :final, rg = 1:2, lb = xf, ub = xf, label = :cf)
-        @test __constraint(ocp, :c0)(x, ∅) == x[1:2]
-        @test __constraint(ocp, :cf)(∅, x) == x[1:2]
+        @test constraint(ocp, :c0)(x, ∅) == x[1:2]
+        @test constraint(ocp, :cf)(∅, x) == x[1:2]
 
         # constraint already exists
         ocp = Model()
@@ -874,10 +874,10 @@ function test_model() # 30 55 185
         constraint!(ocp, :final, lb = 1, ub = 2, label = :cf)
         constraint!(ocp, :control, lb = 0, ub = 1, label = :cu)
         constraint!(ocp, :state, lb = 0, ub = 1, label = :cs)
-        @test __constraint(ocp, :c0)(12, ∅) == 12
-        @test __constraint(ocp, :cf)(∅, 12) == 12
-        @test __constraint(ocp, :cu)(12) == 12
-        @test __constraint(ocp, :cs)(12) == 12
+        @test constraint(ocp, :c0)(12, ∅) == 12
+        @test constraint(ocp, :cf)(∅, 12) == 12
+        @test constraint(ocp, :cu)(12) == 12
+        @test constraint(ocp, :cs)(12) == 12
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
@@ -887,10 +887,10 @@ function test_model() # 30 55 185
         constraint!(ocp, :final, lb = [1, 2], ub = [2, 3], label = :cf)
         constraint!(ocp, :control, lb = [0, 1], ub = [1, 2], label = :cu)
         constraint!(ocp, :state, lb = [0, 1], ub = [1, 2], label = :cs)
-        @test __constraint(ocp, :c0)([12, 13], ∅) == [12, 13]
-        @test __constraint(ocp, :cf)(∅, [12, 13]) == [12, 13]
-        @test __constraint(ocp, :cu)([12, 13]) == [12, 13]
-        @test __constraint(ocp, :cs)([12, 13]) == [12, 13]
+        @test constraint(ocp, :c0)([12, 13], ∅) == [12, 13]
+        @test constraint(ocp, :cf)(∅, [12, 13]) == [12, 13]
+        @test constraint(ocp, :cu)([12, 13]) == [12, 13]
+        @test constraint(ocp, :cs)([12, 13]) == [12, 13]
 
         # constraint already exists
         ocp = Model()
@@ -910,10 +910,10 @@ function test_model() # 30 55 185
         constraint!(ocp, :final, rg = 1, lb = 1, ub = 2, label = :cf)
         constraint!(ocp, :control, rg = 1, lb = 0, ub = 1, label = :cu)
         constraint!(ocp, :state, rg = 1, lb = 0, ub = 1, label = :cs)
-        @test __constraint(ocp, :c0)(12, ∅) == 12
-        @test __constraint(ocp, :cf)(∅, 12) == 12
-        @test __constraint(ocp, :cu)(12) == 12
-        @test __constraint(ocp, :cs)(12) == 12
+        @test constraint(ocp, :c0)(12, ∅) == 12
+        @test constraint(ocp, :cf)(∅, 12) == 12
+        @test constraint(ocp, :cu)(12) == 12
+        @test constraint(ocp, :cs)(12) == 12
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
@@ -960,10 +960,10 @@ function test_model() # 30 55 185
         constraint!(ocp, :final, rg = 1:2, lb = [1, 2], ub = [2, 3], label = :cf)
         constraint!(ocp, :control, rg = 1:2, lb = [0, 1], ub = [1, 2], label = :cu)
         constraint!(ocp, :state, rg = 1:2, lb = [0, 1], ub = [1, 2], label = :cs)
-        @test __constraint(ocp, :c0)([12, 13], ∅) == [12, 13]
-        @test __constraint(ocp, :cf)(∅, [12, 13]) == [12, 13]
-        @test __constraint(ocp, :cu)([12, 13]) == [12, 13]
-        @test __constraint(ocp, :cs)([12, 13]) == [12, 13]
+        @test constraint(ocp, :c0)([12, 13], ∅) == [12, 13]
+        @test constraint(ocp, :cf)(∅, [12, 13]) == [12, 13]
+        @test constraint(ocp, :cu)([12, 13]) == [12, 13]
+        @test constraint(ocp, :cs)([12, 13]) == [12, 13]
 
         # constraint already exists
         ocp = Model()
@@ -983,10 +983,10 @@ function test_model() # 30 55 185
         constraint!(ocp, :control, f = u -> u, lb = 0, ub = 0, label = :cu)
         constraint!(ocp, :state, f = x -> x, lb = 0, ub = 0, label = :cs)
         constraint!(ocp, :mixed, f = (x, u) -> x + u, lb = 1, ub = 1, label = :cm)
-        @test __constraint(ocp, :cb)(12, 13) == 12 + 13
-        @test __constraint(ocp, :cu)(12) == 12
-        @test __constraint(ocp, :cs)(12) == 12
-        @test __constraint(ocp, :cm)(12, 13) == 12 + 13
+        @test constraint(ocp, :cb)(12, 13) == 12 + 13
+        @test constraint(ocp, :cu)(12) == 12
+        @test constraint(ocp, :cs)(12) == 12
+        @test constraint(ocp, :cm)(12, 13) == 12 + 13
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
@@ -996,10 +996,10 @@ function test_model() # 30 55 185
         constraint!(ocp, :control, f = u -> u[1], lb = 0, ub = 0, label = :cu)
         constraint!(ocp, :state, f = x -> x[1], lb = 0, ub = 0, label = :cs)
         constraint!(ocp, :mixed, f = (x, u) -> x[1] + u[1], lb = 1, ub = 1, label = :cm)
-        @test __constraint(ocp, :cb)([13, 14], [16, 17]) == 13 + 16
-        @test __constraint(ocp, :cu)([12, 13]) == 12
-        @test __constraint(ocp, :cs)([12, 13]) == 12
-        @test __constraint(ocp, :cm)([12, 13], [14, 15]) == 12 + 14
+        @test constraint(ocp, :cb)([13, 14], [16, 17]) == 13 + 16
+        @test constraint(ocp, :cu)([12, 13]) == 12
+        @test constraint(ocp, :cs)([12, 13]) == 12
+        @test constraint(ocp, :cm)([12, 13], [14, 15]) == 12 + 14
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
@@ -1023,10 +1023,10 @@ function test_model() # 30 55 185
             ub = [0, 0],
             label = :cm,
         )
-        @test __constraint(ocp, :cb)([13, 14, 15], [17, 18, 19]) == [13 + 17, 14 + 18]
-        @test __constraint(ocp, :cu)([12, 13, 14]) == [12, 13]
-        @test __constraint(ocp, :cs)([12, 13, 14]) == [12, 13]
-        @test __constraint(ocp, :cm)([12, 13, 14], [15, 16, 17]) == [12 + 15, 13 + 16]
+        @test constraint(ocp, :cb)([13, 14, 15], [17, 18, 19]) == [13 + 17, 14 + 18]
+        @test constraint(ocp, :cu)([12, 13, 14]) == [12, 13]
+        @test constraint(ocp, :cs)([12, 13, 14]) == [12, 13]
+        @test constraint(ocp, :cm)([12, 13, 14], [15, 16, 17]) == [12 + 15, 13 + 16]
 
         # constraint already exists
         ocp = Model()
@@ -1053,10 +1053,10 @@ function test_model() # 30 55 185
         constraint!(ocp, :control, f = u -> u, lb = 0, ub = 1, label = :cu)
         constraint!(ocp, :state, f = x -> x, lb = 0, ub = 1, label = :cs)
         constraint!(ocp, :mixed, f = (x, u) -> x + u, lb = 1, ub = 1, label = :cm)
-        @test __constraint(ocp, :cb)(12, 13) == 12 + 13
-        @test __constraint(ocp, :cu)(12) == 12
-        @test __constraint(ocp, :cs)(12) == 12
-        @test __constraint(ocp, :cm)(12, 13) == 12 + 13
+        @test constraint(ocp, :cb)(12, 13) == 12 + 13
+        @test constraint(ocp, :cu)(12) == 12
+        @test constraint(ocp, :cs)(12) == 12
+        @test constraint(ocp, :cm)(12, 13) == 12 + 13
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
@@ -1066,10 +1066,10 @@ function test_model() # 30 55 185
         constraint!(ocp, :control, f = u -> u[1], lb = 0, ub = 1, label = :cu)
         constraint!(ocp, :state, f = x -> x[1], lb = 0, ub = 1, label = :cs)
         constraint!(ocp, :mixed, f = (x, u) -> x[1] + u[1], lb = 1, ub = 1, label = :cm)
-        @test __constraint(ocp, :cb)([13, 14], [16, 17]) == 13 + 16
-        @test __constraint(ocp, :cu)([12, 13]) == 12
-        @test __constraint(ocp, :cs)([12, 13]) == 12
-        @test __constraint(ocp, :cm)([12, 13], [14, 15]) == 12 + 14
+        @test constraint(ocp, :cb)([13, 14], [16, 17]) == 13 + 16
+        @test constraint(ocp, :cu)([12, 13]) == 12
+        @test constraint(ocp, :cs)([12, 13]) == 12
+        @test constraint(ocp, :cm)([12, 13], [14, 15]) == 12 + 14
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
@@ -1093,10 +1093,10 @@ function test_model() # 30 55 185
             ub = [1, 1],
             label = :cm,
         )
-        @test __constraint(ocp, :cb)([13, 14, 15], [17, 18, 19]) == [13 + 17, 14 + 18]
-        @test __constraint(ocp, :cu)([12, 13, 14]) == [12, 13]
-        @test __constraint(ocp, :cs)([12, 13, 14]) == [12, 13]
-        @test __constraint(ocp, :cm)([12, 13, 14], [15, 16, 17]) == [12 + 15, 13 + 16]
+        @test constraint(ocp, :cb)([13, 14, 15], [17, 18, 19]) == [13 + 17, 14 + 18]
+        @test constraint(ocp, :cu)([12, 13, 14]) == [12, 13]
+        @test constraint(ocp, :cs)([12, 13, 14]) == [12, 13]
+        @test constraint(ocp, :cm)([12, 13, 14], [15, 16, 17]) == [12 + 15, 13 + 16]
     end
 
     @testset "constraint! 7" begin
@@ -1122,11 +1122,11 @@ function test_model() # 30 55 185
             ub = [1, 0, 1, 0],
             label = :eq5,
         )
-        @test __constraint(ocp, :eq1)(v) == v
-        @test __constraint(ocp, :eq2)(v) == v[1]
-        @test __constraint(ocp, :eq3)(v) == v[1:2]
-        @test __constraint(ocp, :eq4)(v) == v[1:2:4]
-        @test __constraint(ocp, :eq5)(v) == v .^ 2
+        @test constraint(ocp, :eq1)(v) == v
+        @test constraint(ocp, :eq2)(v) == v[1]
+        @test constraint(ocp, :eq3)(v) == v[1:2]
+        @test constraint(ocp, :eq4)(v) == v[1:2:4]
+        @test constraint(ocp, :eq5)(v) == v .^ 2
     end
 
     @testset "constraint! 8" begin
@@ -1156,15 +1156,15 @@ function test_model() # 30 55 185
             ub = [0, 0, 0, 0],
             label = :eq5,
         )
-        @test __constraint(ocp, :cb)(x0, xf, v) == x0 + xf + v[1]
-        @test __constraint(ocp, :cu)(u, v) == u + v[1]
-        @test __constraint(ocp, :cs)(x, v) == x + v[1]
-        @test __constraint(ocp, :cm)(x, u, v) == x + u + v[1]
-        @test __constraint(ocp, :eq1)(v) == v
-        @test __constraint(ocp, :eq2)(v) == v[1]
-        @test __constraint(ocp, :eq3)(v) == v[1:2]
-        @test __constraint(ocp, :eq4)(v) == v[1:2:4]
-        @test __constraint(ocp, :eq5)(v) == v .^ 2
+        @test constraint(ocp, :cb)(x0, xf, v) == x0 + xf + v[1]
+        @test constraint(ocp, :cu)(u, v) == u + v[1]
+        @test constraint(ocp, :cs)(x, v) == x + v[1]
+        @test constraint(ocp, :cm)(x, u, v) == x + u + v[1]
+        @test constraint(ocp, :eq1)(v) == v
+        @test constraint(ocp, :eq2)(v) == v[1]
+        @test constraint(ocp, :eq3)(v) == v[1:2]
+        @test constraint(ocp, :eq4)(v) == v[1:2:4]
+        @test constraint(ocp, :eq5)(v) == v .^ 2
     end
 
     @testset "constraint! 9" begin
@@ -1173,28 +1173,28 @@ function test_model() # 30 55 185
         state!(ocp, 1)
         control!(ocp, 1)
         dynamics!(ocp, (x, u) -> x + u)
-        @test __dynamics(ocp)(1, 2) == 3
+        @test dynamics(ocp)(1, 2) == 3
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
         state!(ocp, 2)
         control!(ocp, 2)
         dynamics!(ocp, (x, u) -> x[1] + u[1])
-        @test __dynamics(ocp)([1, 2], [3, 4]) == 4
+        @test dynamics(ocp)([1, 2], [3, 4]) == 4
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
         state!(ocp, 2)
         control!(ocp, 2)
         dynamics!(ocp, (x, u) -> [x[1] + u[1], x[2] + u[2]])
-        @test __dynamics(ocp)([1, 2], [3, 4]) == [4, 6]
+        @test dynamics(ocp)([1, 2], [3, 4]) == [4, 6]
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
         state!(ocp, 2)
         control!(ocp, 1)
         dynamics!(ocp, (x, u) -> [x[1] + u, x[2] + u])
-        @test __dynamics(ocp)([1, 2], 3) == [4, 5]
+        @test dynamics(ocp)([1, 2], 3) == [4, 5]
     end
 
     @testset "constraint! 10" begin
@@ -1204,8 +1204,8 @@ function test_model() # 30 55 185
         control!(ocp, 1)
         constraint!(ocp, :initial, lb = 0, ub = 1, label = :c0)
         constraint!(ocp, :final, lb = 1, ub = 2, label = :cf)
-        @test __constraint(ocp, :c0)(12, ∅) == 12
-        @test __constraint(ocp, :cf)(∅, 12) == 12
+        @test constraint(ocp, :c0)(12, ∅) == 12
+        @test constraint(ocp, :cf)(∅, 12) == 12
     end
 
     @testset "constraint! 11" begin
@@ -1888,60 +1888,60 @@ function test_model() # 30 55 185
         state!(ocp, 1)
         control!(ocp, 1)
         objective!(ocp, :lagrange, (x, u) -> 0.5u^2)
-        @test __lagrange(ocp)(1, 2) == 2
+        @test lagrange(ocp)(1, 2) == 2
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
         state!(ocp, 2)
         control!(ocp, 2)
         objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2)
-        @test __lagrange(ocp)([1, 2], [3, 4]) == 4.5
+        @test lagrange(ocp)([1, 2], [3, 4]) == 4.5
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
         state!(ocp, 1)
         control!(ocp, 1)
         objective!(ocp, :mayer, (x0, xf) -> 0.5x0^2)
-        @test __mayer(ocp)(2, 3) == 2
+        @test mayer(ocp)(2, 3) == 2
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
         state!(ocp, 2)
         control!(ocp, 2)
         objective!(ocp, :mayer, (x0, xf) -> 0.5x0[1]^2)
-        @test __mayer(ocp)([2, 3], [5, 6]) == 2
+        @test mayer(ocp)([2, 3], [5, 6]) == 2
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
         state!(ocp, 1)
         control!(ocp, 1)
         objective!(ocp, :bolza, (x0, xf) -> 0.5x0^2, (x, u) -> 0.5u^2)
-        @test __mayer(ocp)(2, 3) == 2
-        @test __lagrange(ocp)(1, 2) == 2
+        @test mayer(ocp)(2, 3) == 2
+        @test lagrange(ocp)(1, 2) == 2
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
         state!(ocp, 2)
         control!(ocp, 2)
         objective!(ocp, :bolza, (x0, xf) -> 0.5x0[1]^2, (x, u) -> 0.5u[1]^2)
-        @test __mayer(ocp)([2, 3], [5, 6]) == 2
-        @test __lagrange(ocp)([1, 2], [3, 4]) == 4.5
+        @test mayer(ocp)([2, 3], [5, 6]) == 2
+        @test lagrange(ocp)([1, 2], [3, 4]) == 4.5
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
         state!(ocp, 2)
         control!(ocp, 2)
         objective!(ocp, :lagrange, (x, u) -> 0.5u[1]^2)
-        @test __lagrange(ocp)([1, 2], [3, 4]) == 4.5
-        @test isnothing(__mayer(ocp))
+        @test lagrange(ocp)([1, 2], [3, 4]) == 4.5
+        @test isnothing(mayer(ocp))
 
         ocp = Model()
         time!(ocp; t0 = 0, tf = 1)
         state!(ocp, 2)
         control!(ocp, 2)
         objective!(ocp, :mayer, (x0, xf) -> 0.5x0[1]^2)
-        @test __mayer(ocp)([2, 3], [5, 6]) == 2
-        @test isnothing(__lagrange(ocp))
+        @test mayer(ocp)([2, 3], [5, 6]) == 2
+        @test isnothing(lagrange(ocp))
     end
 
     @testset "redeclarations" begin
diff --git a/test/test_plot_joseph.jl b/test/test_plot_joseph.jl
deleted file mode 100644
index 0f40155..0000000
--- a/test/test_plot_joseph.jl
+++ /dev/null
@@ -1,54 +0,0 @@
-# using OptimalControl
-
-# @def ocp begin
-
-#     t ∈ [ 0, 1 ], time
-#     x ∈ R, state
-#     u ∈ R, control
-#     x(0) == 0
-#     -1 ≤ u(t) ≤ 1
-#     ẋ(t) == -x(t) + u(t)
-#     x(1) → min
-
-# end
-
-# sol = solve(ocp, grid_size=20, print_level=5)
-# plt = plot(sol)
-
-# a la main
-using CTBase
-using PlotUtils
-
-n = 1
-m = 1
-t0 = 0
-tf = 1
-x0 = 0
-x = t -> -1 + 1e-8 * rand()
-p = t -> 0 + 1e-8 * rand()
-u = t -> 0
-objective = 1
-#
-N = 201
-times = range(t0, tf, N)
-#
-
-sol = OptimalControlSolution()
-sol.state_dimension = n
-sol.control_dimension = m
-sol.time_grid = times
-sol.time_name = "t"
-sol.state = x
-sol.state_name = "x"
-sol.state_components_names = ["x"]
-sol.costate = p
-sol.control = u
-sol.control_name = "u"
-sol.control_components_names = ["u"]
-sol.objective = objective
-sol.iterations = 0
-sol.stopping = :dummy
-sol.message = "ceci est un test"
-sol.success = true
-
-plt = plot(sol)
diff --git a/test/test_plot_manual.jl b/test/test_plot_manual.jl
deleted file mode 100644
index a5ce33a..0000000
--- a/test/test_plot_manual.jl
+++ /dev/null
@@ -1,211 +0,0 @@
-using CTBase
-using Plots
-
-layout = :split
-size = (900, 600)
-control_plt = :all
-
-#
-do_plot_2 = true
-
-# ----------------------------------------
-# SOL 1
-n = 2
-m = 1
-t0 = 0.0
-tf = 1.0
-x0 = [-1.0, 0.0]
-xf = [0.0, 0.0]
-a = x0[1]
-b = x0[2]
-C = [
-    -(tf - t0)^3/6.0 (tf - t0)^2/2.0
-    -(tf - t0)^2/2.0 (tf-t0)
-]
-D = [-a - b * (tf - t0), -b] + xf
-p0 = C \ D
-α = p0[1]
-β = p0[2]
-x =
-    t -> [
-        a + b * (t - t0) + β * (t - t0)^2 / 2.0 - α * (t - t0)^3 / 6.0,
-        b + β * (t - t0) - α * (t - t0)^2 / 2.0,
-    ]
-p = t -> [α, -α * (t - t0) + β]
-u = t -> [p(t)[2]]
-objective = 0.5 * (α^2 * (tf - t0)^3 / 3 + β^2 * (tf - t0) - α * β * (tf - t0)^2)
-#
-N = 201
-times = range(t0, tf, N)
-#
-
-sol = OptimalControlSolution()
-sol.state_dimension = n
-sol.control_dimension = m
-sol.time_grid = times
-sol.time_name = "t"
-sol.state = x
-sol.state_name = "x"
-sol.state_components_names = ["x" * ctindices(i) for i ∈ range(1, n)]
-sol.costate = p
-sol.control = u
-sol.control_name = "u"
-sol.control_components_names = ["u"]
-sol.objective = objective
-sol.iterations = 0
-sol.stopping = :dummy
-sol.message = "ceci est un test"
-sol.success = true
-
-#
-plt = plot(
-    sol,
-    layout = layout,
-    control = control_plt,
-    size = size,
-    flip = true,
-    linewidth = 5,
-    solution_label = "sol1",
-)
-#plot(sol, layout=:group)
-#ps=plot(sol, :time, (:state, 1))
-#plot!(ps, sol, :time, (:control, 1))
-
-# ----------------------------------------
-# SOL 2
-n = 2
-m = 1
-t0 = 0.0
-tf = 1.0
-x0 = [-1.0, -1.0]
-xf = [0.0, 0.0]
-a = x0[1]
-b = x0[2]
-C = [
-    -(tf - t0)^3/6.0 (tf - t0)^2/2.0
-    -(tf - t0)^2/2.0 (tf-t0)
-]
-D = [-a - b * (tf - t0), -b] + xf
-p0 = C \ D
-α = p0[1]
-β = p0[2]
-x =
-    t -> [
-        a + b * (t - t0) + β * (t - t0)^2 / 2.0 - α * (t - t0)^3 / 6.0,
-        b + β * (t - t0) - α * (t - t0)^2 / 2.0,
-    ]
-p = t -> [α, -α * (t - t0) + β]
-u = t -> [p(t)[2]]
-objective = 0.5 * (α^2 * (tf - t0)^3 / 3 + β^2 * (tf - t0) - α * β * (tf - t0)^2)
-#
-N = 201
-times = range(t0, tf, N)
-#
-
-sol = OptimalControlSolution()
-sol.state_dimension = n
-sol.control_dimension = m
-sol.time_grid = times
-sol.time_name = "t"
-sol.state = x
-sol.state_name = "y"
-sol.state_components_names = ["y" * ctindices(i) for i ∈ range(1, n)]
-sol.costate = p
-sol.control = u
-sol.control_name = "v"
-sol.control_components_names = ["v"]
-sol.objective = objective
-sol.iterations = 0
-sol.stopping = :dummy
-sol.message = "ceci est un test"
-sol.success = true
-
-if do_plot_2
-    plot!(plt, sol, layout = layout, size = size, control = control_plt, solution_label = "sol2")
-else
-    plt
-end
-
-# ----------------------------------------
-# Orbital transfer consumption
-#=
-n=4
-m=2
-
-x0 = [-42272.67, 0, 0, -5796.72]
-μ      = 5.1658620912*1e12
-rf     = 42165.0 ;
-rf3    = rf^3  ;
-m0     = 2000.0
-F_max = 100.0
-γ_max  = F_max*3600.0^2/(m0*10^3)
-t0     = 0.0
-α      = sqrt(μ/rf3);
-β      = 0.0
-
-tol    = 1e-9;
-
-F_max_100  = 100.0
-
-tf_min = 13.40318195708344 # minimal time for Fmax = 100 N
-tf = 1.5*tf_min
-
-Th(F) = F*3600.0^2/(10^3)
-u_max = Th(F_max)
-
-# the solution
-x0 = [x0; 0]
-
-u0(x,p) = [0, 0]
-u1(x,p) = p[3:4]/norm(p[3:4])
-
-Hc(x,p) = p[1]*x[3] + p[2]*x[4] + p[3]*(-μ*x[1]/norm(x[1:2])^3) + p[4]*(-μ*x[2]/norm(x[1:2])^3)
-H(x,p,u) = -norm(u) + Hc(x,p) + u[1]*p[3]*γ_max + u[2]*p[4]*γ_max + p[5]*norm(u)
-H0(x,p) = H(x,p,u0(x,p)) 
-H1(x,p) = H(x,p,u1(x,p))
-
-# Flow
-f0 = Flow(Hamiltonian(H0));
-f1 = Flow(Hamiltonian(H1));
-
-# Initial guess
-p0_guess = [0.02698412111231433, 0.006910835140705538, 0.050397371862031096, -0.0032972040120747836, -1.0076835239866583e-23]
-ti_guess = [0.4556797711668658, 3.6289692721936913, 11.683607683450061, 12.505465498856514]
-ξ_guess  = [p0_guess;ti_guess]
-
-p0 = p0_guess
-t1, t2, t3, t4 = ti_guess
-
-# computing x, p, u
-f = f1 * (t1, f0) * (t2, f1) * (t3, f0) * (t4, f1)
-ode_sol  = f((t0, tf), x0, p0)
-
-x = t -> ode_sol(t)[1:4]
-p = t -> ode_sol(t)[6:9]
-u = t -> [0,0]*(t ∈ Interval(t1,t2)∪Interval(t3,t4)) +
-         p(t)[3:4]/norm(p(t)[3:4])*(t ∈ Interval(t0,t1)∪Interval(t2,t3)∪Interval(t4,tf))
-objective = ode_sol(tf)[5]
-
-#
-N=201
-times = range(t0, tf, N)
-#
-sol = OptimalControlSolution() #n, m, times, x, p, u)
-sol.state_dimension = n
-sol.control_dimension = m
-sol.time_grid = times
-sol.state = x
-sol.state_name = "y"
-sol.state_components_names = [ "x" * ctindices(1), "x" * ctindices(2), "v" * ctindices(1), "v" * ctindices(2)]
-sol.adjoint = p
-sol.control = u
-sol.control_name = "u"
-sol.control_components_names = [ "u" * ctindices(i) for i ∈ range(1, m)]
-sol.objective = objective
-sol.iterations = 0
-sol.message = "structure: B+B0B+B0B+"
-sol.success = true
-sol.infos[:resolution] = :numerical
-
-plt_transfert = plot(sol, layout=:split, size=(900, 600))
-=#