diff --git a/benchmark/prof.jl b/benchmark/prof.jl
index b366beb..12de115 100644
--- a/benchmark/prof.jl
+++ b/benchmark/prof.jl
@@ -33,7 +33,7 @@ function init(;in_place, grid_size, disc_method)
 end
 
 
-function test_unit(;test_get=false, test_dyn=false, test_unit_cons=false, test_mayer=false, test_obj=false, test_block=true, test_cons=false, test_trans=false, test_solve=false, warntype=false, jet=false, profile=false, grid_size=100, disc_method=:trapeze, in_place=true)
+function test_unit(;test_get=false, test_dyn=false, test_unit_cons=false, test_mayer=false, test_obj=false, test_block=false, test_cons=false, test_trans=false, test_solve=false, warntype=false, jet=false, profile=false, grid_size=100, disc_method=:trapeze, in_place=true)
     
     # define problem and variables
     prob, docp, xu = init(in_place=in_place, grid_size=grid_size, disc_method=disc_method)
@@ -45,24 +45,18 @@ function test_unit(;test_get=false, test_dyn=false, test_unit_cons=false, test_m
     if test_get
         println("Getters")
         print("t"); @btime CTDirect.get_final_time($xu, $docp)
-        print("v"); @btime CTDirect.get_optim_variable($xu, $docp)
-        print("vv"); @btime CTDirect.get_OCP_variable($xu, $docp)
-        print("x"); @btime CTDirect.get_state_at_time_step($xu, $docp, $docp.dim_NLP_steps)
-        print("xx"); @btime CTDirect.get_OCP_state_at_time_step($xu, $docp, $docp.dim_NLP_steps)        
-        print("u"); @btime CTDirect.get_control_at_time_step($xu, $docp, $docp.dim_NLP_steps) 
+        print("x"); @btime CTDirect.get_OCP_state_at_time_step($xu, $docp, 1)
+        print("u"); @btime CTDirect.get_OCP_control_at_time_step($xu, $docp, 1)
+        print("v"); @btime CTDirect.get_OCP_variable($xu, $docp)
         if warntype
             @code_warntype CTDirect.get_final_time(xu, docp)
-            @code_warntype CTDirect.get_optim_variable(xu, docp)
-            @code_warntype CTDirect.get_state_at_time_step(xu, docp, docp.dim_NLP_steps)
-            @code_warntype CTDirect.get_control_at_time_step(xu, docp, docp.dim_NLP_steps)
             @code_warntype CTDirect.get_time_grid!(xu, docp)
+            @code_warntype CTDirect.get_OCP_state_at_time_step(xu, docp, 1)
+            @code_warntype CTDirect.get_OCP_control_at_time_step(xu, docp, 1)
+            @code_warntype CTDirect.get_OCP_variable(xu, docp)
         end
     end
 
-    t = CTDirect.get_final_time(xu, docp)
-    v = CTDirect.get_optim_variable(xu, docp)
-    x = CTDirect.get_state_at_time_step(xu, docp, docp.dim_NLP_steps)
-    u = CTDirect.get_control_at_time_step(xu, docp, docp.dim_NLP_steps)
     f = similar(xu, docp.dim_NLP_x)
 
     if test_dyn
@@ -107,73 +101,30 @@ function test_unit(;test_get=false, test_dyn=false, test_unit_cons=false, test_m
         nx = docp.dim_NLP_x
         m = docp.dim_NLP_u
         N = docp.dim_NLP_steps
-        x0 = CTDirect.get_state_at_time_step(xu, docp, 1)
-        xx0 = CTDirect.get_OCP_state_at_time_step(xu, docp, 1)
-        xf = CTDirect.get_state_at_time_step(xu, docp, N+1)
-        xxf = CTDirect.get_OCP_state_at_time_step(xu, docp, N+1)
-        v = CTDirect.get_optim_variable(xu, docp)
-        vv = CTDirect.get_OCP_variable(xu, docp)
+        x0 = CTDirect.get_OCP_state_at_time_step(xu, docp, 1)
+        xf = CTDirect.get_OCP_state_at_time_step(xu, docp, N+1)
+        v = CTDirect.get_OCP_variable(xu, docp)
         obj = similar(xu,1)
         
-        #=println("")
-        print("Local Mayer: views for x0/xf and scalar v"); @btime local_mayer($obj, (@view $xu[1:$n]), (@view $xu[($nx + $m) * $N + 1: ($nx + $m) * $N + $n]), $xu[end])
-
-        print("Local Mayer: raw getters for x0/xf and v"); @btime local_mayer($obj, $x0, $xf, $v)
-        print("Local Mayer: getters with scalarization"); @btime local_mayer($obj, $docp._x($x0), $docp._x($xf), $docp._v($v))
-        
-        print("Local Mayer: scal/vec getters"); @btime local_mayer($obj, $xx0, $xxf, $vv)
-        =#
-
+        # local mayer
         println("")
+        print("Local Mayer: views for x0/xf and scalar v"); @btime local_mayer($obj, (@view $xu[1:$n]), (@view $xu[($nx + $m) * $N + 1: ($nx + $m) * $N + $n]), $xu[end]) # OK
+        print("Local Mayer: param scal/vec getters"); @btime local_mayer($obj, $x0, $xf, $v) # OK
+        print("OCP Mayer: param scal/vec getters"); @btime $docp.ocp.mayer($obj, $x0, $xf, $v) # 3 allocs (112)
 
-        #print("OCP Mayer: views for x0/xf and scalar v"); @btime $docp.ocp.mayer($obj, (@view $xu[1:$n]), (@view $xu[($nx + $m) * $N + 1: ($nx + $m) * $N + $n]), $xu[end])
-
-        #print("OCP Mayer: raw getters for x0/xf and v"); @btime $docp.ocp.mayer($obj, $x0, $xf, $v)
-
-        print("OCP Mayer: getters with scalarization"); @btime $docp.ocp.mayer($obj, $docp._x($x0), $docp._x($xf), $docp._v($v))
-
-        print("OCP Mayer: scal/vec getters"); @btime $docp.ocp.mayer($obj, $xx0, $xxf, $vv)
-
-        if warntype
-            #println("code warntype local mayer")
-            #@code_warntype local_mayer(obj, (@view xu[1:n]), (@view xu[(nx + m) * N + 1: (nx + m) * N + n]), xu[end])
-            #println("code warntype end")
-            println("code warntype ocp mayer getters + scalarization")
-            @code_warntype docp.ocp.mayer(obj, docp._x(x0), docp._x(xf), docp._v(v))
-            println("code warntype end")
-            println("code warntype ocp mayer scal/vect getters")
-            @code_warntype docp.ocp.mayer(obj, xx0, xxf, vv)
-            println("code warntype end")
+        warntype && @code_warntype docp.ocp.mayer(obj, x0, xf, v)
+        jet && display(@report_opt docp.ocp.mayer(obj, x0, xf, v))
+        if profile
+            Profile.Allocs.@profile sample_rate=1.0 docp.ocp.mayer(obj, x0, xf, v)
+            results = Profile.Allocs.fetch()
+            PProf.Allocs.pprof()
         end
-        if jet 
-            println("JET ocp mayer getters + scalarization")
-            display(@report_opt docp.ocp.mayer(obj, docp._x(x0), docp._x(xf), docp._v(v)))
-            println("JET end")
-            println("JET ocp mayer scal/vect getters")
-            display(@report_opt docp.ocp.mayer(obj, xx0, xxf, vv))
-            println("JET end")
-        end             
     end
 
     if test_obj
         print("Objective"); @btime CTDirect.DOCP_objective($xu, $docp)
-        print("Objective_param"); @btime CTDirect.DOCP_objective_param($xu, $docp)
-        if warntype 
-            println("code warntype Objective")
-            @code_warntype CTDirect.DOCP_objective(xu, docp)
-            println("code warntype end")
-            println("code warntype Objective_param")
-            @code_warntype CTDirect.DOCP_objective_param(xu, docp)
-            println("code warntype end")
-        end            
-        if jet
-            println("JET Objective")
-            display(@report_opt CTDirect.DOCP_objective(xu, docp))
-            println("JET end")
-            println("JET Objective_param")
-            display(@report_opt CTDirect.DOCP_objective_param(xu, docp))
-            println("JET end")
-        end
+        warntype && @code_warntype CTDirect.DOCP_objective(xu, docp) # quasi OK (inplace/outplace for ocp.mayer return ?)
+        jet && display(@report_opt CTDirect.DOCP_objective(xu, docp))
         if profile
             Profile.Allocs.@profile sample_rate=1.0 CTDirect.DOCP_objective(xu, docp)
             results = Profile.Allocs.fetch()
@@ -182,34 +133,29 @@ function test_unit(;test_get=false, test_dyn=false, test_unit_cons=false, test_m
     end
 
     if test_block
-        print("Constraints block")
-        CTDirect.get_time_grid!(xu, docp)
-        v = CTDirect.get_OCP_variable_param(xu, docp)
-        work = CTDirect.setWorkArray_param(docp, xu, docp.NLP_time_grid, v)
+        CTDirect.get_time_grid!(xu, docp) # type OK
         i = 1
-        @btime CTDirect.setConstraintBlock_param!($docp, $c, $xu, $v, $docp.NLP_time_grid, $i, $work)
-        +warntype
-        +profile
-
+        v = CTDirect.get_OCP_variable(xu, docp)
+        work = CTDirect.setWorkArray(docp, xu, docp.NLP_time_grid, v)
+        print("Constraints block")
+        @btime CTDirect.setConstraintBlock!($docp, $c, $xu, $v, $docp.NLP_time_grid, $i, $work)
+        warntype && @code_warntype CTDirect.setConstraintBlock!(docp, c, xu, v, docp.NLP_time_grid, i, work)
+        jet && display(@report_opt CTDirect.setConstraintBlock!(docp, c, xu, v, docp.NLP_time_grid, i, work))
+        if profile
+            Profile.Allocs.@profile sample_rate=1.0 CTDirect.setConstraintBlock!(docp, c, xu, v, docp.NLP_time_grid, i, work)
+            results = Profile.Allocs.fetch()
+            PProf.Allocs.pprof()
+        end        
     end
 
     # DOCP_constraints
     if test_cons
-        #print("Constraints"); @btime CTDirect.DOCP_constraints!($c, $xu, $docp)
-        print("Constraints param"); @btime CTDirect.DOCP_constraints_param!($c, $xu, $docp)
-        if any(c.==666.666)
-            error("undefined values in constraints ",c)
-        end
-        if warntype 
-            #@code_warntype CTDirect.DOCP_constraints!(c, xu, docp)
-            @code_warntype CTDirect.DOCP_constraints_param!(c, xu, docp)
-        end
-        if jet
-            #display(@report_opt CTDirect.DOCP_constraints!(c, xu, docp))
-            display(@report_opt CTDirect.DOCP_constraints_param!(c, xu, docp))
-        end
+        print("Constraints"); @btime CTDirect.DOCP_constraints!($c, $xu, $docp)
+        any(c.==666.666) && error("undefined values in constraints ",c)
+        warntype && @code_warntype CTDirect.DOCP_constraints!(c, xu, docp)
+        jet && display(@report_opt CTDirect.DOCP_constraints!(c, xu, docp))
         if profile
-            Profile.Allocs.@profile sample_rate=1.0 CTDirect.DOCP_constraints_param!(c, xu, docp)
+            Profile.Allocs.@profile sample_rate=1.0 CTDirect.DOCP_constraints!(c, xu, docp)
             results = Profile.Allocs.fetch()
             PProf.Allocs.pprof()
         end
@@ -231,54 +177,3 @@ function test_unit(;test_get=false, test_dyn=false, test_unit_cons=false, test_m
 
 end
 
-
-
-#= constraints profile
-  Total:           0       6998 (flat, cum) 89.08%
-     95            .          .               ui = get_control_at_time_step(xu, docp, i) 
-     96            .          .            
-     97            .          .               #1. state equation 
-     98            .          .               if i <= docp.dim_NLP_steps 
-     99            .          .                   # more variables 
-    100            .        100                   fi = copy(work) # create new copy, not just a reference 
-    101            .          .                   tip1 = time_grid[i+1] 
-    102            .          .                   xip1 = get_state_at_time_step(xu, docp, i+1) 
-    103            .          .                   uip1 = get_control_at_time_step(xu, docp, i+1) 
-    104            .          .                   if docp.has_inplace 
-    105            .          .                       docp.dynamics_ext(work, tip1, xip1, uip1, v) 
-    106            .          .                   else 
-    107            .          .                       # copy, do not create a new variable ! 
-    108            .       1000                       work[:] = docp.dynamics_ext(tip1, xip1, uip1, v) #+++ jet runtime dispatch here 
-    109            .          .                   end 
-    110            .          .            
-    111            .          .                   # trapeze rule with 'smart' update for dynamics (similar with @.) 
-    112            .        800                   c[offset+1:offset+docp.dim_NLP_x] = xip1 - (xi + 0.5 * (tip1 - ti) * (fi + work)) #+++ jet runtime dispatch here even with explicit index ranges 
-    113            .          .                   offset += docp.dim_NLP_x 
-    114            .          .               end 
-    115            .          .            
-    116            .          .               # 2. path constraints 
-    117            .          .               # Notes on allocations:.= seems similar 
-    118            .          .               if docp.dim_u_cons > 0 
-    119            .          .                   if docp.has_inplace 
-    120            .          .                       docp.control_constraints[2]((@view c[offset+1:offset+docp.dim_u_cons]),ti, docp._u(ui), v) 
-    121            .          .                   else 
-    122            .       1632                       c[offset+1:offset+docp.dim_u_cons] = docp.control_constraints[2](ti, docp._u(ui), v) 
-    123            .          .                   end 
-    124            .          .               end 
-    125            .          .               if docp.dim_x_cons > 0  
-    126            .          .                   if docp.has_inplace 
-    127            .          .                       docp.state_constraints[2]((@view c[offset+docp.dim_u_cons+1:offset+docp.dim_u_cons+docp.dim_x_cons]),ti, docp._x(xi), v) 
-    128            .          .                   else 
-    129            .       1531                       c[offset+docp.dim_u_cons+1:offset+docp.dim_u_cons+docp.dim_x_cons] = docp.state_constraints[2](ti, docp._x(xi), v) 
-    130            .          .                   end 
-    131            .          .               end 
-    132            .          .               if docp.dim_mixed_cons > 0  
-    133            .          .                   if docp.has_inplace 
-    134            .          .                       docp.mixed_constraints[2]((@view c[offset+docp.dim_u_cons+docp.dim_x_cons+1:offset+docp.dim_u_cons+docp.dim_x_cons+docp.dim_mixed_cons]), ti, docp._x(xi), docp._u(ui), v) 
-    135            .          .                   else 
-    136            .       1935                       c[offset+docp.dim_u_cons+docp.dim_x_cons+1:offset+docp.dim_u_cons+docp.dim_x_cons+docp.dim_mixed_cons] = docp.mixed_constraints[2](ti, docp._x(xi), docp._u(ui), v) 
-    137            .          .                   end 
-    138            .          .               end 
-    139            .          .            
-    140            .          .           end 
-    =#
\ No newline at end of file
diff --git a/src/problem.jl b/src/problem.jl
index e4b87bb..5ed1892 100644
--- a/src/problem.jl
+++ b/src/problem.jl
@@ -24,9 +24,6 @@ struct DOCP{T <: Discretization, X <: ScalVect, U <: ScalVect, V <: ScalVect}
     ## OCP
     ocp::OptimalControlModel # remove at some point ?
 
-    # functions
-    dynamics_ext!::Function
-
     # constraints and their bounds
     control_constraints::Any
     state_constraints::Any
@@ -70,7 +67,7 @@ struct DOCP{T <: Discretization, X <: ScalVect, U <: ScalVect, V <: ScalVect}
     dim_OCP_x::Int  # original OCP state
     dim_NLP_steps::Int
     NLP_normalized_time_grid::Vector{Float64}
-    NLP_time_grid::Vector{Any}
+    NLP_time_grid::Vector{Float64}
     dim_NLP_variables::Int
     dim_NLP_constraints::Int
 
@@ -84,10 +81,6 @@ struct DOCP{T <: Discretization, X <: ScalVect, U <: ScalVect, V <: ScalVect}
     discretization::T
 
     # scalar / vector aux function
-    _vec::Function
-    _x::Function
-    _u::Function
-    _v::Function
     _type_x::X
     _type_u::U  
     _type_v::V
@@ -141,14 +134,6 @@ struct DOCP{T <: Discretization, X <: ScalVect, U <: ScalVect, V <: ScalVect}
         dim_OCP_x = ocp.state_dimension
 
         # times
-        if is_free_initial_time || is_free_final_time
-            # time grid will be recomputed at each NLP iteration
-            NLP_time_grid = Vector{Any}(undef, dim_NLP_steps+1)
-        else
-            # compute time grid once for all 
-            NLP_time_grid = @. ocp.initial_time + (NLP_normalized_time_grid * (ocp.final_time - ocp.initial_time))
-        end
-
         # use 2 different variables for value / index (type stability)
         if is_free_initial_time
             index_initial_time = ocp.initial_time
@@ -165,6 +150,14 @@ struct DOCP{T <: Discretization, X <: ScalVect, U <: ScalVect, V <: ScalVect}
             index_final_time = Index(1) # unused
         end
 
+        if is_free_initial_time || is_free_final_time
+            # time grid will be recomputed at each NLP iteration
+            NLP_time_grid = Vector{Float64}(undef, dim_NLP_steps+1)
+        else
+            # compute time grid once for all 
+            NLP_time_grid = @. fixed_initial_time + (NLP_normalized_time_grid * (fixed_final_time - fixed_initial_time))
+        end
+
         # NLP constraints 
         # parse NLP constraints (and initialize dimensions)
         control_constraints,
@@ -187,49 +180,6 @@ struct DOCP{T <: Discretization, X <: ScalVect, U <: ScalVect, V <: ScalVect}
         dim_mixed_cons = dim_mixed_constraints(ocp)
         dim_boundary_cons = dim_boundary_constraints(ocp)
 
-        # encapsulated OCP functions
-        _vec(f::AbstractVector) = f
-        _vec(f::Number) = [f]
-
-
-        # x[1:dim_OCP_x] adds allocations even with a view
-        # returning just x also allocates, although a bit less
-        _x(x) = (dim_OCP_x == 1) ? x[1] : x[1:dim_OCP_x]
-        _u(u) = (dim_NLP_u == 1) ? u[1] : u
-        _v(v) = (dim_NLP_v == 1) ? v[1] : v
-        # NB. defining 2 functions inside the if is not better
-
-
-
-        # extended dynamics with lagrange cost
-        # remove ?
-        if is_inplace
-            if is_lagrange
-                dynamics_ext! = function (f, t, x, u, v)
-                    ocp.dynamics((@view f[1:dim_OCP_x]), t, _x(x), _u(u), _v(v))
-                    ocp.lagrange((@view f[dim_NLP_x:dim_NLP_x]), t, _x(x), _u(u), _v(v))
-                    return
-                end
-
-            else
-                dynamics_ext! = (f, t, x, u, v) -> ocp.dynamics(f, t, _x(x), _u(u), _v(v))
-            end
-        else
-            if is_lagrange
-                dynamics_ext! = function (f, t, x, u, v)
-                    f[1:docp.dim_OCP_x] .= ocp.dynamics(t, _x(x), _u(u), _v(v)) # .= required for scalar case
-                    f[docp.dim_NLP_x] = ocp.lagrange(t, _x(x), _u(u), _v(v))
-                    return
-                end                
-            else
-                dynamics_ext! = function (f, t, x, u, v)
-                    f[:] = ocp.dynamics(t, _x(x), _u(u), _v(v))
-                    return
-                end
-            end
-        end
-
-
         # parameter: discretization method
         if disc_method == "midpoint"
             discretization = CTDirect.Midpoint(dim_NLP_x, dim_NLP_u, dim_NLP_steps)
@@ -259,7 +209,6 @@ struct DOCP{T <: Discretization, X <: ScalVect, U <: ScalVect, V <: ScalVect}
         # call constructor with const fields
         docp = new{typeof(discretization), typeof(_type_x), typeof(_type_u), typeof(_type_v)}(
             ocp,
-            dynamics_ext!,
             control_constraints,
             state_constraints,
             mixed_constraints,
@@ -302,10 +251,6 @@ struct DOCP{T <: Discretization, X <: ScalVect, U <: ScalVect, V <: ScalVect}
             zeros(dim_NLP_constraints),
             zeros(dim_NLP_constraints),
             discretization,
-            _vec,
-            _x,
-            _u,
-            _v,
             _type_x,
             _type_u,
             _type_v
@@ -392,45 +337,7 @@ $(TYPEDSIGNATURES)
 
 Compute the objective for the DOCP problem.
 """
-function DOCP_objective_param(xu, docp::DOCP)
-
-    obj = similar(xu, 1)
-    N = docp.dim_NLP_steps
-
-    # optimization variables
-    v = get_OCP_variable_param(xu, docp)
-
-    # final state is always needed since lagrange cost is there
-    xf = get_OCP_state_at_time_step_param(xu, docp, N+1)
-
-    # mayer cost
-    if docp.is_mayer
-        x0 = get_OCP_state_at_time_step_param(xu, docp, 1)
-        if docp.is_inplace
-            docp.ocp.mayer(obj, x0, xf, v)
-        else
-            obj[1] = docp.ocp.mayer(x0, xf, v)
-        end
-    end
-
-    # lagrange cost
-    if docp.is_lagrange
-        if docp.is_mayer # NB can this actually happen in OCP (cf bolza) ?
-            obj[1] = obj[1] + get_lagrange_state_at_time_step(xu, docp, docp.dim_NLP_steps+1)
-        else
-            obj[1] = get_lagrange_state_at_time_step(xu, docp, docp.dim_NLP_steps+1)
-        end
-    end
-
-    # maximization problem
-    if docp.is_maximization
-        obj[1] = -obj[1]
-    end
-    
-    return obj[1]
-end
-
-function DOCP_objective_OCP(xu, docp::DOCP)
+function DOCP_objective(xu, docp::DOCP)
 
     obj = similar(xu, 1)
     N = docp.dim_NLP_steps
@@ -453,49 +360,10 @@ function DOCP_objective_OCP(xu, docp::DOCP)
 
     # lagrange cost
     if docp.is_lagrange
-        error("get lagrange state")
         if docp.is_mayer # NB can this actually happen in OCP (cf bolza) ?
-            obj[1] = obj[1] + xf[end]
-        else
-            obj[1] = xf[end]
-        end
-    end
-
-    # maximization problem
-    if docp.is_maximization
-        obj[1] = -obj[1]
-    end
-    
-    return obj[1]
-end
-
-function DOCP_objective(xu, docp::DOCP)
-
-    obj = similar(xu, 1)
-    N = docp.dim_NLP_steps
-
-    # optimization variables
-    v = get_optim_variable(xu, docp)
-
-    # final state is always needed since lagrange cost is there
-    xf = get_state_at_time_step(xu, docp, N+1)
-
-    # mayer cost
-    if docp.is_mayer
-        x0 = get_state_at_time_step(xu, docp, 1)
-        if docp.is_inplace
-            docp.ocp.mayer(obj, docp._x(x0), docp._x(xf), docp._v(v))
-        else
-            obj[1] = docp.ocp.mayer(docp._x(x0), docp._x(xf), docp._v(v))
-        end
-    end
-
-    # lagrange cost
-    if docp.is_lagrange
-        if docp.is_mayer # NB can this actually happen in OCP (cf bolza) ?
-            obj[1] = obj[1] + xf[end]
+            obj[1] = obj[1] + get_lagrange_state_at_time_step(xu, docp, docp.dim_NLP_steps+1)
         else
-            obj[1] = xf[end]
+            obj[1] = get_lagrange_state_at_time_step(xu, docp, docp.dim_NLP_steps+1)
         end
     end
 
@@ -513,41 +381,16 @@ $(TYPEDSIGNATURES)
 
 Compute the constraints C for the DOCP problem (modeled as LB <= C(X) <= UB).
 """
-function DOCP_constraints_param!(c, xu, docp::DOCP)
-
-    # initialization
-    if docp.is_free_initial_time || docp.is_free_final_time
-        get_time_grid!(xu, docp)
-    end
-    v = get_OCP_variable_param(xu, docp)
-    # NB using inplace here did not seem to give better results
-    work = setWorkArray_param(docp, xu, docp.NLP_time_grid, v)
-
-    # main loop on time steps 
-    for i = 1:docp.dim_NLP_steps+1
-        setConstraintBlock_param!(docp, c, xu, v, docp.NLP_time_grid, i, work)
-    end
-
-    # point constraints
-    setPointConstraints_param!(docp, c, xu, v)
-
-    # NB. the function *needs* to return c for AD...
-    return c
-end
-
 function DOCP_constraints!(c, xu, docp::DOCP)
 
     # initialization
-    if docp.is_free_initial_time || docp.is_free_final_time
-        get_time_grid!(xu, docp)
-    end
-    v = get_optim_variable(xu, docp)
-    # NB using inplace here did not seem to give better results
-    work = setWorkArray(docp, xu, docp.NLP_time_grid, v)
+    time_grid = get_time_grid!(xu, docp)
+    v = get_OCP_variable(xu, docp)
+    work = setWorkArray(docp, xu, time_grid, v)
 
     # main loop on time steps 
     for i = 1:docp.dim_NLP_steps+1
-        setConstraintBlock!(docp, c, xu, v, docp.NLP_time_grid, i, work)
+        setConstraintBlock!(docp, c, xu, v, time_grid, i, work)
     end
 
     # point constraints
@@ -594,14 +437,14 @@ $(TYPEDSIGNATURES)
 
 Set the boundary and variable constraints
 """
-function setPointConstraints_param!(docp::DOCP, c, xu, v)
+function setPointConstraints!(docp::DOCP, c, xu, v)
 
     # offset
     offset = docp.dim_NLP_steps * (docp.dim_NLP_x * (1+docp.discretization.stage) + docp.dim_path_cons) + docp.dim_path_cons
 
     # variables
-    x0 = get_OCP_state_at_time_step_param(xu, docp, 1)
-    xf = get_OCP_state_at_time_step_param(xu, docp, docp.dim_NLP_steps+1)
+    x0 = get_OCP_state_at_time_step(xu, docp, 1)
+    xf = get_OCP_state_at_time_step(xu, docp, docp.dim_NLP_steps+1)
 
     # boundary constraints
     if docp.dim_boundary_cons > 0
@@ -627,38 +470,6 @@ function setPointConstraints_param!(docp::DOCP, c, xu, v)
     end
 end
 
-function setPointConstraints!(docp::DOCP, c, xu, v)
-
-    # offset
-    offset = docp.dim_NLP_steps * (docp.dim_NLP_x * (1+docp.discretization.stage) + docp.dim_path_cons) + docp.dim_path_cons
-
-    # variables
-    x0 = get_state_at_time_step(xu, docp, 1)
-    xf = get_state_at_time_step(xu, docp, docp.dim_NLP_steps+1)
-
-    # boundary constraints
-    if docp.dim_boundary_cons > 0
-        if docp.is_inplace
-            docp.boundary_constraints[2]((@view c[offset+1:offset+docp.dim_boundary_cons]), docp._x(x0), docp._x(xf), v)
-        else
-            c[offset+1:offset+docp.dim_boundary_cons] = docp.boundary_constraints[2](docp._x(x0), docp._x(xf), v)
-        end
-    end
-
-    # variable constraints
-    if docp.dim_v_cons > 0
-        if docp.is_inplace
-            docp.variable_constraints[2]((@view c[offset+docp.dim_boundary_cons+1:offset+docp.dim_boundary_cons+docp.dim_v_cons]), v)
-        else
-            c[offset+docp.dim_boundary_cons+1:offset+docp.dim_boundary_cons+docp.dim_v_cons] = docp.variable_constraints[2](v)
-        end
-    end
-
-    # null initial condition for lagrangian cost state
-    if docp.is_lagrange
-        c[offset+docp.dim_boundary_cons+docp.dim_v_cons+1] = x0[end]
-    end
-end
 
 """
 $(TYPEDSIGNATURES)
diff --git a/src/solution.jl b/src/solution.jl
index 3c55b82..f3de098 100644
--- a/src/solution.jl
+++ b/src/solution.jl
@@ -45,9 +45,8 @@ function CTBase.OptimalControlSolution(
     mult_UB = nothing,
 )
 
-    # time grid (NB. we get a Vector{Any} that requires explicit conversion)
-    get_time_grid!(primal, docp) 
-    T = convert(Vector{Float64},docp.NLP_time_grid)
+    # time grid
+    T = get_time_grid!(primal, docp)
     
     # recover primal variables
     X, U, v, box_multipliers =
@@ -98,7 +97,7 @@ function parse_DOCP_solution_primal(docp, solution; mult_LB = nothing, mult_UB =
 
     # state and control variables
     N = docp.dim_NLP_steps
-    X = zeros(N + 1, docp.dim_NLP_x)
+    X = zeros(N + 1, docp.dim_OCP_x)
     U = zeros(N + 1, docp.dim_NLP_u)
     v = Float64[]
 
@@ -109,8 +108,8 @@ function parse_DOCP_solution_primal(docp, solution; mult_LB = nothing, mult_UB =
     if isnothing(mult_UB) || length(mult_UB) == 0
         mult_UB = zeros(docp.dim_NLP_variables)
     end
-    mult_state_box_lower = zeros(N + 1, docp.dim_NLP_x)
-    mult_state_box_upper = zeros(N + 1, docp.dim_NLP_x)
+    mult_state_box_lower = zeros(N + 1, docp.dim_OCP_x)
+    mult_state_box_upper = zeros(N + 1, docp.dim_OCP_x)
     mult_control_box_lower = zeros(N + 1, docp.dim_NLP_u)
     mult_control_box_upper = zeros(N + 1, docp.dim_NLP_u)
     mult_variable_box_lower = zeros(N + 1, docp.dim_NLP_v)
@@ -118,23 +117,23 @@ function parse_DOCP_solution_primal(docp, solution; mult_LB = nothing, mult_UB =
 
     # retrieve optimization variables
     if docp.is_variable
-        v = get_optim_variable(solution, docp)
-        mult_variable_box_lower = get_optim_variable(mult_LB, docp)
-        mult_variable_box_upper = get_optim_variable(mult_UB, docp)
+        v = get_OCP_variable(solution, docp)
+        mult_variable_box_lower = get_OCP_variable(mult_LB, docp)
+        mult_variable_box_upper = get_OCP_variable(mult_UB, docp)
     end
 
     # loop over time steps
     disc = docp.discretization
     for i = 1:(N + 1)
         # state and control variables at current step
-        X[i,:] = get_state_at_time_step(solution, docp, i)
-        U[i,:] = get_control_at_time_step(solution, docp, i)
+        X[i,:] .= get_OCP_state_at_time_step(solution, docp, i)
+        U[i,:] .= get_OCP_control_at_time_step(solution, docp, i)
 
         # box multipliers
-        mult_state_box_lower[i, :] = get_state_at_time_step(mult_LB, docp, i)
-        mult_state_box_upper[i, :] = get_state_at_time_step(mult_UB, docp, i)
-        mult_control_box_lower[i, :] = get_control_at_time_step(mult_LB, docp, i)
-        mult_control_box_upper[i, :] = get_control_at_time_step(mult_UB, docp, i)
+        mult_state_box_lower[i, :] .= get_OCP_state_at_time_step(mult_LB, docp, i)
+        mult_state_box_upper[i, :] .= get_OCP_state_at_time_step(mult_UB, docp, i)
+        mult_control_box_lower[i, :] .= get_OCP_control_at_time_step(mult_LB, docp, i)
+        mult_control_box_upper[i, :] .= get_OCP_control_at_time_step(mult_UB, docp, i)
 
     end
 
diff --git a/src/trapeze.jl b/src/trapeze.jl
index 51fd6f5..f9d7e2b 100644
--- a/src/trapeze.jl
+++ b/src/trapeze.jl
@@ -17,10 +17,10 @@ struct Trapeze <: Discretization
 end
 
 # not only for Trapeze, but may be redefined if needed
-function get_OCP_variable_param(xu, docp::DOCP{<: Discretization, <: ScalVect, <: ScalVect, ScalVariable})
+function get_OCP_variable(xu, docp::DOCP{<: Discretization, <: ScalVect, <: ScalVect, ScalVariable})
     return xu[docp.dim_NLP_variables]
 end
-function get_OCP_variable_param(xu, docp::DOCP{<: Discretization, <: ScalVect, <: ScalVect, VectVariable})
+function get_OCP_variable(xu, docp::DOCP{<: Discretization, <: ScalVect, <: ScalVect, VectVariable})
     return @view xu[(docp.dim_NLP_variables - docp.dim_NLP_v + 1):docp.dim_NLP_variables]
 end
 
@@ -31,25 +31,11 @@ $(TYPEDSIGNATURES)
 Retrieve state and control variables at given time step from the NLP variables.
 Convention: 1 <= i <= dim_NLP_steps+1
 """
-function get_state_at_time_step(xu, docp::DOCP{Trapeze}, i)
-    offset = (i-1) * (docp.dim_NLP_x + docp.dim_NLP_u)
-    return @view xu[(offset + 1):(offset + docp.dim_NLP_x)]
-end
-   
-function get_OCP_state_at_time_step(xu, docp::DOCP{Trapeze}, i)
-    offset = (i-1) * (docp.dim_NLP_x + docp.dim_NLP_u)
-    if docp.dim_OCP_x == 1
-        return xu[offset+1]
-    else
-        return @view xu[(offset + 1):(offset + docp.dim_OCP_x)]
-    end
-end
-
-function get_OCP_state_at_time_step_param(xu, docp::DOCP{Trapeze, ScalVariable, <: ScalVect, <: ScalVect}, i)
+function get_OCP_state_at_time_step(xu, docp::DOCP{Trapeze, ScalVariable, <: ScalVect, <: ScalVect}, i)
     offset = (i-1) * (docp.dim_NLP_x + docp.dim_NLP_u)
     return xu[offset+1]
 end
-function get_OCP_state_at_time_step_param(xu, docp::DOCP{Trapeze, VectVariable, <: ScalVect, <: ScalVect}, i)
+function get_OCP_state_at_time_step(xu, docp::DOCP{Trapeze, VectVariable, <: ScalVect, <: ScalVect}, i)
     offset = (i-1) * (docp.dim_NLP_x + docp.dim_NLP_u)
     return @view xu[(offset + 1):(offset + docp.dim_OCP_x)]
 end
@@ -59,25 +45,12 @@ function get_lagrange_state_at_time_step(xu, docp, i)
     return xu[offset + docp.dim_NLP_x]
 end
 
-function get_control_at_time_step(xu, docp::DOCP{Trapeze}, i)
-    offset = (i-1) * (docp.dim_NLP_x + docp.dim_NLP_u) + docp.dim_NLP_x
-    return @view xu[(offset + 1):(offset + docp.dim_NLP_u)]
-end
-
-function get_OCP_control_at_time_step(xu, docp::DOCP{Trapeze}, i)
-    offset = (i-1) * (docp.dim_NLP_x + docp.dim_NLP_u) + docp.dim_NLP_x
-    if docp.dim_NLP_u == 1
-        return xu[offset+1]
-    else
-        return @view xu[(offset + 1):(offset + docp.dim_NLP_u)]
-    end
-end
 
-function get_OCP_control_at_time_step_param(xu, docp::DOCP{Trapeze, <: ScalVect, ScalVariable, <: ScalVect}, i)
+function get_OCP_control_at_time_step(xu, docp::DOCP{Trapeze, <: ScalVect, ScalVariable, <: ScalVect}, i)
     offset = (i-1) * (docp.dim_NLP_x + docp.dim_NLP_u) + docp.dim_NLP_x
     return xu[offset+1]
 end
-function get_OCP_control_at_time_step_param(xu, docp::DOCP{Trapeze, <: ScalVect, VectVariable, <: ScalVect}, i)
+function get_OCP_control_at_time_step(xu, docp::DOCP{Trapeze, <: ScalVect, VectVariable, <: ScalVect}, i)
     offset = (i-1) * (docp.dim_NLP_x + docp.dim_NLP_u) + docp.dim_NLP_x
     return @view xu[(offset + 1):(offset + docp.dim_NLP_u)]
 end
@@ -104,11 +77,11 @@ function set_control_at_time_step!(xu, u_init, docp::DOCP{Trapeze}, i)
 end
 
 
-function setWorkArray_param(docp::DOCP{Trapeze}, xu, time_grid, v)
+function setWorkArray(docp::DOCP{Trapeze}, xu, time_grid, v)
     work = similar(xu, docp.dim_NLP_x)
     t0 = time_grid[1]
-    x0 = get_OCP_state_at_time_step_param(xu, docp, 1)
-    u0 = get_OCP_control_at_time_step_param(xu, docp, 1)
+    x0 = get_OCP_state_at_time_step(xu, docp, 1)
+    u0 = get_OCP_control_at_time_step(xu, docp, 1)
     if docp.is_inplace
         docp.ocp.dynamics((@view work[1:docp.dim_OCP_x]), t0, x0, u0, v)
     else
@@ -123,14 +96,7 @@ function setWorkArray_param(docp::DOCP{Trapeze}, xu, time_grid, v)
     end
     return work
 end
-function setWorkArray(docp::DOCP{Trapeze}, xu, time_grid, v)
-    work = similar(xu, docp.dim_NLP_x)
-    t0 = time_grid[1]
-    x0 = get_state_at_time_step(xu, docp, 1)
-    u0 = get_control_at_time_step(xu, docp, 1)
-    docp.dynamics_ext!(work, t0, x0, u0, v)
-    return work
-end
+
 
 """
 $(TYPEDSIGNATURES)
@@ -138,23 +104,23 @@ $(TYPEDSIGNATURES)
 Set the constraints corresponding to the state equation
 Convention: 1 <= i <= dim_NLP_steps (+1)
 """
-function setConstraintBlock_param!(docp::DOCP{Trapeze}, c, xu, v, time_grid, i, work)
+function setConstraintBlock!(docp::DOCP{Trapeze}, c, xu, v, time_grid, i, work)
 
     # offset for previous steps
     offset = (i-1)*(docp.dim_NLP_x + docp.dim_path_cons)
 
     # 0. variables
     ti = time_grid[i]
-    xi = get_OCP_state_at_time_step_param(xu, docp, i)
-    ui = get_OCP_control_at_time_step_param(xu, docp, i)
+    xi = get_OCP_state_at_time_step(xu, docp, i)
+    ui = get_OCP_control_at_time_step(xu, docp, i)
 
     #1. state equation
     if i <= docp.dim_NLP_steps
         # more variables
         fi = copy(work) # create new copy, not just a reference
         tip1 = time_grid[i+1]
-        xip1 = get_OCP_state_at_time_step_param(xu, docp, i+1)
-        uip1 = get_OCP_control_at_time_step_param(xu, docp, i+1)
+        xip1 = get_OCP_state_at_time_step(xu, docp, i+1)
+        uip1 = get_OCP_control_at_time_step(xu, docp, i+1)
         if docp.is_inplace
             docp.ocp.dynamics((@view work[1:docp.dim_OCP_x]), tip1, xip1, uip1, v)
         else
@@ -171,9 +137,9 @@ function setConstraintBlock_param!(docp::DOCP{Trapeze}, c, xu, v, time_grid, i,
         # trapeze rule with 'smart' update for dynamics (@. allocs more and is slower -_-)
         # or split dyn and lag ?
         if docp.dim_OCP_x == 1
-            c[offset+1] = xip1 - (xi + 0.5 * (tip1 - ti) * (fi + work))
+            c[offset+1] = xip1 - (xi + 0.5 * (tip1 - ti) * (fi[1] + work[1]))
         else
-            c[offset+1:offset+docp.dim_OCP_x] = xip1 - (xi + 0.5 * (tip1 - ti) * (fi + work))
+            c[offset+1:offset+docp.dim_OCP_x] = xip1 - (xi + 0.5 * (tip1 - ti) * (fi[1:docp.dim_OCP_x] + work[1:docp.dim_OCP_x]))
         end
         if docp.is_lagrange
             c[offset+docp.dim_NLP_x] = get_lagrange_state_at_time_step(xu, docp, i+1) - (get_lagrange_state_at_time_step(xu, docp, i) + 0.5 * (tip1 - ti) * (fi[docp.dim_NLP_x] + work[docp.dim_NLP_x]))
@@ -206,56 +172,3 @@ function setConstraintBlock_param!(docp::DOCP{Trapeze}, c, xu, v, time_grid, i,
     end
 
 end
-
-
-function setConstraintBlock!(docp::DOCP{Trapeze}, c, xu, v, time_grid, i, work)
-
-    disc = docp.discretization
-
-    # offset for previous steps
-    offset = (i-1)*(docp.dim_NLP_x + docp.dim_path_cons)
-
-    # 0. variables
-    ti = time_grid[i]
-    xi = get_state_at_time_step(xu, docp, i)
-    ui = get_control_at_time_step(xu, docp, i)
-
-    #1. state equation
-    if i <= docp.dim_NLP_steps
-        # more variables
-        fi = copy(work) # create new copy, not just a reference
-        tip1 = time_grid[i+1]
-        xip1 = get_state_at_time_step(xu, docp, i+1)
-        uip1 = get_control_at_time_step(xu, docp, i+1)
-        docp.dynamics_ext!(work, tip1, xip1, uip1, v)
-
-        # trapeze rule with 'smart' update for dynamics (similar with @.)
-        c[offset+1:offset+docp.dim_NLP_x] = xip1 - (xi + 0.5 * (tip1 - ti) * (fi + work)) #+++ jet runtime dispatch here even with explicit index ranges
-        offset += docp.dim_NLP_x
-    end
-
-    # 2. path constraints
-    # Notes on allocations:.= seems similar
-    if docp.dim_u_cons > 0
-        if docp.is_inplace
-            docp.control_constraints[2]((@view c[offset+1:offset+docp.dim_u_cons]),ti, docp._u(ui), v)
-        else
-            c[offset+1:offset+docp.dim_u_cons] = docp.control_constraints[2](ti, docp._u(ui), v)
-        end
-    end
-    if docp.dim_x_cons > 0 
-        if docp.is_inplace
-            docp.state_constraints[2]((@view c[offset+docp.dim_u_cons+1:offset+docp.dim_u_cons+docp.dim_x_cons]),ti, docp._x(xi), v)
-        else
-            c[offset+docp.dim_u_cons+1:offset+docp.dim_u_cons+docp.dim_x_cons] = docp.state_constraints[2](ti, docp._x(xi), v)
-        end
-    end
-    if docp.dim_mixed_cons > 0 
-        if docp.is_inplace
-            docp.mixed_constraints[2]((@view c[offset+docp.dim_u_cons+docp.dim_x_cons+1:offset+docp.dim_u_cons+docp.dim_x_cons+docp.dim_mixed_cons]), ti, docp._x(xi), docp._u(ui), v)
-        else
-            c[offset+docp.dim_u_cons+docp.dim_x_cons+1:offset+docp.dim_u_cons+docp.dim_x_cons+docp.dim_mixed_cons] = docp.mixed_constraints[2](ti, docp._x(xi), docp._u(ui), v)
-        end
-    end
-
-end
diff --git a/src/utils.jl b/src/utils.jl
index 2b73f93..c13e380 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -5,32 +5,12 @@ Internal layout for NLP variables:
 with the convention u([t_i,t_i+1[) = U_i and u(tf) = U_N-1
 =#
 
-# getter for optimization variables
-function get_optim_variable(xu, docp)
-    # empty for dim 0
-    return @view xu[(end - docp.dim_NLP_v + 1):end]
-end
-
-function get_OCP_variable(xu, docp)
-    # empty for dim 0
-    if docp.dim_NLP_v == 1
-        return xu[end]
-    else
-        return @view xu[(end - docp.dim_NLP_v + 1):end]
-    end
-end
-#=function get_OCP_variable2(xu, docp, dim_v::Val{1})
-    return xu[end]
-end
-function get_OCP_variable2(xu, docp, dim_v)
-    return @view xu[(end - docp.dim_NLP_v + 1):end]
-end=#
 
 # +++ in problem.jl ?
 # getters for initial and final time
 function get_initial_time(xu, docp)
     if docp.is_free_initial_time
-        return get_optim_variable(xu, docp)[docp.index_initial_time]
+        return get_OCP_variable(xu, docp)[docp.index_initial_time]
     else
         return docp.fixed_initial_time
     end
@@ -38,7 +18,7 @@ end
 
 function get_final_time(xu, docp)
     if docp.is_free_final_time
-        return get_optim_variable(xu, docp)[docp.index_final_time]
+        return get_OCP_variable(xu, docp)[docp.index_final_time]
     else
         return docp.fixed_final_time
     end
@@ -46,10 +26,13 @@ end
 
 # time grid
 function get_time_grid!(xu, docp)
-    t0 = get_initial_time(xu, docp)
-    tf = get_final_time(xu, docp)
-    @. docp.NLP_time_grid = t0 + docp.NLP_normalized_time_grid * (tf - t0)
-    return
+    if docp.is_free_initial_time || docp.is_free_final_time
+        t0 = get_initial_time(xu, docp)
+        tf = get_final_time(xu, docp)
+        return @. t0 + docp.NLP_normalized_time_grid * (tf - t0)
+    else
+        return docp.NLP_time_grid
+    end
 end