diff --git a/Project.toml b/Project.toml
index 9e939875..80c2ce1c 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,6 +5,7 @@ version = "0.2.1"
 
 [deps]
 CTBase = "54762871-cc72-4466-b8e8-f6c8b58076cd"
+CTDirect = "790bbbee-bee9-49ee-8912-a9de031322d5"
 CTModels = "34c4fa32-2049-4079-8329-de33c2a22e2d"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2"
@@ -20,6 +21,7 @@ OptimalControlModels = "OptimalControl"
 
 [compat]
 CTBase = "0.16"
+CTDirect = "0.16"
 CTModels = "0.6"
 DocStringExtensions = "0.9"
 ExaModels = "0.9"
diff --git a/docs/make.jl b/docs/make.jl
index ec58b8b7..52603e78 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -69,7 +69,7 @@ cp(
 repo_url = "github.com/control-toolbox/OptimalControlProblems.jl"
 
 #
-draft = true
+draft = false
 exclude_from_draft=Symbol[
 #    :beam
 ]
diff --git a/docs/problems.jl b/docs/problems.jl
index 89c32973..ba61045c 100644
--- a/docs/problems.jl
+++ b/docs/problems.jl
@@ -63,11 +63,12 @@ function generate_documentation(
 
     ```@example main
     # import model
-    docp, model_oc = $PROBLEM(OptimalControlBackend())
+    docp = $PROBLEM(OptimalControlBackend())
+    nlp_oc = nlp_model(docp)
 
     # solve
     nlp_sol = NLPModelsIpopt.ipopt(
-        model_oc;
+        nlp_oc;
         print_level=4,
         tol=1e-8,
         mu_strategy="adaptive",
@@ -100,8 +101,8 @@ function generate_documentation(
         (
             Problem=:$PROBLEM,
             Grid_Size=OptimalControlProblems.metadata[:$PROBLEM][:N],
-            Variables=get_nvar(model_oc),
-            Constraints=get_ncon(model_oc),
+            Variables=get_nvar(nlp_oc),
+            Constraints=get_ncon(nlp_oc),
         )
     )
     ```
@@ -155,27 +156,27 @@ function generate_documentation(
 
     ```@example main
     # import model
-    model_jp = $PROBLEM(JuMPBackend())
+    nlp_jp = $PROBLEM(JuMPBackend())
 
     # solve
-    set_optimizer(model_jp, Ipopt.Optimizer)
-    set_optimizer_attribute(model_jp, "print_level", 4)
-    set_optimizer_attribute(model_jp, "tol", 1e-8)
-    set_optimizer_attribute(model_jp, "mu_strategy", "adaptive")
-    set_optimizer_attribute(model_jp, "linear_solver", "mumps")
-    set_optimizer_attribute(model_jp, "sb", "yes")
-    optimize!(model_jp)
+    set_optimizer(nlp_jp, Ipopt.Optimizer)
+    set_optimizer_attribute(nlp_jp, "print_level", 4)
+    set_optimizer_attribute(nlp_jp, "tol", 1e-8)
+    set_optimizer_attribute(nlp_jp, "mu_strategy", "adaptive")
+    set_optimizer_attribute(nlp_jp, "linear_solver", "mumps")
+    set_optimizer_attribute(nlp_jp, "sb", "yes")
+    optimize!(nlp_jp)
     ```
 
     For the numerical comparison, define:
 
     ```@example main
-    t_jp = time_grid(:$PROBLEM, model_jp)
-    x_jp = state(:$PROBLEM, model_jp).(t_jp)
-    u_jp = control(:$PROBLEM, model_jp).(t_jp)
-    o_jp = objective_value(model_jp)
-    v_jp = variable(:$PROBLEM, model_jp)
-    i_jp = barrier_iterations(model_jp)
+    t_jp = time_grid(:$PROBLEM, nlp_jp)
+    x_jp = state(:$PROBLEM, nlp_jp).(t_jp)
+    u_jp = control(:$PROBLEM, nlp_jp).(t_jp)
+    o_jp = objective_value(nlp_jp)
+    v_jp = variable(:$PROBLEM, nlp_jp)
+    i_jp = barrier_iterations(nlp_jp)
     nothing # hide
     ```
 
@@ -187,9 +188,9 @@ function generate_documentation(
     push!(data_re,
         (
             Model=:JuMP,
-            Flag=termination_status(model_jp),
-            Iterations=barrier_iterations(model_jp),
-            Objective=objective_value(model_jp),
+            Flag=termination_status(nlp_jp),
+            Iterations=barrier_iterations(nlp_jp),
+            Objective=objective_value(nlp_jp),
         )
     )
     ```
@@ -199,10 +200,10 @@ function generate_documentation(
     Overlay the JuMP solution on the previous plots:
 
     ```@example main
-    t = time_grid(:$PROBLEM, model_jp)     # t0, ..., tN = tf
-    x = state(:$PROBLEM, model_jp)         # function of time
-    u = control(:$PROBLEM, model_jp)       # function of time
-    p = costate(:$PROBLEM, model_jp)       # function of time
+    t = time_grid(:$PROBLEM, nlp_jp)     # t0, ..., tN = tf
+    x = state(:$PROBLEM, nlp_jp)         # function of time
+    u = control(:$PROBLEM, nlp_jp)       # function of time
+    p = costate(:$PROBLEM, nlp_jp)       # function of time
 
     for i in 1:n # state
         label = i == 1 ? "JuMP" : :none
@@ -231,11 +232,12 @@ function generate_documentation(
     function plot_initial_guess()
 
         # import OptimalControl model
-        docp, model_oc = $PROBLEM(OptimalControlBackend())
+        docp = $PROBLEM(OptimalControlBackend())
+        nlp_oc = nlp_model(docp)
 
         # solve
         nlp_sol = NLPModelsIpopt.ipopt(
-            model_oc;
+            nlp_oc;
             max_iter=0,
             print_level=5,
             tol=1e-8,
@@ -263,23 +265,23 @@ function generate_documentation(
         end
 
         # import JuMP model
-        model_jp = $PROBLEM(JuMPBackend())
+        nlp_jp = $PROBLEM(JuMPBackend())
 
         # solve
-        set_optimizer(model_jp, Ipopt.Optimizer)
-        set_optimizer_attribute(model_jp, "max_iter", 0)
-        set_optimizer_attribute(model_jp, "print_level", 5)
-        set_optimizer_attribute(model_jp, "tol", 1e-8)
-        set_optimizer_attribute(model_jp, "mu_strategy", "adaptive")
-        set_optimizer_attribute(model_jp, "linear_solver", "mumps")
-        set_optimizer_attribute(model_jp, "sb", "yes")
-        optimize!(model_jp)
+        set_optimizer(nlp_jp, Ipopt.Optimizer)
+        set_optimizer_attribute(nlp_jp, "max_iter", 0)
+        set_optimizer_attribute(nlp_jp, "print_level", 5)
+        set_optimizer_attribute(nlp_jp, "tol", 1e-8)
+        set_optimizer_attribute(nlp_jp, "mu_strategy", "adaptive")
+        set_optimizer_attribute(nlp_jp, "linear_solver", "mumps")
+        set_optimizer_attribute(nlp_jp, "sb", "yes")
+        optimize!(nlp_jp)
 
         # plot
-        t = time_grid(:$PROBLEM, model_jp)     # t0, ..., tN = tf
-        x = state(:$PROBLEM, model_jp)         # function of time
-        u = control(:$PROBLEM, model_jp)       # function of time
-        p = costate(:$PROBLEM, model_jp)       # function of time
+        t = time_grid(:$PROBLEM, nlp_jp)     # t0, ..., tN = tf
+        x = state(:$PROBLEM, nlp_jp)         # function of time
+        u = control(:$PROBLEM, nlp_jp)       # function of time
+        p = costate(:$PROBLEM, nlp_jp)       # function of time
 
         for i in 1:n # state
             label = i == 1 ? "JuMP" : :none
diff --git a/docs/src/assets/Manifest.toml b/docs/src/assets/Manifest.toml
index 988215c7..fbc589ca 100644
--- a/docs/src/assets/Manifest.toml
+++ b/docs/src/assets/Manifest.toml
@@ -1042,7 +1042,7 @@ uuid = "5f98b655-cc9a-415a-b60e-744165666948"
 version = "1.1.1"
 
 [[deps.OptimalControlProblems]]
-deps = ["CTBase", "CTModels", "DocStringExtensions", "ExaModels", "OrderedCollections"]
+deps = ["CTBase", "CTDirect", "CTModels", "DocStringExtensions", "ExaModels", "OrderedCollections"]
 path = "/Users/ocots/Research/logiciels/dev/control-toolbox/OptimalControlProblems"
 uuid = "59046045-fb9c-4c23-964f-ff0a25704f96"
 version = "0.2.1"
diff --git a/docs/src/dev-add.md b/docs/src/dev-add.md
index dc48c437..eb31e7b8 100644
--- a/docs/src/dev-add.md
+++ b/docs/src/dev-add.md
@@ -42,9 +42,9 @@ function OptimalControlProblems.new_problem(::OptimalControlBackend; N::Int=step
     init = () 
 
     # DOCP and NLP
-    docp, nlp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
+    docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
 
-    return docp, nlp
+    return docp
 
 end
 ```
@@ -61,11 +61,11 @@ function OptimalControlProblems.new_problem(::JuMPBackend; N::Int=steps_number_d
     tf = final_time_data(:new_problem)
 
     # model
-    model = JuMP.Model()
+    nlp = JuMP.Model()
 
     # Define the problem here
     # ...
 
-    return model
+    return nlp
 end
 ```
diff --git a/docs/src/problems-introduction.md b/docs/src/problems-introduction.md
index 63fadd52..049e2b1d 100644
--- a/docs/src/problems-introduction.md
+++ b/docs/src/problems-introduction.md
@@ -136,8 +136,9 @@ data_nlp = DataFrame(           # to store data of the NLPs
 for problem in problems()
 
     #
-    docp, model = eval(problem)(OptimalControlBackend())
-    ocp = docp.ocp
+    docp = eval(problem)(OptimalControlBackend())
+    nlp = nlp_model(docp)
+    ocp = ocp_model(docp)
 
     #
     cost = if has_mayer_cost(ocp) && has_lagrange_cost(ocp)
@@ -188,8 +189,8 @@ for problem in problems()
         (
             Problem=problem,
             Steps=N,
-            Variables=get_nvar(model),
-            Constraints=get_ncon(model),
+            Variables=get_nvar(nlp),
+            Constraints=get_ncon(nlp),
         )
     )
 end
diff --git a/docs/src/tutorial-get.md b/docs/src/tutorial-get.md
index 067db642..8d6a39f4 100644
--- a/docs/src/tutorial-get.md
+++ b/docs/src/tutorial-get.md
@@ -4,7 +4,7 @@ Each problem in the **OptimalControlProblems** package is modelled both in JuMP
 
 ## Get an OptimalControl model
 
-### NLP Model
+### DOCP and NLP Model
 
 To get an OptimalControl model, first install [OptimalControl](https://control-toolbox.org/OptimalControl.jl/stable/#Installation) and import the packages:
 
@@ -16,44 +16,48 @@ using OptimalControlProblems
 Then, to obtain the OptimalControl model of the beam problem, run:
 
 ```@example main_oc
-_, model = beam(OptimalControlBackend())
-model # hide
+docp = beam(OptimalControlBackend())
+nlp = nlp_model(docp)
 ```
 
-The model represents the nonlinear programming problem (NLP) obtained after discretising the optimal control problem (OCP). See the [Introduction](@ref problems-introduction) page for details. The model is an [`ADNLPModels.ADNLPModel`](@extref), which provides automatic differentiation (AD)-based models that follow the [NLPModels.jl](https://github.com/JuliaSmoothOptimizers/NLPModels.jl) API.
+The `nlp` model represents the nonlinear programming problem (NLP) obtained after discretising the optimal control problem (OCP). See the [Introduction](@ref problems-introduction) page for details. The model is an [`ADNLPModels.ADNLPModel`](@extref), which provides automatic differentiation (AD)-based models that follow the [NLPModels.jl](https://github.com/JuliaSmoothOptimizers/NLPModels.jl) API.
 
-### DOCP
+!!! note
+
+    You also have access to the DOCP model, which corresponds to the discretised optimal control problem. Roughly speaking, the DOCP model is the union of the NLP and OCP models. For more details, see this [tutorial](@extref Tutorials Discretization-and-NLP-problem) or the documentation of [`CTDirect.DOCP`](@extref).
+
+### OCP
 
-You also have access to the DOCP model, which corresponds to the discretised optimal control problem. For more details, see this [tutorial](@extref Tutorials Discretization-and-NLP-problem) or the documentation of [`CTDirect.DOCP`](@extref).
+You also have access to the OCP model, which corresponds to the optimal control problem.
 
 ```@example main_oc
-docp, model = beam(OptimalControlBackend())
+ocp = ocp_model(docp)
 nothing # hide
 ```
 
 ### Number of variables, constraints, and nonzeros
 
-The model follows the [NLPModels.jl](https://github.com/JuliaSmoothOptimizers/NLPModels.jl) API. See the existing [Attributes](https://jso.dev/NLPModels.jl/stable/#Attributes) and the available getter functions (`get_X`) [here](https://jso.dev/NLPModels.jl/stable/reference).  
+The `nlp` model follows the [NLPModels.jl](https://github.com/JuliaSmoothOptimizers/NLPModels.jl) API. See the existing [Attributes](https://jso.dev/NLPModels.jl/stable/#Attributes) and the available getter functions (`get_X`) [here](https://jso.dev/NLPModels.jl/stable/reference).  
 
 To get the number of variables:
 
 ```@example main_oc
 using NLPModels
-get_nvar(model)
+get_nvar(nlp)
 ```
 
 To get the number of constraints:
 
 ```@example main_oc
-get_ncon(model)
+get_ncon(nlp)
 ```
 
 To get the number of nonzeros:
 
 ```@example main_oc
-nnzo = get_nnzo(model) # Gradient of the objective
-nnzj = get_nnzj(model) # Jacobian of the constraints
-nnzh = get_nnzh(model) # Hessian of the Lagrangian
+nnzo = get_nnzo(nlp) # Gradient of the objective
+nnzj = get_nnzj(nlp) # Jacobian of the constraints
+nnzh = get_nnzh(nlp) # Hessian of the Lagrangian
 
 println("nnzo = ", nnzo)
 println("nnzj = ", nnzj)
@@ -77,13 +81,15 @@ OptimalControlProblems.metadata[:beam][:N]
 Each problem can be parameterised by the number of steps:
 
 ```@example main_oc
-docp, model = beam(OptimalControlBackend(); N=100)
-get_nvar(model)
+docp = beam(OptimalControlBackend(); N=100)
+nlp = nlp_model(docp)
+get_nvar(nlp)
 ```
 
 ```@example main_oc
-docp, model = beam(OptimalControlBackend(); N=200)
-get_nvar(model)
+docp = beam(OptimalControlBackend(); N=200)
+nlp = nlp_model(docp)
+get_nvar(nlp)
 ```
 
 ## Get a JuMP model
@@ -98,7 +104,7 @@ using OptimalControlProblems
 Then, to obtain the JuMP model of the beam problem, run:
 
 ```@example main_jp
-model = beam(JuMPBackend())
+nlp = beam(JuMPBackend())
 ```
 
 !!! note
diff --git a/docs/src/tutorial-solve.md b/docs/src/tutorial-solve.md
index 1ce448f9..84657786 100644
--- a/docs/src/tutorial-solve.md
+++ b/docs/src/tutorial-solve.md
@@ -12,11 +12,12 @@ First, import the problem:
 
 ```@example main
 using OptimalControl
-docp, model = beam(OptimalControlBackend())
+docp = beam(OptimalControlBackend())
+nlp = nlp_model(docp)
 nothing # hide
 ```
 
-The model represents the nonlinear programming problem (NLP) obtained after discretising the optimal control problem (OCP). See the [Introduction](@ref problems-introduction) page for details. The model is an [`ADNLPModels.ADNLPModel`](@extref), which provides automatic differentiation (AD)-based models that follow the [NLPModels.jl](https://github.com/JuliaSmoothOptimizers/NLPModels.jl) API. 
+The `nlp` model represents the nonlinear programming problem (NLP) obtained after discretising the optimal control problem (OCP). See the [Introduction](@ref problems-introduction) page for details. The model is an [`ADNLPModels.ADNLPModel`](@extref), which provides automatic differentiation (AD)-based models that follow the [NLPModels.jl](https://github.com/JuliaSmoothOptimizers/NLPModels.jl) API. 
 
 We can then solve the problem using, for instance, [`NLPModelsIpopt.ipopt`](@extref):
 
@@ -25,7 +26,7 @@ using NLPModelsIpopt
 
 # Solve the model
 sol = NLPModelsIpopt.ipopt(
-    model;
+    nlp;
     print_level=5,
     tol=1e-8,
     mu_strategy="adaptive",
@@ -84,7 +85,7 @@ First, import the JuMP model:
 
 ```@example main
 using JuMP
-model = beam(JuMPBackend())
+nlp = beam(JuMPBackend())
 ```
 
 We can then solve the problem using the [`JuMP.optimize!`](@extref) function:
@@ -93,28 +94,28 @@ We can then solve the problem using the [`JuMP.optimize!`](@extref) function:
 using Ipopt
 
 # Set the optimiser
-set_optimizer(model, Ipopt.Optimizer)
+set_optimizer(nlp, Ipopt.Optimizer)
 
 # Set optimiser attributes
-set_optimizer_attribute(model, "tol", 1e-8)
-set_optimizer_attribute(model, "mu_strategy", "adaptive")
-set_optimizer_attribute(model, "linear_solver", "mumps")
-set_optimizer_attribute(model, "sb", "yes")
+set_optimizer_attribute(nlp, "tol", 1e-8)
+set_optimizer_attribute(nlp, "mu_strategy", "adaptive")
+set_optimizer_attribute(nlp, "linear_solver", "mumps")
+set_optimizer_attribute(nlp, "sb", "yes")
 
 # Solve the model
-optimize!(model)
+optimize!(nlp)
 ```
 
 To get the number of iterations:
 
 ```@example main
-barrier_iterations(model)
+barrier_iterations(nlp)
 ```
 
 To get the objective value:
 
 ```@example main
-objective_value(model)
+objective_value(nlp)
 ```
 
 To get the time grid, state, control, and costate, OptimalControlProblems provides the following getters:
@@ -122,10 +123,10 @@ To get the time grid, state, control, and costate, OptimalControlProblems provid
 ```@example main
 problem = :beam
 
-t = time_grid(problem, model)    # t0, ..., tN = tf
-x = state(problem, model)        # function of time
-u = control(problem, model)      # function of time
-p = costate(problem, model)      # function of time
+t = time_grid(problem, nlp)    # t0, ..., tN = tf
+x = state(problem, nlp)        # function of time
+u = control(problem, nlp)      # function of time
+p = costate(problem, nlp)      # function of time
 
 tf = t[end]
 println("tf = ", tf)
@@ -138,7 +139,7 @@ println("p(tf) = ", p(tf))
     If the problem includes additional optimisation variables, such as the final time, you can retrieve them with:
 
     ```julia
-    v = variable(problem, model)
+    v = variable(problem, nlp)
     ```
 
 We can add the state, costate, and control to the plot:
diff --git a/ext/Descriptions/quadrotor.md b/ext/Descriptions/quadrotor.md
index b83edd68..517946d2 100644
--- a/ext/Descriptions/quadrotor.md
+++ b/ext/Descriptions/quadrotor.md
@@ -1,5 +1,5 @@
 This problem models the **minimum-time maneuvering of a quadrotor** from an initial position to a target position while respecting thrust and tilt constraints.  
-The quadrotor is modeled as a rigid body with translational and rotational dynamics, subject to actuator limits and tilt angle restrictions.
+The quadrotor is modelled as a rigid body with translational and rotational dynamics, subject to actuator limits and tilt angle restrictions.
 
 ### System Description
 
diff --git a/ext/OptimalControlModels/beam.jl b/ext/OptimalControlModels/beam.jl
index 70238c0b..1d1675ca 100644
--- a/ext/OptimalControlModels/beam.jl
+++ b/ext/OptimalControlModels/beam.jl
@@ -20,13 +20,7 @@ It then performs direct transcription to generate a discrete optimal control pro
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.beam(OptimalControlBackend(); N=100);
-
-julia> docp
-DOCP object with 100 discretisation points
-
-julia> nlp
-NLP model corresponding to the beam problem
+julia> docp = OptimalControlProblems.beam(OptimalControlBackend(); N=100);
 ```
 
 # References
@@ -56,6 +50,6 @@ function OptimalControlProblems.beam(::OptimalControlBackend; N::Int=steps_numbe
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
-    return docp, nlp
+
+    return docp
 end
diff --git a/ext/OptimalControlModels/bioreactor.jl b/ext/OptimalControlModels/bioreactor.jl
index 6244bc0c..e8935933 100644
--- a/ext/OptimalControlModels/bioreactor.jl
+++ b/ext/OptimalControlModels/bioreactor.jl
@@ -20,13 +20,7 @@ It performs direct transcription to create a discretised optimal control problem
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.bioreactor(OptimalControlBackend(); N=100);
-
-julia> docp
-DOCP object with 100 discretisation points
-
-julia> nlp
-NLP model corresponding to the Bioreactor problem
+julia> docp = OptimalControlProblems.bioreactor(OptimalControlBackend(); N=100);
 ```
 
 # References
@@ -94,6 +88,6 @@ function OptimalControlProblems.bioreactor(::OptimalControlBackend; N::Int=steps
 
     # NLPModel + DOCP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
-    return docp, nlp
+
+    return docp
 end
diff --git a/ext/OptimalControlModels/cart_pendulum.jl b/ext/OptimalControlModels/cart_pendulum.jl
index 37cb53a5..612790ff 100644
--- a/ext/OptimalControlModels/cart_pendulum.jl
+++ b/ext/OptimalControlModels/cart_pendulum.jl
@@ -20,13 +20,7 @@ It performs direct transcription to produce a discretised optimal control proble
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.cart_pendulum(OptimalControlBackend(); N=100);
-
-julia> docp
-DOCP object with 100 discretisation points
-
-julia> nlp
-NLP model corresponding to the Cart-Pendulum problem
+julia> docp = OptimalControlProblems.cart_pendulum(OptimalControlBackend(); N=100);
 ```
 
 # References
@@ -107,6 +101,6 @@ function OptimalControlProblems.cart_pendulum(::OptimalControlBackend; N::Int=st
 
     # NLPModel + DOCP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
-    return docp, nlp
+
+    return docp
 end
diff --git a/ext/OptimalControlModels/chain.jl b/ext/OptimalControlModels/chain.jl
index d7510e14..e2f2a8f6 100644
--- a/ext/OptimalControlModels/chain.jl
+++ b/ext/OptimalControlModels/chain.jl
@@ -20,13 +20,7 @@ It performs direct transcription to produce a discretised optimal control proble
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.chain(OptimalControlBackend(); N=100);
-
-julia> docp
-DOCP object with 100 discretisation points
-
-julia> nlp
-NLP model corresponding to the Hanging Chain problem
+julia> docp = OptimalControlProblems.chain(OptimalControlBackend(); N=100);
 ```
 
 # References
@@ -85,6 +79,6 @@ function OptimalControlProblems.chain(::OptimalControlBackend; N::Int=steps_numb
 
     # NLPModel + DOCP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
-    return docp, nlp
+
+    return docp
 end
diff --git a/ext/OptimalControlModels/dielectrophoretic_particle.jl b/ext/OptimalControlModels/dielectrophoretic_particle.jl
index 15e7a42c..1c1c2460 100644
--- a/ext/OptimalControlModels/dielectrophoretic_particle.jl
+++ b/ext/OptimalControlModels/dielectrophoretic_particle.jl
@@ -20,13 +20,7 @@ It performs direct transcription to produce a discretised optimal control proble
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.dielectrophoretic_particle(OptimalControlBackend(); N=100);
-
-julia> docp
-DOCP object with 100 discretisation points
-
-julia> nlp
-NLP model corresponding to the dielectrophoretic particle problem
+julia> docp = OptimalControlProblems.dielectrophoretic_particle(OptimalControlBackend(); N=100);
 ```
 
 # References
@@ -71,7 +65,6 @@ function OptimalControlProblems.dielectrophoretic_particle(
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/double_oscillator.jl b/ext/OptimalControlModels/double_oscillator.jl
index f4bcbf05..cb02b6f8 100644
--- a/ext/OptimalControlModels/double_oscillator.jl
+++ b/ext/OptimalControlModels/double_oscillator.jl
@@ -20,13 +20,7 @@ It uses direct transcription to produce a discretised optimal control problem (D
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.double_oscillator(OptimalControlBackend(); N=100);
-
-julia> docp
-DOCP object with 100 discretisation points
-
-julia> nlp
-NLP model corresponding to the double oscillator problem
+julia> docp = OptimalControlProblems.double_oscillator(OptimalControlBackend(); N=100);
 ```
 
 # References
@@ -80,7 +74,6 @@ function OptimalControlProblems.double_oscillator(::OptimalControlBackend; N::In
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/ducted_fan.jl b/ext/OptimalControlModels/ducted_fan.jl
index 91906a1d..b7086081 100644
--- a/ext/OptimalControlModels/ducted_fan.jl
+++ b/ext/OptimalControlModels/ducted_fan.jl
@@ -20,13 +20,7 @@ It returns both a discretised direct optimal control problem (DOCP) and the corr
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.ducted_fan(OptimalControlBackend(); N=250);
-
-julia> docp
-DOCP object with 250 discretisation points
-
-julia> nlp
-NLP model corresponding to the planar ducted fan problem
+julia> docp = OptimalControlProblems.ducted_fan(OptimalControlBackend(); N=250);
 ```
 
 # References
@@ -101,7 +95,6 @@ function OptimalControlProblems.ducted_fan(::OptimalControlBackend; N::Int=steps
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/electric_vehicle.jl b/ext/OptimalControlModels/electric_vehicle.jl
index 9ef6ff72..713cbaed 100644
--- a/ext/OptimalControlModels/electric_vehicle.jl
+++ b/ext/OptimalControlModels/electric_vehicle.jl
@@ -20,13 +20,7 @@ It returns both a discretised direct optimal control problem (DOCP) and the corr
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.electric_vehicle(OptimalControlBackend(); N=500);
-
-julia> docp
-DOCP object with 500 discretisation points
-
-julia> nlp
-NLP model corresponding to the electric vehicle problem
+julia> docp = OptimalControlProblems.electric_vehicle(OptimalControlBackend(); N=500);
 ```
 
 # References
@@ -72,7 +66,6 @@ function OptimalControlProblems.electric_vehicle(::OptimalControlBackend; N::Int
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/glider.jl b/ext/OptimalControlModels/glider.jl
index 15ae0a4e..d1109de9 100644
--- a/ext/OptimalControlModels/glider.jl
+++ b/ext/OptimalControlModels/glider.jl
@@ -20,13 +20,7 @@ It returns both a discretised direct optimal control problem (DOCP) and the corr
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.glider(OptimalControlBackend(); N=500);
-
-julia> docp
-DOCP object with 500 discretisation points
-
-julia> nlp
-NLP model corresponding to the hang glider problem
+julia> docp = OptimalControlProblems.glider(OptimalControlBackend(); N=500);
 ```
 
 # References
@@ -112,7 +106,6 @@ function OptimalControlProblems.glider(::OptimalControlBackend; N::Int=steps_num
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/insurance.jl b/ext/OptimalControlModels/insurance.jl
index 28e5701c..47d13865 100644
--- a/ext/OptimalControlModels/insurance.jl
+++ b/ext/OptimalControlModels/insurance.jl
@@ -20,13 +20,7 @@ It returns both a discretised direct optimal control problem (DOCP) and the corr
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.insurance(OptimalControlBackend(); N=500);
-
-julia> docp
-DOCP object with 500 discretisation points
-
-julia> nlp
-NLP model corresponding to the insurance optimisation problem
+julia> docp = OptimalControlProblems.insurance(OptimalControlBackend(); N=500);
 ```
 
 # References
@@ -105,7 +99,6 @@ function OptimalControlProblems.insurance(::OptimalControlBackend; N::Int=steps_
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/jackson.jl b/ext/OptimalControlModels/jackson.jl
index 7af13c6f..bcfb71d8 100644
--- a/ext/OptimalControlModels/jackson.jl
+++ b/ext/OptimalControlModels/jackson.jl
@@ -20,13 +20,7 @@ It returns both a discretised direct optimal control problem (DOCP) and the corr
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.jackson(OptimalControlBackend(); N=500);
-
-julia> docp
-DOCP object with 500 discretisation points
-
-julia> nlp
-NLP model corresponding to the Jackson problem
+julia> docp = OptimalControlProblems.jackson(OptimalControlBackend(); N=500);
 ```
 
 # References
@@ -71,7 +65,6 @@ function OptimalControlProblems.jackson(::OptimalControlBackend; N::Int=steps_nu
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/moonlander.jl b/ext/OptimalControlModels/moonlander.jl
index d082b54b..c603eeae 100644
--- a/ext/OptimalControlModels/moonlander.jl
+++ b/ext/OptimalControlModels/moonlander.jl
@@ -20,13 +20,7 @@ It returns both a discretised direct optimal control problem (DOCP) and the corr
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.moonlander(OptimalControlBackend(); N=500);
-
-julia> docp
-DOCP object with 500 discretisation points
-
-julia> nlp
-NLP model corresponding to the Moonlander problem
+julia> docp = OptimalControlProblems.moonlander(OptimalControlBackend(); N=500);
 ```
 """
 function OptimalControlProblems.moonlander(::OptimalControlBackend; N::Int=steps_number_data(:moonlander))
@@ -102,7 +96,6 @@ function OptimalControlProblems.moonlander(::OptimalControlBackend; N::Int=steps
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/quadrotor.jl b/ext/OptimalControlModels/quadrotor.jl
index fa9fd577..911f3474 100644
--- a/ext/OptimalControlModels/quadrotor.jl
+++ b/ext/OptimalControlModels/quadrotor.jl
@@ -20,13 +20,7 @@ It returns both a discretised direct optimal control problem (DOCP) and the corr
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.quadrotor(OptimalControlBackend(); N=50);
-
-julia> docp
-DOCP object with 50 discretisation points
-
-julia> nlp
-NLP model corresponding to the Quadrotor problem
+julia> docp = OptimalControlProblems.quadrotor(OptimalControlBackend(); N=50);
 ```
 """
 function OptimalControlProblems.quadrotor(::OptimalControlBackend; N::Int=steps_number_data(:quadrotor))
@@ -117,7 +111,6 @@ function OptimalControlProblems.quadrotor(::OptimalControlBackend; N::Int=steps_
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/robbins.jl b/ext/OptimalControlModels/robbins.jl
index 4ea2b498..08a4edc4 100644
--- a/ext/OptimalControlModels/robbins.jl
+++ b/ext/OptimalControlModels/robbins.jl
@@ -21,13 +21,7 @@ Reference: [Robbins Problem on BOCOP](https://github.com/control-toolbox/bocop/t
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.robbins(OptimalControlBackend(); N=500);
-
-julia> docp
-DOCP object with 500 discretisation points
-
-julia> nlp
-NLP model corresponding to the Robbins problem
+julia> docp = OptimalControlProblems.robbins(OptimalControlBackend(); N=500);
 ```
 """
 function OptimalControlProblems.robbins(::OptimalControlBackend; N::Int=steps_number_data(:robbins))
@@ -61,7 +55,6 @@ function OptimalControlProblems.robbins(::OptimalControlBackend; N::Int=steps_nu
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/robot.jl b/ext/OptimalControlModels/robot.jl
index c3371853..22ce02ad 100644
--- a/ext/OptimalControlModels/robot.jl
+++ b/ext/OptimalControlModels/robot.jl
@@ -20,13 +20,7 @@ Reference: Robot arm problem on BOCOP [here](https://github.com/control-toolbox/
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.robot(OptimalControlBackend(); N=250);
-
-julia> docp
-DOCP object with 250 discretisation points
-
-julia> nlp
-NLP model corresponding to the robot arm problem
+julia> docp = OptimalControlProblems.robot(OptimalControlBackend(); N=250);
 ```
 """
 function OptimalControlProblems.robot(::OptimalControlBackend; N::Int=steps_number_data(:robot))
@@ -99,7 +93,6 @@ function OptimalControlProblems.robot(::OptimalControlBackend; N::Int=steps_numb
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/rocket.jl b/ext/OptimalControlModels/rocket.jl
index 569f67c1..ee37b45f 100644
--- a/ext/OptimalControlModels/rocket.jl
+++ b/ext/OptimalControlModels/rocket.jl
@@ -20,13 +20,7 @@ Reference: Goddard Rocket Problem [here](https://github.com/control-toolbox/boco
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.rocket(OptimalControlBackend(); N=500);
-
-julia> docp
-DOCP object with 500 discretisation points
-
-julia> nlp
-NLP model corresponding to the Goddard rocket problem
+julia> docp = OptimalControlProblems.rocket(OptimalControlBackend(); N=500);
 ```
 """
 function OptimalControlProblems.rocket(::OptimalControlBackend; N::Int=steps_number_data(:rocket))
@@ -92,7 +86,6 @@ function OptimalControlProblems.rocket(::OptimalControlBackend; N::Int=steps_num
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/space_shuttle.jl b/ext/OptimalControlModels/space_shuttle.jl
index 7dacf1d6..9c3abe74 100644
--- a/ext/OptimalControlModels/space_shuttle.jl
+++ b/ext/OptimalControlModels/space_shuttle.jl
@@ -21,13 +21,7 @@ Note: No heating limit path constraint is included.
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.space_shuttle(OptimalControlBackend(); N=500);
-
-julia> docp
-DOCP object with 500 discretisation points
-
-julia> nlp
-NLP model corresponding to the Space Shuttle reentry trajectory problem
+julia> docp = OptimalControlProblems.space_shuttle(OptimalControlBackend(); N=500);
 ```
 """
 function OptimalControlProblems.space_shuttle(::OptimalControlBackend; N::Int=steps_number_data(:space_shuttle))
@@ -163,7 +157,6 @@ function OptimalControlProblems.space_shuttle(::OptimalControlBackend; N::Int=st
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/steering.jl b/ext/OptimalControlModels/steering.jl
index d2cbfef3..18ea723c 100644
--- a/ext/OptimalControlModels/steering.jl
+++ b/ext/OptimalControlModels/steering.jl
@@ -20,13 +20,7 @@ The state vector has four components, and the control is a single scalar input.
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.steering(OptimalControlBackend(); N=500);
-
-julia> docp
-DOCP object with 500 discretisation points
-
-julia> nlp
-NLP model corresponding to the particle steering problem
+julia> docp = OptimalControlProblems.steering(OptimalControlBackend(); N=500);
 ```
 """
 function OptimalControlProblems.steering(::OptimalControlBackend; N::Int=steps_number_data(:steering))
@@ -75,7 +69,6 @@ function OptimalControlProblems.steering(::OptimalControlBackend; N::Int=steps_n
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/truck_trailer.jl b/ext/OptimalControlModels/truck_trailer.jl
index 9291a506..99c81136 100644
--- a/ext/OptimalControlModels/truck_trailer.jl
+++ b/ext/OptimalControlModels/truck_trailer.jl
@@ -20,13 +20,7 @@ The problem includes path constraints for articulation angles between the traile
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.truck_trailer(OptimalControlBackend(); N=200);
-
-julia> docp
-DOCP object with 200 discretisation points
-
-julia> nlp
-NLP model corresponding to the truck-trailer parking problem
+julia> docp = OptimalControlProblems.truck_trailer(OptimalControlBackend(); N=200);
 ```
 """
 function OptimalControlProblems.truck_trailer(::OptimalControlBackend; N::Int=steps_number_data(:truck_trailer))
@@ -135,7 +129,6 @@ function OptimalControlProblems.truck_trailer(::OptimalControlBackend; N::Int=st
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/ext/OptimalControlModels/vanderpol.jl b/ext/OptimalControlModels/vanderpol.jl
index c4e43098..3d84a763 100644
--- a/ext/OptimalControlModels/vanderpol.jl
+++ b/ext/OptimalControlModels/vanderpol.jl
@@ -20,13 +20,7 @@ The problem formulation can be found [here](https://github.com/control-toolbox/b
 ```julia-repl
 julia> using OptimalControlProblems
 
-julia> docp, nlp = OptimalControlProblems.vanderpol(OptimalControlBackend(); N=500);
-
-julia> docp
-DOCP object for the Van der Pol problem with 500 discretisation points
-
-julia> nlp
-NLP model corresponding to the Van der Pol control problem
+julia> docp = OptimalControlProblems.vanderpol(OptimalControlBackend(); N=500);
 ```
 """
 function OptimalControlProblems.vanderpol(::OptimalControlBackend; N::Int=steps_number_data(:vanderpol))
@@ -56,7 +50,6 @@ function OptimalControlProblems.vanderpol(::OptimalControlBackend; N::Int=steps_
 
     # DOCP and NLP
     docp = direct_transcription(ocp; init=init, grid_size=N, disc_method=:trapeze)
-    nlp = model(docp)
 
-    return docp, nlp
+    return docp
 end
diff --git a/src/OptimalControlProblems.jl b/src/OptimalControlProblems.jl
index 851ee663..4361b76c 100644
--- a/src/OptimalControlProblems.jl
+++ b/src/OptimalControlProblems.jl
@@ -1,11 +1,22 @@
 module OptimalControlProblems
 
 using CTBase
+using CTDirect
 import CTModels: CTModels, time_grid, state, control, costate
 import ExaModels: ExaModels, variable
 using DocStringExtensions
 using OrderedCollections: OrderedDict
 
+# -----------------
+# SHOULD NO BE HERE
+nlp_model(docp::CTDirect.DOCP) = docp.nlp
+ocp_model(docp::CTDirect.DOCP) = docp.ocp
+
+export nlp_model, ocp_model
+
+# -----------------
+
+
 """
 $(TYPEDEF)
 
diff --git a/test/figures/init/beam.pdf b/test/figures/init/beam.pdf
index a6f22eb7..c6eefb8a 100644
Binary files a/test/figures/init/beam.pdf and b/test/figures/init/beam.pdf differ
diff --git a/test/figures/init/bioreactor.pdf b/test/figures/init/bioreactor.pdf
index 983cdb0d..91453442 100644
Binary files a/test/figures/init/bioreactor.pdf and b/test/figures/init/bioreactor.pdf differ
diff --git a/test/figures/init/cart_pendulum.pdf b/test/figures/init/cart_pendulum.pdf
index fde41947..dd7f82ac 100644
Binary files a/test/figures/init/cart_pendulum.pdf and b/test/figures/init/cart_pendulum.pdf differ
diff --git a/test/figures/init/chain.pdf b/test/figures/init/chain.pdf
index fab21016..b79fb09b 100644
Binary files a/test/figures/init/chain.pdf and b/test/figures/init/chain.pdf differ
diff --git a/test/figures/init/dielectrophoretic_particle.pdf b/test/figures/init/dielectrophoretic_particle.pdf
index 3d36e7bc..35df5ce3 100644
Binary files a/test/figures/init/dielectrophoretic_particle.pdf and b/test/figures/init/dielectrophoretic_particle.pdf differ
diff --git a/test/figures/init/double_oscillator.pdf b/test/figures/init/double_oscillator.pdf
index f1914e3b..e38c7972 100644
Binary files a/test/figures/init/double_oscillator.pdf and b/test/figures/init/double_oscillator.pdf differ
diff --git a/test/figures/init/ducted_fan.pdf b/test/figures/init/ducted_fan.pdf
index c4471106..467c6a83 100644
Binary files a/test/figures/init/ducted_fan.pdf and b/test/figures/init/ducted_fan.pdf differ
diff --git a/test/figures/init/electric_vehicle.pdf b/test/figures/init/electric_vehicle.pdf
index 2d0ca809..1759cc5b 100644
Binary files a/test/figures/init/electric_vehicle.pdf and b/test/figures/init/electric_vehicle.pdf differ
diff --git a/test/figures/init/glider.pdf b/test/figures/init/glider.pdf
index fd3f85cc..4880815e 100644
Binary files a/test/figures/init/glider.pdf and b/test/figures/init/glider.pdf differ
diff --git a/test/figures/init/insurance.pdf b/test/figures/init/insurance.pdf
index d1359885..e991b0e3 100644
Binary files a/test/figures/init/insurance.pdf and b/test/figures/init/insurance.pdf differ
diff --git a/test/figures/init/jackson.pdf b/test/figures/init/jackson.pdf
index f1719b49..f1210f12 100644
Binary files a/test/figures/init/jackson.pdf and b/test/figures/init/jackson.pdf differ
diff --git a/test/figures/init/moonlander.pdf b/test/figures/init/moonlander.pdf
index ee76adc0..b834b6e7 100644
Binary files a/test/figures/init/moonlander.pdf and b/test/figures/init/moonlander.pdf differ
diff --git a/test/figures/init/quadrotor.pdf b/test/figures/init/quadrotor.pdf
index e6602e45..1ddd9c63 100644
Binary files a/test/figures/init/quadrotor.pdf and b/test/figures/init/quadrotor.pdf differ
diff --git a/test/figures/init/robbins.pdf b/test/figures/init/robbins.pdf
index 0f503f10..017c198c 100644
Binary files a/test/figures/init/robbins.pdf and b/test/figures/init/robbins.pdf differ
diff --git a/test/figures/init/robot.pdf b/test/figures/init/robot.pdf
index 25766717..92b276cc 100644
Binary files a/test/figures/init/robot.pdf and b/test/figures/init/robot.pdf differ
diff --git a/test/figures/init/rocket.pdf b/test/figures/init/rocket.pdf
index b433e9f5..f837c62b 100644
Binary files a/test/figures/init/rocket.pdf and b/test/figures/init/rocket.pdf differ
diff --git a/test/figures/init/space_shuttle.pdf b/test/figures/init/space_shuttle.pdf
index 6d8f776e..5066837a 100644
Binary files a/test/figures/init/space_shuttle.pdf and b/test/figures/init/space_shuttle.pdf differ
diff --git a/test/figures/init/steering.pdf b/test/figures/init/steering.pdf
index 43964600..5b807d7b 100644
Binary files a/test/figures/init/steering.pdf and b/test/figures/init/steering.pdf differ
diff --git a/test/figures/init/truck_trailer.pdf b/test/figures/init/truck_trailer.pdf
index a1cac653..52207f74 100644
Binary files a/test/figures/init/truck_trailer.pdf and b/test/figures/init/truck_trailer.pdf differ
diff --git a/test/figures/init/vanderpol.pdf b/test/figures/init/vanderpol.pdf
index a82e5272..4ab4f43b 100644
Binary files a/test/figures/init/vanderpol.pdf and b/test/figures/init/vanderpol.pdf differ
diff --git a/test/figures/solution/beam.pdf b/test/figures/solution/beam.pdf
index e70048e7..95dd1f84 100644
Binary files a/test/figures/solution/beam.pdf and b/test/figures/solution/beam.pdf differ
diff --git a/test/figures/solution/bioreactor.pdf b/test/figures/solution/bioreactor.pdf
index 67c7793f..92afa41d 100644
Binary files a/test/figures/solution/bioreactor.pdf and b/test/figures/solution/bioreactor.pdf differ
diff --git a/test/figures/solution/cart_pendulum.pdf b/test/figures/solution/cart_pendulum.pdf
index 47dc122a..5dacd8ab 100644
Binary files a/test/figures/solution/cart_pendulum.pdf and b/test/figures/solution/cart_pendulum.pdf differ
diff --git a/test/figures/solution/chain.pdf b/test/figures/solution/chain.pdf
index 0103b4f5..e4c7870b 100644
Binary files a/test/figures/solution/chain.pdf and b/test/figures/solution/chain.pdf differ
diff --git a/test/figures/solution/dielectrophoretic_particle.pdf b/test/figures/solution/dielectrophoretic_particle.pdf
index 6b145270..ad48dfb1 100644
Binary files a/test/figures/solution/dielectrophoretic_particle.pdf and b/test/figures/solution/dielectrophoretic_particle.pdf differ
diff --git a/test/figures/solution/double_oscillator.pdf b/test/figures/solution/double_oscillator.pdf
index 182ad6df..7eca11ab 100644
Binary files a/test/figures/solution/double_oscillator.pdf and b/test/figures/solution/double_oscillator.pdf differ
diff --git a/test/figures/solution/ducted_fan.pdf b/test/figures/solution/ducted_fan.pdf
index c2a1bef4..3e5466eb 100644
Binary files a/test/figures/solution/ducted_fan.pdf and b/test/figures/solution/ducted_fan.pdf differ
diff --git a/test/figures/solution/electric_vehicle.pdf b/test/figures/solution/electric_vehicle.pdf
index fe3fdd9b..277b8e65 100644
Binary files a/test/figures/solution/electric_vehicle.pdf and b/test/figures/solution/electric_vehicle.pdf differ
diff --git a/test/figures/solution/glider.pdf b/test/figures/solution/glider.pdf
index ac343608..38612566 100644
Binary files a/test/figures/solution/glider.pdf and b/test/figures/solution/glider.pdf differ
diff --git a/test/figures/solution/insurance.pdf b/test/figures/solution/insurance.pdf
index 8bf206a6..33c0cd5a 100644
Binary files a/test/figures/solution/insurance.pdf and b/test/figures/solution/insurance.pdf differ
diff --git a/test/figures/solution/jackson.pdf b/test/figures/solution/jackson.pdf
index f33b1112..8e12caa2 100644
Binary files a/test/figures/solution/jackson.pdf and b/test/figures/solution/jackson.pdf differ
diff --git a/test/figures/solution/moonlander.pdf b/test/figures/solution/moonlander.pdf
index 3e2a12a1..1baec84f 100644
Binary files a/test/figures/solution/moonlander.pdf and b/test/figures/solution/moonlander.pdf differ
diff --git a/test/figures/solution/quadrotor.pdf b/test/figures/solution/quadrotor.pdf
index 88da7f88..4b3287ca 100644
Binary files a/test/figures/solution/quadrotor.pdf and b/test/figures/solution/quadrotor.pdf differ
diff --git a/test/figures/solution/robbins.pdf b/test/figures/solution/robbins.pdf
index 4db5d228..4143b32b 100644
Binary files a/test/figures/solution/robbins.pdf and b/test/figures/solution/robbins.pdf differ
diff --git a/test/figures/solution/robot.pdf b/test/figures/solution/robot.pdf
index 2cadd14f..dc2bfbd2 100644
Binary files a/test/figures/solution/robot.pdf and b/test/figures/solution/robot.pdf differ
diff --git a/test/figures/solution/rocket.pdf b/test/figures/solution/rocket.pdf
index 3f6c891a..abfd8608 100644
Binary files a/test/figures/solution/rocket.pdf and b/test/figures/solution/rocket.pdf differ
diff --git a/test/figures/solution/space_shuttle.pdf b/test/figures/solution/space_shuttle.pdf
index 148c11ee..2857f17e 100644
Binary files a/test/figures/solution/space_shuttle.pdf and b/test/figures/solution/space_shuttle.pdf differ
diff --git a/test/figures/solution/steering.pdf b/test/figures/solution/steering.pdf
index e9fee2d7..6fb9688f 100644
Binary files a/test/figures/solution/steering.pdf and b/test/figures/solution/steering.pdf differ
diff --git a/test/figures/solution/truck_trailer.pdf b/test/figures/solution/truck_trailer.pdf
index dd529559..b05ac1ed 100644
Binary files a/test/figures/solution/truck_trailer.pdf and b/test/figures/solution/truck_trailer.pdf differ
diff --git a/test/figures/solution/vanderpol.pdf b/test/figures/solution/vanderpol.pdf
index 23563c32..b4287d69 100644
Binary files a/test/figures/solution/vanderpol.pdf and b/test/figures/solution/vanderpol.pdf differ
diff --git a/test/test_JuMP.jl b/test/test_JuMP.jl
index 710076a5..be8a29f6 100644
--- a/test/test_JuMP.jl
+++ b/test/test_JuMP.jl
@@ -12,35 +12,35 @@ function test_JuMP()
             DEBUG && println("│")
 
             # Set up the model
-            model = OptimalControlProblems.eval(f)(JuMPBackend(); N=N)
-            set_optimizer(model, Ipopt.Optimizer)
-            set_silent(model)
-            set_optimizer_attribute(model, "tol", TOL)
-            set_optimizer_attribute(model, "max_iter", MAX_ITER)
-            set_optimizer_attribute(model, "mu_strategy", MU_STRATEGY)
-            set_optimizer_attribute(model, "linear_solver", "mumps")
-            set_optimizer_attribute(model, "max_wall_time", MAX_WALL_TIME)
-            set_optimizer_attribute(model, "sb", SB)
+            nlp = OptimalControlProblems.eval(f)(JuMPBackend(); N=N)
+            set_optimizer(nlp, Ipopt.Optimizer)
+            set_silent(nlp)
+            set_optimizer_attribute(nlp, "tol", TOL)
+            set_optimizer_attribute(nlp, "max_iter", MAX_ITER)
+            set_optimizer_attribute(nlp, "mu_strategy", MU_STRATEGY)
+            set_optimizer_attribute(nlp, "linear_solver", "mumps")
+            set_optimizer_attribute(nlp, "max_wall_time", MAX_WALL_TIME)
+            set_optimizer_attribute(nlp, "sb", SB)
 
             # Solve the model
             DEBUG && println("├─  Solve")
             DEBUG && println("│")
             print("  First solve:  ");
-            @time optimize!(model)
+            @time optimize!(nlp)
             print("  Second solve: ");
-            @time optimize!(model)
+            @time optimize!(nlp)
             DEBUG && println("│")
 
             # Infos
             DEBUG && println("├─  Infos")
             DEBUG && println("│")
-            DEBUG && println("│     termination_status: ", termination_status(model))
-            DEBUG && println("│     objective: ", objective_value(model))
-            DEBUG && println("│     iterations: ", barrier_iterations(model))
+            DEBUG && println("│     termination_status: ", termination_status(nlp))
+            DEBUG && println("│     objective: ", objective_value(nlp))
+            DEBUG && println("│     iterations: ", barrier_iterations(nlp))
             DEBUG && println("│")
 
             # Test
-            res = @my_test_broken termination_status(model) == MOI.LOCALLY_SOLVED
+            res = @my_test_broken termination_status(nlp) == MOI.LOCALLY_SOLVED
             keep_problem = keep_problem && (typeof(res) == Test.Pass)
             DEBUG &&
                 (typeof(res) == Test.Pass) &&
diff --git a/test/test_OptimalControl.jl b/test/test_OptimalControl.jl
index 68da5046..a77bd684 100644
--- a/test/test_OptimalControl.jl
+++ b/test/test_OptimalControl.jl
@@ -21,15 +21,16 @@ function test_OptimalControl()
             DEBUG && println("│")
 
             # Set up the model
-            _, model = OptimalControlProblems.eval(f)(OptimalControlBackend(); N=N)
+            docp = OptimalControlProblems.eval(f)(OptimalControlBackend(); N=N)
+            nlp = nlp_model(docp)
 
             # Solve the model
             DEBUG && println("├─  Solve")
             DEBUG && println("│")
             print("  First solve:  ");
-            @time sol = NLPModelsIpopt.ipopt(model; kwargs...)
+            @time sol = NLPModelsIpopt.ipopt(nlp; kwargs...)
             print("  Second solve: ");
-            @time sol = NLPModelsIpopt.ipopt(model; kwargs...)
+            @time sol = NLPModelsIpopt.ipopt(nlp; kwargs...)
             DEBUG && println("│")
 
             # Infos
diff --git a/test/test_quick.jl b/test/test_quick.jl
index 957fd124..5741e429 100644
--- a/test/test_quick.jl
+++ b/test/test_quick.jl
@@ -26,7 +26,8 @@ function test_quick()
             DEBUG && println("│")
 
             ########## OptimalControl ##########
-            docp, nlp = OptimalControlProblems.eval(f)(OptimalControlBackend(); N=N)
+            docp = OptimalControlProblems.eval(f)(OptimalControlBackend(); N=N)
+            nlp = nlp_model(docp)
             nlp_sol = NLPModelsIpopt.ipopt(nlp; kwargs...)
             sol = build_OCP_solution(
                 docp;
@@ -37,17 +38,17 @@ function test_quick()
             o_oc = objective(sol)
 
             ############### JuMP ###############
-            model = OptimalControlProblems.eval(f)(JuMPBackend(); N=N)
-            set_optimizer(model, Ipopt.Optimizer)
-            set_silent(model)
-            set_optimizer_attribute(model, "tol", TOL)
-            set_optimizer_attribute(model, "max_iter", MAX_ITER)
-            set_optimizer_attribute(model, "mu_strategy", MU_STRATEGY)
-            set_optimizer_attribute(model, "linear_solver", "mumps")
-            set_optimizer_attribute(model, "max_wall_time", MAX_WALL_TIME)
-            set_optimizer_attribute(model, "sb", SB)
-            optimize!(model)
-            o_jp = objective_value(model)
+            nlp = OptimalControlProblems.eval(f)(JuMPBackend(); N=N)
+            set_optimizer(nlp, Ipopt.Optimizer)
+            set_silent(nlp)
+            set_optimizer_attribute(nlp, "tol", TOL)
+            set_optimizer_attribute(nlp, "max_iter", MAX_ITER)
+            set_optimizer_attribute(nlp, "mu_strategy", MU_STRATEGY)
+            set_optimizer_attribute(nlp, "linear_solver", "mumps")
+            set_optimizer_attribute(nlp, "max_wall_time", MAX_WALL_TIME)
+            set_optimizer_attribute(nlp, "sb", SB)
+            optimize!(nlp)
+            o_jp = objective_value(nlp)
 
             ############### TEST ###############
             # objective
diff --git a/test/utils.jl b/test/utils.jl
index ed4b8fec..9a34000d 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -121,7 +121,8 @@ function comparison(; max_iter, test_name)
             ########## OptimalControl ##########
 
             # set up the OptimalControl model 
-            docp, nlp = OptimalControlProblems.eval(f)(OptimalControlBackend(); N=N)
+            docp = OptimalControlProblems.eval(f)(OptimalControlBackend(); N=N)
+            nlp = nlp_model(docp)
 
             # solve the problem
             nlp_sol = NLPModelsIpopt.ipopt(nlp; Options...)
@@ -147,27 +148,27 @@ function comparison(; max_iter, test_name)
             ############### JuMP ###############
 
             # set up the JuMP model
-            model = OptimalControlProblems.eval(f)(JuMPBackend(); N=N)
-            set_optimizer(model, Ipopt.Optimizer)
-            set_silent(model)
-            set_optimizer_attribute(model, "tol", Options[:tol])
-            set_optimizer_attribute(model, "max_iter", Options[:max_iter])
-            set_optimizer_attribute(model, "mu_strategy", Options[:mu_strategy])
-            set_optimizer_attribute(model, "linear_solver", "mumps")
-            set_optimizer_attribute(model, "max_wall_time", Options[:max_wall_time])
-            set_optimizer_attribute(model, "sb", Options[:sb])
+            nlp = OptimalControlProblems.eval(f)(JuMPBackend(); N=N)
+            set_optimizer(nlp, Ipopt.Optimizer)
+            set_silent(nlp)
+            set_optimizer_attribute(nlp, "tol", Options[:tol])
+            set_optimizer_attribute(nlp, "max_iter", Options[:max_iter])
+            set_optimizer_attribute(nlp, "mu_strategy", Options[:mu_strategy])
+            set_optimizer_attribute(nlp, "linear_solver", "mumps")
+            set_optimizer_attribute(nlp, "max_wall_time", Options[:max_wall_time])
+            set_optimizer_attribute(nlp, "sb", Options[:sb])
 
             # solve the model
-            optimize!(model)
+            optimize!(nlp)
 
             # retrieve values
-            t_jp = time_grid(f, model)
-            x_jp = state(f, model).(t_jp)
-            u_jp = control(f, model).(t_jp)
-            o_jp = objective_value(model)
-            i_jp = barrier_iterations(model)
-            p_jp = costate(f, model).(t_jp)
-            v_jp = variable(f, model)
+            t_jp = time_grid(f, nlp)
+            x_jp = state(f, nlp).(t_jp)
+            u_jp = control(f, nlp).(t_jp)
+            o_jp = objective_value(nlp)
+            i_jp = barrier_iterations(nlp)
+            p_jp = costate(f, nlp).(t_jp)
+            v_jp = variable(f, nlp)
 
             ############ TEST ############