diff --git a/src/CTBase.jl b/src/CTBase.jl
index 58488d3e..bdccaf49 100644
--- a/src/CTBase.jl
+++ b/src/CTBase.jl
@@ -11,7 +11,6 @@ $(EXPORTS)
 """
 module CTBase
 
-# using
 import Base
 using DocStringExtensions
 using DifferentiationInterface:
@@ -23,18 +22,18 @@ using DifferentiationInterface:
     prepare_gradient,
     prepare_jacobian
 import ForwardDiff
-using Interpolations: linear_interpolation, Line, Interpolations # for default interpolation
-using MLStyle # pattern matching
+using Interpolations: linear_interpolation, Line, Interpolations # For default interpolation
+using MLStyle # Pattern matching
 using Parameters # @with_kw: to have default values in struct
-using Printf # to print an Opt imalControlModel
+using Printf # To print an OptimalControlModel
 using DataStructures # OrderedDict for aliases
-using Unicode # unicode primitives
-using PrettyTables # to print a table
+using Unicode # Unicode primitives
+using PrettyTables # To print a table
 using ReplMaker
 using MacroTools: @capture, postwalk, striplines
 using LinearAlgebra
 
-# to suppress ambiguities 
+# To suppress ambiguities 
 using SparseArrays, StaticArrays
 
 # --------------------------------------------------------------------------------------------------
@@ -229,7 +228,7 @@ include("repl.jl")
 include("init.jl")
 include("print.jl")
 
-# numeric types
+# Numeric types
 export ctNumber, ctVector, Time, Times, TimesDisc
 
 export States, Costates, Controls, State, Costate, Control, Variable, Dimension, Index
@@ -237,10 +236,10 @@ export DState, DCostate
 export TimeDependence, Autonomous, NonAutonomous
 export VariableDependence, NonFixed, Fixed
 
-# description
+# Description
 export Description, add, getFullDescription, remove
 
-# exceptions
+# Exceptions
 export CTException, ParsingError, AmbiguousDescription, IncorrectMethod
 export IncorrectArgument, IncorrectOutput, NotImplemented, UnauthorizedCall
 export ExtensionError
@@ -248,26 +247,24 @@ export ExtensionError
 # AD
 export set_AD_backend
 
-# functions
+# Functions
 export Hamiltonian, HamiltonianVectorField, VectorField
 export Mayer, Lagrange, Dynamics, ControlLaw, FeedbackControl, Multiplier
+export Mayer!, Lagrange!, Dynamics!
 export BoundaryConstraint, StateConstraint, ControlConstraint, MixedConstraint, VariableConstraint
+export BoundaryConstraint!, StateConstraint!, ControlConstraint!, MixedConstraint!, VariableConstraint!
 
-# model
+# Model
 export OptimalControlModel
 export Model
-export __OCPModel # redirection to Model to avoid confusion with other Model functions from other packages. Due to @def macro
+export __OCPModel # todo: to be updated, redirection to Model to avoid confusion with other Model functions from other packages. Due to @def macro
 export variable!,
-    time!,
-    constraint!,
-    dynamics!,
-    objective!,
-    state!,
-    control!,
-    remove_constraint!,
-    model_expression!
-export is_autonomous, is_fixed, is_time_independent, is_time_dependent, is_min, is_max
+    time!, constraint!, dynamics!, objective!, state!, control!, remove_constraint!, model_expression!
+export is_autonomous, is_fixed
+export is_time_independent, is_time_dependent
+export is_min, is_max
 export is_variable_dependent, is_variable_independent
+export is_in_place
 export nlp_constraints!, constraints, constraints_labels, constraint
 export has_free_final_time, has_free_initial_time, has_lagrange_cost, has_mayer_cost
 export dim_control_constraints, dim_state_constraints, dim_mixed_constraints, dim_path_constraints
@@ -278,8 +275,9 @@ export control_dimension, control_components_names, control_name
 export state_dimension, state_components_names, state_name
 export variable_dimension, variable_components_names, variable_name
 export lagrange, mayer, criterion, dynamics
+export __constraint, __lagrange, __mayer, __dynamics # todo: remove after in place tests
 
-# solution
+# Solution
 export OptimalControlSolution
 export time_grid, control, state, variable, costate, objective
 export state_discretized, control_discretized, costate_discretized
@@ -301,13 +299,13 @@ export control_constraints!,
 export state_constraints!, mult_state_constraints!, mult_state_box_lower!, mult_state_box_upper!
 export mixed_constraints!, mult_mixed_constraints!
 
-# initialization
+# Initialization
 export OptimalControlInit
 
-# utils
+# Utils
 export ctgradient, ctjacobian, ctinterpolate, ctindices, ctupperscripts
 
-# differential geometry
+# Differential geometry
 export Lie, @Lie, Poisson, HamiltonianLift, AbstractHamiltonian, Lift, ⋅, ∂ₜ
 
 # ctparser_utils
@@ -320,4 +318,4 @@ export @def
 export ct_repl, ct_repl_update_model
 isdefined(Base, :active_repl) && ct_repl()
 
-end
+end
\ No newline at end of file
diff --git a/src/functions.jl b/src/functions.jl
index 21fe3c23..0a5bb4fc 100644
--- a/src/functions.jl
+++ b/src/functions.jl
@@ -17,6 +17,11 @@ function BoundaryConstraint(f::Function; variable::Bool = false)
     return BoundaryConstraint{variable_dependence}(f)
 end
 
+function BoundaryConstraint!(f!::Function; variable::Bool = false)
+    variable_dependence = variable ? NonFixed : Fixed
+    return BoundaryConstraint!{variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -36,6 +41,12 @@ function BoundaryConstraint(f::Function, dependencies::DataType...)
     return BoundaryConstraint{variable_dependence}(f)
 end
 
+function BoundaryConstraint!(f!::Function, dependencies::DataType...)
+    __check_dependencies(dependencies)
+    variable_dependence = NonFixed ∈ dependencies ? NonFixed : Fixed
+    return BoundaryConstraint!{variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -66,6 +77,18 @@ function (F::BoundaryConstraint{NonFixed})(x0::State, xf::State, v::Variable)::c
     return F.f(x0, xf, v)
 end
 
+function (F::BoundaryConstraint!{Fixed})(r::ctVector, x0::State, xf::State)::Nothing
+    return F.f!(r, x0, xf)
+end
+
+function (F::BoundaryConstraint!{Fixed})(r::ctVector, x0::State, xf::State, v::Variable)::Nothing
+    return F.f!(r, x0, xf)
+end
+
+function (F::BoundaryConstraint!{NonFixed})(r::ctVector, x0::State, xf::State, v::Variable)::Nothing
+    return F.f!(r, x0, xf, v)
+end
+
 # --------------------------------------------------------------------------------------------------
 """
 
@@ -85,6 +108,11 @@ function Mayer(f::Function; variable::Bool = false)
     return Mayer{variable_dependence}(f)
 end
 
+function Mayer!(f!::Function; variable::Bool = false)
+    variable_dependence = variable ? NonFixed : Fixed
+    return Mayer!{variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -104,6 +132,12 @@ function Mayer(f::Function, dependencies::DataType...)
     return Mayer{variable_dependence}(f)
 end
 
+function Mayer!(f!::Function, dependencies::DataType...)
+    __check_dependencies(dependencies)
+    variable_dependence = NonFixed ∈ dependencies ? NonFixed : Fixed
+    return Mayer!{variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -134,6 +168,18 @@ function (F::Mayer{NonFixed})(x0::State, xf::State, v::Variable)::ctNumber
     return F.f(x0, xf, v)
 end
 
+function (F::Mayer!{Fixed})(r::ctVector, x0::State, xf::State)::Nothing
+    return F.f!(r, x0, xf)
+end
+
+function (F::Mayer!{Fixed})(r::ctVector, x0::State, xf::State, v::Variable)::Nothing
+    return F.f!(r, x0, xf)
+end
+
+function (F::Mayer!{NonFixed})(r::ctVector, x0::State, xf::State, v::Variable)::Nothing
+    return F.f!(r, x0, xf, v)
+end
+
 # --------------------------------------------------------------------------------------------------
 # Hamiltonian
 """
@@ -672,6 +718,12 @@ function Lagrange(f::Function; autonomous::Bool = true, variable::Bool = false)
     return Lagrange{time_dependence, variable_dependence}(f)
 end
 
+function Lagrange!(f!::Function; autonomous::Bool = true, variable::Bool = false)
+    time_dependence = autonomous ? Autonomous : NonAutonomous
+    variable_dependence = variable ? NonFixed : Fixed
+    return Lagrange!{time_dependence, variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -700,6 +752,13 @@ function Lagrange(f::Function, dependencies::DataType...)
     return Lagrange{time_dependence, variable_dependence}(f)
 end
 
+function Lagrange!(f!::Function, dependencies::DataType...)
+    __check_dependencies(dependencies)
+    time_dependence = NonAutonomous ∈ dependencies ? NonAutonomous : Autonomous
+    variable_dependence = NonFixed ∈ dependencies ? NonFixed : Fixed
+    return Lagrange!{time_dependence, variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -775,6 +834,40 @@ function (F::Lagrange{NonAutonomous, NonFixed})(
     return F.f(t, x, u, v)
 end
 
+function (F::Lagrange!{Autonomous, Fixed})(r::ctVector, x::State, u::Control)::Nothing
+    return F.f!(r, x, u)
+end
+
+function (F::Lagrange!{Autonomous, Fixed})(r::ctVector, t::Time, x::State, u::Control, v::Variable)::Nothing
+    return F.f!(r, x, u)
+end
+
+function (F::Lagrange!{Autonomous, NonFixed})(r::ctVector, x::State, u::Control, v::Variable)::Nothing
+    return F.f!(r, x, u, v)
+end
+
+function (F::Lagrange!{Autonomous, NonFixed})(r::ctVector, t::Time, x::State, u::Control, v::Variable)::Nothing
+    return F.f!(r, x, u, v)
+end
+
+function (F::Lagrange!{NonAutonomous, Fixed})(r::ctVector, t::Time, x::State, u::Control)::Nothing
+    return F.f!(r, t, x, u)
+end
+
+function (F::Lagrange!{NonAutonomous, Fixed})(r::ctVector, t::Time, x::State, u::Control, v::Variable)::Nothing
+    return F.f!(r, t, x, u)
+end
+
+function (F::Lagrange!{NonAutonomous, NonFixed})(
+    r::ctVector,
+    t::Time,
+    x::State,
+    u::Control,
+    v::Variable,
+)::Nothing
+    return F.f!(r, t, x, u, v)
+end
+
 # --------------------------------------------------------------------------------------------------
 """
 
@@ -796,6 +889,12 @@ function Dynamics(f::Function; autonomous::Bool = true, variable::Bool = false)
     return Dynamics{time_dependence, variable_dependence}(f)
 end
 
+function Dynamics!(f!::Function; autonomous::Bool = true, variable::Bool = false)
+    time_dependence = autonomous ? Autonomous : NonAutonomous
+    variable_dependence = variable ? NonFixed : Fixed
+    return Dynamics!{time_dependence, variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -821,6 +920,13 @@ function Dynamics(f::Function, dependencies::DataType...)
     return Dynamics{time_dependence, variable_dependence}(f)
 end
 
+function Dynamics!(f!::Function, dependencies::DataType...)
+    __check_dependencies(dependencies)
+    time_dependence = NonAutonomous ∈ dependencies ? NonAutonomous : Autonomous
+    variable_dependence = NonFixed ∈ dependencies ? NonFixed : Fixed
+    return Dynamics!{time_dependence, variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -883,6 +989,40 @@ function (F::Dynamics{NonAutonomous, NonFixed})(
     return F.f(t, x, u, v)
 end
 
+function (F::Dynamics!{Autonomous, Fixed})(r::ctVector, x::State, u::Control)::Nothing
+    return F.f!(r, x, u)
+end
+
+function (F::Dynamics!{Autonomous, Fixed})(r::ctVector, t::Time, x::State, u::Control, v::Variable)::Nothing
+    return F.f!(r, x, u)
+end
+
+function (F::Dynamics!{Autonomous, NonFixed})(r::ctVector, x::State, u::Control, v::Variable)::Nothing
+    return F.f!(r, x, u, v)
+end
+
+function (F::Dynamics!{Autonomous, NonFixed})(r::ctVector, t::Time, x::State, u::Control, v::Variable)::Nothing
+    return F.f!(r, x, u, v)
+end
+
+function (F::Dynamics!{NonAutonomous, Fixed})(r::ctVector, t::Time, x::State, u::Control)::Nothing
+    return F.f!(r, t, x, u)
+end
+
+function (F::Dynamics!{NonAutonomous, Fixed})(r::ctVector, t::Time, x::State, u::Control, v::Variable)::Nothing
+    return F.f!(r, t, x, u)
+end
+
+function (F::Dynamics!{NonAutonomous, NonFixed})(
+    r::ctVector,
+    t::Time,
+    x::State,
+    u::Control,
+    v::Variable,
+)::Nothing
+    return F.f!(r, t, x, u, v)
+end
+
 # --------------------------------------------------------------------------------------------------
 """
 
@@ -904,6 +1044,12 @@ function StateConstraint(f::Function; autonomous::Bool = true, variable::Bool =
     return StateConstraint{time_dependence, variable_dependence}(f)
 end
 
+function StateConstraint!(f!::Function; autonomous::Bool = true, variable::Bool = false)
+    time_dependence = autonomous ? Autonomous : NonAutonomous
+    variable_dependence = variable ? NonFixed : Fixed
+    return StateConstraint!{time_dependence, variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -929,6 +1075,13 @@ function StateConstraint(f::Function, dependencies::DataType...)
     return StateConstraint{time_dependence, variable_dependence}(f)
 end
 
+function StateConstraint!(f!::Function, dependencies::DataType...)
+    __check_dependencies(dependencies)
+    time_dependence = NonAutonomous ∈ dependencies ? NonAutonomous : Autonomous
+    variable_dependence = NonFixed ∈ dependencies ? NonFixed : Fixed
+    return StateConstraint!{time_dependence, variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -990,6 +1143,34 @@ function (F::StateConstraint{NonAutonomous, NonFixed})(t::Time, x::State, v::Var
     return F.f(t, x, v)
 end
 
+function (F::StateConstraint!{Autonomous, Fixed})(r::ctVector, x::State)::Nothing
+    return F.f!(r, x)
+end
+
+function (F::StateConstraint!{Autonomous, Fixed})(r::ctVector, t::Time, x::State, v::Variable)::Nothing
+    return F.f!(r, x)
+end
+
+function (F::StateConstraint!{Autonomous, NonFixed})(r::ctVector, x::State, v::Variable)::Nothing
+    return F.f!(r, x, v)
+end
+
+function (F::StateConstraint!{Autonomous, NonFixed})(r::ctVector, t::Time, x::State, v::Variable)::Nothing
+    return F.f!(r, x, v)
+end
+
+function (F::StateConstraint!{NonAutonomous, Fixed})(r::ctVector, t::Time, x::State)::Nothing
+    return F.f!(r, t, x)
+end
+
+function (F::StateConstraint!{NonAutonomous, Fixed})(r::ctVector, t::Time, x::State, v::Variable)::Nothing
+    return F.f!(r, t, x)
+end
+
+function (F::StateConstraint!{NonAutonomous, NonFixed})(r::ctVector, t::Time, x::State, v::Variable)::Nothing
+    return F.f!(r, t, x, v)
+end
+
 # --------------------------------------------------------------------------------------------------
 
 """
@@ -1012,6 +1193,12 @@ function ControlConstraint(f::Function; autonomous::Bool = true, variable::Bool
     return ControlConstraint{time_dependence, variable_dependence}(f)
 end
 
+function ControlConstraint!(f!::Function; autonomous::Bool = true, variable::Bool = false)
+    time_dependence = autonomous ? Autonomous : NonAutonomous
+    variable_dependence = variable ? NonFixed : Fixed
+    return ControlConstraint!{time_dependence, variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -1035,6 +1222,13 @@ function ControlConstraint(f::Function, dependencies::DataType...)
     return ControlConstraint{time_dependence, variable_dependence}(f)
 end
 
+function ControlConstraint!(f!::Function, dependencies::DataType...) # todo: useless? only use version above? 
+    __check_dependencies(dependencies)
+    time_dependence = NonAutonomous ∈ dependencies ? NonAutonomous : Autonomous
+    variable_dependence = NonFixed ∈ dependencies ? NonFixed : Fixed
+    return ControlConstraint!{time_dependence, variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -1094,6 +1288,34 @@ function (F::ControlConstraint{NonAutonomous, NonFixed})(t::Time, u::Control, v:
     return F.f(t, u, v)
 end
 
+function (F::ControlConstraint!{Autonomous, Fixed})(r::ctVector, u::Control)::Nothing
+    return F.f!(r, u)
+end
+
+function (F::ControlConstraint!{Autonomous, Fixed})(r::ctVector, t::Time, u::Control, v::Variable)::Nothing
+    return F.f!(r, u)
+end
+
+function (F::ControlConstraint!{Autonomous, NonFixed})(r::ctVector, u::Control, v::Variable)::Nothing
+    return F.f!(r, u, v)
+end
+
+function (F::ControlConstraint!{Autonomous, NonFixed})(r::ctVector, t::Time, u::Control, v::Variable)::Nothing
+    return F.f!(r, u, v)
+end
+
+function (F::ControlConstraint!{NonAutonomous, Fixed})(r::ctVector, t::Time, u::Control)::Nothing
+    return F.f!(r, t, u)
+end
+
+function (F::ControlConstraint!{NonAutonomous, Fixed})(r::ctVector, t::Time, u::Control, v::Variable)::Nothing
+    return F.f!(r, t, u)
+end
+
+function (F::ControlConstraint!{NonAutonomous, NonFixed})(r::ctVector, t::Time, u::Control, v::Variable)::Nothing
+    return F.f!(r, t, u, v)
+end
+
 # --------------------------------------------------------------------------------------------------
 """
 
@@ -1115,6 +1337,12 @@ function MixedConstraint(f::Function; autonomous::Bool = true, variable::Bool =
     return MixedConstraint{time_dependence, variable_dependence}(f)
 end
 
+function MixedConstraint!(f!::Function; autonomous::Bool = true, variable::Bool = false)
+    time_dependence = autonomous ? Autonomous : NonAutonomous
+    variable_dependence = variable ? NonFixed : Fixed
+    return MixedConstraint!{time_dependence, variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -1140,6 +1368,13 @@ function MixedConstraint(f::Function, dependencies::DataType...)
     return MixedConstraint{time_dependence, variable_dependence}(f)
 end
 
+function MixedConstraint!(f!::Function, dependencies::DataType...)
+    __check_dependencies(dependencies)
+    time_dependence = NonAutonomous ∈ dependencies ? NonAutonomous : Autonomous
+    variable_dependence = NonFixed ∈ dependencies ? NonFixed : Fixed
+    return MixedConstraint!{time_dependence, variable_dependence}(f!)
+end
+
 """
 
 $(TYPEDSIGNATURES)
@@ -1223,6 +1458,58 @@ function (F::MixedConstraint{NonAutonomous, NonFixed})(
     return F.f(t, x, u, v)
 end
 
+function (F::MixedConstraint!{Autonomous, Fixed})(r::ctVector, x::State, u::Control)::Nothing
+    return F.f!(r, x, u)
+end
+
+function (F::MixedConstraint!{Autonomous, Fixed})(
+    r::ctVector,
+    t::Time,
+    x::State,
+    u::Control,
+    v::Variable,
+)::Nothing
+    return F.f!(r, x, u)
+end
+
+function (F::MixedConstraint!{Autonomous, NonFixed})(r::ctVector, x::State, u::Control, v::Variable)::Nothing
+    return F.f!(r, x, u, v)
+end
+
+function (F::MixedConstraint!{Autonomous, NonFixed})(
+    r::ctVector,
+    t::Time,
+    x::State,
+    u::Control,
+    v::Variable,
+)::Nothing
+    return F.f!(r, x, u, v)
+end
+
+function (F::MixedConstraint!{NonAutonomous, Fixed})(r::ctVector, t::Time, x::State, u::Control)::Nothing
+    return F.f!(r, t, x, u)
+end
+
+function (F::MixedConstraint!{NonAutonomous, Fixed})(
+    r::ctVector,
+    t::Time,
+    x::State,
+    u::Control,
+    v::Variable,
+)::Nothing
+    return F.f!(r, t, x, u)
+end
+
+function (F::MixedConstraint!{NonAutonomous, NonFixed})(
+    r::ctVector,
+    t::Time,
+    x::State,
+    u::Control,
+    v::Variable,
+)::Nothing
+    return F.f!(r, t, x, u, v)
+end
+
 # --------------------------------------------------------------------------------------------------
 """
 
@@ -1240,6 +1527,10 @@ function (F::VariableConstraint)(v::Variable)::ctVector
     return F.f(v)
 end
 
+function (F::VariableConstraint!)(r::ctVector, v::Variable)::Nothing
+    return F.f!(r, v)
+end
+
 # --------------------------------------------------------------------------------------------------
 """
 
diff --git a/src/onepass.jl b/src/onepass.jl
index f840fd1e..74f4325b 100644
--- a/src/onepass.jl
+++ b/src/onepass.jl
@@ -210,7 +210,7 @@ p_variable!(p, ocp, v, q; components_names = nothing, log = false) = begin
         p.aliases[Symbol(v, ctupperscripts(i))] = :($v^$i)
     end # make: v¹, v²... if the variable is named v
     if (isnothing(components_names))
-        __wrap(:(variable!($ocp, $q, $vv)), p.lnum, p.line)
+        code = :( variable!($ocp, $q, $vv) )
     else
         qq == length(components_names.args) ||
             return __throw("the number of variable components must be $qq", p.lnum, p.line)
@@ -218,8 +218,9 @@ p_variable!(p, ocp, v, q; components_names = nothing, log = false) = begin
             p.aliases[components_names.args[i]] = :($v[$i])
         end # aliases from names given by the user
         ss = QuoteNode(string.(components_names.args))
-        __wrap(:(variable!($ocp, $q, $vv, $ss)), p.lnum, p.line)
+        code = :( variable!($ocp, $q, $vv, $ss) )
     end
+    return __wrap(code, p.lnum, p.line) 
 end
 
 p_alias!(p, ocp, a, e; log = false) = begin
@@ -231,7 +232,8 @@ p_alias!(p, ocp, a, e; log = false) = begin
         p.aliases[Symbol(a, ctupperscripts(i))] = :($a^$i)
     end
     p.aliases[a] = e
-    __wrap(:(LineNumberNode(0, "alias: " * string($aa) * " = " * string($ee))), p.lnum, p.line)
+    code = :( LineNumberNode(0, "alias: " * string($aa) * " = " * string($ee)) )
+    return __wrap(code, p.lnum, p.line)
 end
 
 p_time!(p, ocp, t, t0, tf; log = false) = begin
@@ -273,7 +275,7 @@ p_time!(p, ocp, t, t0, tf; log = false) = begin
             _ => return __throw("bad time declaration", p.lnum, p.line)
         end
     end
-    __wrap(code, p.lnum, p.line)
+    return __wrap(code, p.lnum, p.line)
 end
 
 p_state!(p, ocp, x, n; components_names = nothing, log = false) = begin
@@ -284,25 +286,27 @@ p_state!(p, ocp, x, n; components_names = nothing, log = false) = begin
     nn = n isa Integer ? n : 9
     for i ∈ 1:nn
         p.aliases[Symbol(x, ctindices(i))] = :($x[$i])
-    end # make: x₁, x₂... if the state is named x
+    end # Make x₁, x₂... if the state is named x
     for i ∈ 1:nn
         p.aliases[Symbol(x, i)] = :($x[$i])
-    end # make: x1, x2... if the state is named x
+    end # Make x1, x2... if the state is named x
     for i ∈ 1:9
         p.aliases[Symbol(x, ctupperscripts(i))] = :($x^$i)
-    end # make: x¹, x²... if the state is named x
+    end # Make x¹, x²... if the state is named x
     p.aliases[Symbol(Unicode.normalize(string(x, "̇")))] = :(∂($x))
     if (isnothing(components_names))
-        __wrap(:(state!($ocp, $n, $xx)), p.lnum, p.line)
+        code = :( state!($ocp, $n, $xx) )
     else
         nn == length(components_names.args) ||
             return __throw("the number of state components must be $nn", p.lnum, p.line)
         for i ∈ 1:nn
             p.aliases[components_names.args[i]] = :($x[$i])
-        end # aliases from names given by the user
+            # todo: add aliases for state components (scalar) derivatives
+        end # Aliases from names given by the user
         ss = QuoteNode(string.(components_names.args))
-        __wrap(:(state!($ocp, $n, $xx, $ss)), p.lnum, p.line)
+        code = :( state!($ocp, $n, $xx, $ss) )
     end
+    return __wrap(code, p.lnum, p.line)
 end
 
 p_control!(p, ocp, u, m; components_names = nothing, log = false) = begin
@@ -321,7 +325,7 @@ p_control!(p, ocp, u, m; components_names = nothing, log = false) = begin
         p.aliases[Symbol(u, ctupperscripts(i))] = :($u^$i)
     end # make: u¹, u²... if the control is named u
     if (isnothing(components_names))
-        __wrap(:(control!($ocp, $m, $uu)), p.lnum, p.line)
+        code = :( control!($ocp, $m, $uu) )
     else
         mm == length(components_names.args) ||
             return __throw("the number of control components must be $mm", p.lnum, p.line)
@@ -329,8 +333,9 @@ p_control!(p, ocp, u, m; components_names = nothing, log = false) = begin
             p.aliases[components_names.args[i]] = :($u[$i])
         end # aliases from names given by the user
         ss = QuoteNode(string.(components_names.args))
-        __wrap(:(control!($ocp, $m, $uu, $ss)), p.lnum, p.line)
+        code = :( control!($ocp, $m, $uu, $ss) )
     end
+    return __wrap(code, p.lnum, p.line)
 end
 
 p_constraint!(p, ocp, e1, e2, e3, label = gensym(); log = false) = begin
@@ -348,13 +353,15 @@ p_constraint!(p, ocp, e1, e2, e3, label = gensym(); log = false) = begin
             gs = gensym()
             x0 = gensym()
             xf = gensym()
+            r = gensym()
             ee2 = replace_call(e2, p.x, p.t0, x0)
             ee2 = replace_call(ee2, p.x, p.tf, xf)
-            args = [x0, xf]
+            args = [r, x0, xf]
             __v_dep(p) && push!(args, p.v)
             quote
                 function $gs($(args...))
-                    $ee2
+                    @views $r[:] .= $ee2
+                    return nothing
                 end
                 constraint!($ocp, :boundary; f = $gs, lb = $e1, ub = $e3, label = $llabel)
             end
@@ -364,15 +371,17 @@ p_constraint!(p, ocp, e1, e2, e3, label = gensym(); log = false) = begin
         :control_fun => begin
             gs = gensym()
             ut = gensym()
+            r = gensym()
             ee2 = replace_call(e2, p.u, p.t, ut)
             p.t_dep = p.t_dep || has(ee2, p.t)
-            args = []
+            args = [r]
             __t_dep(p) && push!(args, p.t)
             push!(args, ut)
             __v_dep(p) && push!(args, p.v)
             quote
                 function $gs($(args...))
-                    $ee2
+                    @views $r[:] .= $ee2
+                    return nothing
                 end
                 constraint!($ocp, :control; f = $gs, lb = $e1, ub = $e3, label = $llabel)
             end
@@ -382,15 +391,17 @@ p_constraint!(p, ocp, e1, e2, e3, label = gensym(); log = false) = begin
         :state_fun => begin
             gs = gensym()
             xt = gensym()
+            r = gensym()
             ee2 = replace_call(e2, p.x, p.t, xt)
             p.t_dep = p.t_dep || has(ee2, p.t)
-            args = []
+            args = [r]
             __t_dep(p) && push!(args, p.t)
             push!(args, xt)
             __v_dep(p) && push!(args, p.v)
             quote
                 function $gs($(args...))
-                    $ee2
+                    @views $r[:] .= $ee2
+                    return nothing
                 end
                 constraint!($ocp, :state; f = $gs, lb = $e1, ub = $e3, label = $llabel)
             end
@@ -399,10 +410,12 @@ p_constraint!(p, ocp, e1, e2, e3, label = gensym(); log = false) = begin
             :(constraint!($ocp, :variable; rg = $rg, lb = $e1, ub = $e3, label = $llabel))
         :variable_fun => begin
             gs = gensym()
-            args = [p.v]
+            r = gensym()
+            args = [r, p.v]
             quote
                 function $gs($(args...))
-                    $e2
+                    @views $r[:] .= $e2
+                    return nothing
                 end
                 constraint!($ocp, :variable; f = $gs, lb = $e1, ub = $e3, label = $llabel)
             end
@@ -411,22 +424,24 @@ p_constraint!(p, ocp, e1, e2, e3, label = gensym(); log = false) = begin
             gs = gensym()
             xt = gensym()
             ut = gensym()
+            r = gensym()
             ee2 = replace_call(e2, [p.x, p.u], p.t, [xt, ut])
             p.t_dep = p.t_dep || has(ee2, p.t)
-            args = []
+            args = [r]
             __t_dep(p) && push!(args, p.t)
             push!(args, xt, ut)
             __v_dep(p) && push!(args, p.v)
             quote
                 function $gs($(args...))
-                    $ee2
+                    @views $r[:] .= $ee2
+                    return nothing
                 end
                 constraint!($ocp, :mixed; f = $gs, lb = $e1, ub = $e3, label = $llabel)
             end
         end
         _ => return __throw("bad constraint declaration", p.lnum, p.line)
     end
-    __wrap(code, p.lnum, p.line)
+    return __wrap(code, p.lnum, p.line)
 end
 
 p_dynamics!(p, ocp, x, t, e, label = nothing; log = false) = begin
@@ -443,16 +458,16 @@ p_dynamics!(p, ocp, x, t, e, label = nothing; log = false) = begin
     e = replace_call(e, [p.x, p.u], p.t, [xt, ut])
     p.t_dep = p.t_dep || has(e, t)
     gs = gensym()
-    args = []
-    __t_dep(p) && push!(args, p.t)
-    push!(args, xt, ut)
-    __v_dep(p) && push!(args, p.v)
-    __wrap(quote
+    r = gensym()
+    args = [r]; __t_dep(p) && push!(args, p.t); push!(args, xt, ut); __v_dep(p) && push!(args, p.v)
+    code = quote
         function $gs($(args...))
-            $e
+            @views $r[:] .= $e
+            return nothing
         end
         dynamics!($ocp, $gs)
-    end, p.lnum, p.line)
+    end
+    return __wrap(code, p.lnum, p.line)
 end
 
 p_lagrange!(p, ocp, e, type; log = false) = begin
@@ -466,16 +481,16 @@ p_lagrange!(p, ocp, e, type; log = false) = begin
     p.t_dep = p.t_dep || has(e, p.t)
     ttype = QuoteNode(type)
     gs = gensym()
-    args = []
-    __t_dep(p) && push!(args, p.t)
-    push!(args, xt, ut)
-    __v_dep(p) && push!(args, p.v)
-    __wrap(quote
+    r = gensym()
+    args = [r]; __t_dep(p) && push!(args, p.t); push!(args, xt, ut); __v_dep(p) && push!(args, p.v)
+    code = quote
         function $gs($(args...))
-            $e
+            @views $r[:] .= $e
+            return nothing
         end
         objective!($ocp, :lagrange, $gs, $ttype)
-    end, p.lnum, p.line)
+    end
+    return __wrap(code, p.lnum, p.line)
 end
 
 p_mayer!(p, ocp, e, type; log = false) = begin
@@ -491,17 +506,19 @@ p_mayer!(p, ocp, e, type; log = false) = begin
     gs = gensym()
     x0 = gensym()
     xf = gensym()
+    r = gensym()
     e = replace_call(e, p.x, p.t0, x0)
     e = replace_call(e, p.x, p.tf, xf)
     ttype = QuoteNode(type)
-    args = [x0, xf]
-    __v_dep(p) && push!(args, p.v)
-    __wrap(quote
+    args = [r, x0, xf]; __v_dep(p) && push!(args, p.v)
+    code = quote
         function $gs($(args...))
-            $e
+            @views $r[:] .= $e
+            return nothing
         end
         objective!($ocp, :mayer, $gs, $ttype)
-    end, p.lnum, p.line)
+    end
+    return __wrap(code, p.lnum, p.line)
 end
 
 p_bolza!(p, ocp, e1, e2, type; log = false) = begin
@@ -514,29 +531,34 @@ p_bolza!(p, ocp, e1, e2, type; log = false) = begin
     gs1 = gensym()
     x0 = gensym()
     xf = gensym()
+    r1 = gensym()
     e1 = replace_call(e1, p.x, p.t0, x0)
     e1 = replace_call(e1, p.x, p.tf, xf)
-    args1 = [x0, xf]
+    args1 = [r1, x0, xf]
     __v_dep(p) && push!(args1, p.v)
     gs2 = gensym()
     xt = gensym()
     ut = gensym()
+    r2 = gensym()
     e2 = replace_call(e2, [p.x, p.u], p.t, [xt, ut])
     p.t_dep = p.t_dep || has(e2, p.t)
-    args2 = []
+    args2 = [r2]
     __t_dep(p) && push!(args2, p.t)
     push!(args2, xt, ut)
     __v_dep(p) && push!(args2, p.v)
     ttype = QuoteNode(type)
-    __wrap(quote
+    code = quote
         function $gs1($(args1...))
-            $e1
+            $r1[:] .= $e1
+            return nothing
         end
         function $gs2($(args2...))
-            $e2
+            $r2[:] .= $e2
+            return nothing
         end
         objective!($ocp, :bolza, $gs1, $gs2, $ttype)
-    end, p.lnum, p.line)
+    end
+    return __wrap(code, p.lnum, p.line)
 end
 
 """
@@ -615,11 +637,12 @@ macro def(ocp, e, log = false)
         p.t_dep = p0.t_dep
         p.v = p0.v
         code = parse!(p, ocp, e; log = log)
+        in_place = true # todo: remove?
         init = @match (__t_dep(p), __v_dep(p)) begin
-            (false, false) => :($ocp = __OCPModel())
-            (true, false) => :($ocp = __OCPModel(autonomous = false))
-            (false, true) => :($ocp = __OCPModel(variable = true))
-            _ => :($ocp = __OCPModel(autonomous = false, variable = true))
+            (false, false) => :($ocp = __OCPModel(; in_place = $in_place))
+            (true, false) => :($ocp = __OCPModel(autonomous = false; in_place = $in_place))
+            (false, true) => :($ocp = __OCPModel(variable = true; in_place = $in_place))
+            _ => :($ocp = __OCPModel(autonomous = false, variable = true; in_place = $in_place))
         end
         ee = QuoteNode(e)
         code = Expr(:block, init, code, :($ocp.model_expression = $ee; $ocp))
@@ -627,4 +650,4 @@ macro def(ocp, e, log = false)
     catch ex
         :(throw($ex)) # can be caught by user
     end
-end
+end
\ No newline at end of file
diff --git a/src/optimal_control_model-getters.jl b/src/optimal_control_model-getters.jl
index b7c74285..7b22695e 100644
--- a/src/optimal_control_model-getters.jl
+++ b/src/optimal_control_model-getters.jl
@@ -1,6 +1,6 @@
 # ----------------------------------------------------------------------
 #
-# Interaction with the Model that get data. Getters.
+# Model getters.
 #
 
 """
@@ -17,7 +17,8 @@ Return a 6-tuple of tuples:
 - `(xl, xind, xu)` are state linear constraints of a subset of indices
 - `(vl, vind, vu)` are variable linear constraints of a subset of indices
 
-and update information about constraints dimensions of  `ocp`.
+and update information about constraints dimensions of  `ocp`. Functions `ξ`, `η`, `ψ`, `ϕ`, `θ` are used to evaluate the constraints at given time and state, control, or variable values. These functions are in place when the problem is
+defined as in place, that is when `is_in_place(ocp)`.
 
 !!! note
 
@@ -36,41 +37,60 @@ function nlp_constraints!(ocp::OptimalControlModel)
     # we check if the dimensions and times have been set
     __check_all_set(ocp)
 
-    ξf = Vector{ControlConstraint}()
+    ControlConstraint_ = is_in_place(ocp) ? ControlConstraint! : ControlConstraint
+    StateConstraint_ = is_in_place(ocp) ? StateConstraint! : StateConstraint
+    MixedConstraint_ = is_in_place(ocp) ? MixedConstraint! : MixedConstraint
+    BoundaryConstraint_ = is_in_place(ocp) ? BoundaryConstraint! : BoundaryConstraint
+    VariableConstraint_ = is_in_place(ocp) ? VariableConstraint! : VariableConstraint
+
+    ξf = Vector{ControlConstraint_}()
+    ξs = Vector{Int64}()
     ξl = Vector{ctNumber}()
     ξu = Vector{ctNumber}()
-    ηf = Vector{StateConstraint}()
+
+    ηf = Vector{StateConstraint_}()
+    ηs = Vector{Int64}()
     ηl = Vector{ctNumber}()
     ηu = Vector{ctNumber}()
-    ψf = Vector{MixedConstraint}()
+
+    ψf = Vector{MixedConstraint_}()
+    ψs = Vector{Int64}()
     ψl = Vector{ctNumber}()
     ψu = Vector{ctNumber}()
-    ϕf = Vector{BoundaryConstraint}()
+
+    ϕf = Vector{BoundaryConstraint_}()
+    ϕs = Vector{Int64}()
     ϕl = Vector{ctNumber}()
     ϕu = Vector{ctNumber}()
-    θf = Vector{VariableConstraint}()
+
+    θf = Vector{VariableConstraint_}()
+    θs = Vector{Int64}()
     θl = Vector{ctNumber}()
     θu = Vector{ctNumber}()
+
     uind = Vector{Int}()
     ul = Vector{ctNumber}()
     uu = Vector{ctNumber}()
+
     xind = Vector{Int}()
     xl = Vector{ctNumber}()
     xu = Vector{ctNumber}()
+
     vind = Vector{Int}()
     vl = Vector{ctNumber}()
     vu = Vector{ctNumber}()
 
     for (_, c) ∈ constraints(ocp)
         @match c begin
-            (type, f::BoundaryConstraint, lb, ub) && if type ∈ [:initial, :final, :boundary]
-            end => begin
+            (type, f::BoundaryConstraint_, lb, ub) && if type ∈ [:initial, :final, :boundary] end => begin
                 push!(ϕf, f)
+                push!(ϕs, length(lb))
                 append!(ϕl, lb)
                 append!(ϕu, ub)
             end
-            (:control, f::ControlConstraint, lb, ub) => begin
+            (:control, f::ControlConstraint_, lb, ub) => begin
                 push!(ξf, f)
+                push!(ξs, length(lb))
                 append!(ξl, lb)
                 append!(ξu, ub)
             end
@@ -79,8 +99,9 @@ function nlp_constraints!(ocp::OptimalControlModel)
                 append!(ul, lb)
                 append!(uu, ub)
             end
-            (:state, f::StateConstraint, lb, ub) => begin
+            (:state, f::StateConstraint_, lb, ub) => begin
                 push!(ηf, f)
+                push!(ηs, length(lb))
                 append!(ηl, lb)
                 append!(ηu, ub)
             end
@@ -89,13 +110,15 @@ function nlp_constraints!(ocp::OptimalControlModel)
                 append!(xl, lb)
                 append!(xu, ub)
             end
-            (:mixed, f::MixedConstraint, lb, ub) => begin
+            (:mixed, f::MixedConstraint_, lb, ub) => begin
                 push!(ψf, f)
+                push!(ψs, length(lb))
                 append!(ψl, lb)
                 append!(ψu, ub)
             end
-            (:variable, f::VariableConstraint, lb, ub) => begin
+            (:variable, f::VariableConstraint_, lb, ub) => begin
                 push!(θf, f)
+                push!(θs, length(lb))
                 append!(θl, lb)
                 append!(θu, ub)
             end
@@ -109,10 +132,15 @@ function nlp_constraints!(ocp::OptimalControlModel)
     end
 
     @assert length(ξl) == length(ξu)
+    @assert length(ξf) == length(ξs)
     @assert length(ηl) == length(ηu)
+    @assert length(ηf) == length(ηs)
     @assert length(ψl) == length(ψu)
+    @assert length(ψf) == length(ψs)
     @assert length(ϕl) == length(ϕu)
+    @assert length(ϕf) == length(ϕs)
     @assert length(θl) == length(θu)
+    @assert length(θf) == length(θs)
     @assert length(ul) == length(uu)
     @assert length(xl) == length(xu)
     @assert length(vl) == length(vu)
@@ -122,66 +150,116 @@ function nlp_constraints!(ocp::OptimalControlModel)
         val = zeros(ctNumber, dim)
         j = 1
         for i ∈ 1:length(ξf)
+            li = ξs[i]
             vali = ξf[i](t, u, v)
-            li = length(vali)
             val[j:(j + li - 1)] .= vali # .= also allows scalar value for vali
             j = j + li
         end
         return val
     end
 
+    function ξ!(val, t, u, v) # nonlinear control constraints (in place)
+        j = 1
+        for i ∈ 1:length(ξf)
+            li = ξs[i]
+            ξf[i](@view(val[j:(j + li - 1)]), t, u, v)
+            j = j + li
+        end
+        return nothing
+    end
+
     function η(t, x, v) # nonlinear state constraints
         dim = length(ηl)
         val = zeros(ctNumber, dim)
         j = 1
         for i ∈ 1:length(ηf)
+            li = ηs[i]
             vali = ηf[i](t, x, v)
-            li = length(vali)
             val[j:(j + li - 1)] .= vali # .= also allows scalar value for vali
             j = j + li
         end
         return val
     end
 
+    function η!(val, t, x, v) # nonlinear state constraints (in place)
+        j = 1
+        for i ∈ 1:length(ηf)
+            li = ηs[i]
+            ηf[i](@view(val[j:(j + li - 1)]), t, x, v)
+            j = j + li
+        end
+        return nothing
+    end
+
     function ψ(t, x, u, v) # nonlinear mixed constraints
         dim = length(ψl)
         val = zeros(ctNumber, dim)
         j = 1
         for i ∈ 1:length(ψf)
+            li = ψs[i]
             vali = ψf[i](t, x, u, v)
-            li = length(vali)
             val[j:(j + li - 1)] .= vali # .= also allows scalar value for vali
             j = j + li
         end
         return val
     end
 
+    function ψ!(val, t, x, u, v) # nonlinear mixed constraints (in place)
+        j = 1
+        for i ∈ 1:length(ψf)
+            li = ψs[i]
+            ψf[i](@view(val[j:(j + li - 1)]), t, x, u, v)
+            j = j + li
+        end
+        return nothing
+    end
+
     function ϕ(x0, xf, v) # nonlinear boundary constraints
         dim = length(ϕl)
         val = zeros(ctNumber, dim)
         j = 1
         for i ∈ 1:length(ϕf)
+            li = ϕs[i]
             vali = ϕf[i](x0, xf, v)
-            li = length(vali)
             val[j:(j + li - 1)] .= vali # .= also allows scalar value for vali
             j = j + li
         end
         return val
     end
 
+    function ϕ!(val, x0, xf, v) # nonlinear boundary constraints
+        j = 1
+        for i ∈ 1:length(ϕf)
+            li = ϕs[i]
+            ϕf[i](@view(val[j:(j + li - 1)]), x0, xf, v)
+            j = j + li
+        end
+        return nothing
+    end
+
     function θ(v) # nonlinear variable constraints
         dim = length(θl)
         val = zeros(ctNumber, dim)
         j = 1
         for i ∈ 1:length(θf)
+            li = θs[i]
             vali = θf[i](v)
-            li = length(vali)
             val[j:(j + li - 1)] .= vali # .= also allows scalar value for vali
             j = j + li
         end
         return val
     end
 
+    function θ!(val, v) # nonlinear variable constraints
+        j = 1
+        for i ∈ 1:length(θf)
+            li = θs[i]
+            θf[i](@view(val[j:(j + li - 1)]), v)
+            j = j + li
+        end
+        return nothing
+    end
+
     ocp.dim_control_constraints = length(ξl)
     ocp.dim_state_constraints = length(ηl)
     ocp.dim_mixed_constraints = length(ψl)
@@ -191,14 +269,20 @@ function nlp_constraints!(ocp::OptimalControlModel)
     ocp.dim_state_range = length(xl)
     ocp.dim_variable_range = length(vl)
 
-    return (ξl, ξ, ξu),
-    (ηl, η, ηu),
-    (ψl, ψ, ψu),
-    (ϕl, ϕ, ϕu),
-    (θl, θ, θu),
-    (ul, uind, uu),
-    (xl, xind, xu),
-    (vl, vind, vu)
+    ξ_ = is_in_place(ocp) ? ξ! : ξ
+    η_ = is_in_place(ocp) ? η! : η
+    ψ_ = is_in_place(ocp) ? ψ! : ψ
+    ϕ_ = is_in_place(ocp) ? ϕ! : ϕ
+    θ_ = is_in_place(ocp) ? θ! : θ
+
+    return (ξl, ξ_, ξu),
+           (ηl, η_, ηu),
+           (ψl, ψ_, ψu),
+           (ϕl, ϕ_, ϕu),
+           (θl, θ_, θu),
+           (ul, uind, uu),
+           (xl, xind, xu),
+           (vl, vind, vu)
 end
 
 """
@@ -241,41 +325,59 @@ function constraint(
     label::Symbol,
 ) where {T <: TimeDependence, V <: VariableDependence}
     con = ocp.constraints[label]
+    BoundaryConstraint_ = is_in_place(ocp) ? BoundaryConstraint! : BoundaryConstraint
+    ControlConstraint_ = is_in_place(ocp) ? ControlConstraint! : ControlConstraint
+    StateConstraint_ = is_in_place(ocp) ? StateConstraint! : StateConstraint
+    MixedConstraint_ = is_in_place(ocp) ? MixedConstraint! : MixedConstraint
+    VariableConstraint_ = is_in_place(ocp) ? VariableConstraint! : VariableConstraint
+
     @match con begin
-        (:initial, f::BoundaryConstraint, _, _) => return f
-        (:final, f::BoundaryConstraint, _, _) => return f
-        (:boundary, f::BoundaryConstraint, _, _) => return f
-        (:control, f::ControlConstraint, _, _) => return f
+        (:initial, f::BoundaryConstraint_, _, _) => return f
+        (:final, f::BoundaryConstraint_, _, _) => return f
+        (:boundary, f::BoundaryConstraint_, _, _) => return f
+        (:control, f::ControlConstraint_, _, _) => return f
         (:control, rg, _, _) => begin
             C = @match ocp begin
-                ::OptimalControlModel{Autonomous, Fixed} => ControlConstraint(u -> u[rg], T, V)
-                ::OptimalControlModel{Autonomous, NonFixed} =>
+                ::OptimalControlModel{Autonomous, Fixed} => is_in_place(ocp) ?
+                    ControlConstraint!((r, u) -> (@views r[:] .= u[rg]; nothing), T, V) : # todo: CC!{T, V}(fun) syntax?
+                    ControlConstraint(u -> u[rg], T, V)
+                ::OptimalControlModel{Autonomous, NonFixed} => is_in_place(ocp) ?
+                    ControlConstraint!((r, u, v) -> (@views r[:] .= u[rg]; nothing), T, V) :
                     ControlConstraint((u, v) -> u[rg], T, V)
-                ::OptimalControlModel{NonAutonomous, Fixed} =>
+                ::OptimalControlModel{NonAutonomous, Fixed} => is_in_place(ocp) ?
+                    ControlConstraint!((r, t, u) -> (@views r[:] .= u[rg]; nothing), T, V) :
                     ControlConstraint((t, u) -> u[rg], T, V)
-                ::OptimalControlModel{NonAutonomous, NonFixed} =>
+                ::OptimalControlModel{NonAutonomous, NonFixed} => is_in_place(ocp) ?
+                    ControlConstraint!((r, t, u, v) -> (@views r[:] .= u[rg]; nothing), T, V) :
                     ControlConstraint((t, u, v) -> u[rg], T, V)
                 _ => nothing
             end
             return C
         end
-        (:state, f::StateConstraint, _, _) => return f
+        (:state, f::StateConstraint_, _, _) => return f
         (:state, rg, _, _) => begin
             S = @match ocp begin
-                ::OptimalControlModel{Autonomous, Fixed} => StateConstraint(x -> x[rg], T, V)
-                ::OptimalControlModel{Autonomous, NonFixed} =>
+                ::OptimalControlModel{Autonomous, Fixed} => is_in_place(ocp) ?
+                    StateConstraint!((r, x) -> (@views r[:] .= x[rg]; nothing), T, V) :
+                    StateConstraint(x -> x[rg], T, V)
+                ::OptimalControlModel{Autonomous, NonFixed} => is_in_place(ocp) ?
+                    StateConstraint!((r, x, v) -> (@views r[:] .= x[rg]; nothing), T, V) :
                     StateConstraint((x, v) -> x[rg], T, V)
-                ::OptimalControlModel{NonAutonomous, Fixed} =>
+                ::OptimalControlModel{NonAutonomous, Fixed} => is_in_place(ocp) ?
+                    StateConstraint!((r, t, x) -> (@views r[:] .= x[rg]; nothing), T, V) :
                     StateConstraint((t, x) -> x[rg], T, V)
-                ::OptimalControlModel{NonAutonomous, NonFixed} =>
+                ::OptimalControlModel{NonAutonomous, NonFixed} => is_in_place(ocp) ?
+                    StateConstraint!((r, t, x, v) -> (@views r[:] .= x[rg]; nothing), T, V) :
                     StateConstraint((t, x, v) -> x[rg], T, V)
                 _ => nothing
             end
             return S
         end
-        (:mixed, f::MixedConstraint, _, _) => return f
-        (:variable, f::VariableConstraint, _, _) => return f
-        (:variable, rg, _, _) => return VariableConstraint(v -> v[rg])
+        (:mixed, f::MixedConstraint_, _, _) => return f
+        (:variable, f::VariableConstraint_, _, _) => return f
+        (:variable, rg, _, _) => return ( is_in_place(ocp) ?
+            VariableConstraint!((r, v) -> (@views r[:] .= v[rg]; nothing)) :
+            VariableConstraint(v -> v[rg]) )
         _ => error("Internal error")
     end
 end
@@ -570,6 +672,13 @@ is_variable_independent(ocp::OptimalControlModel) = !is_variable_dependent(ocp)
 """
 $(TYPEDSIGNATURES)
 
+Return `true` if functions defining the ocp are in-place. Return nothing if this information has not yet been set.
+"""
+is_in_place(ocp::OptimalControlModel) = ocp.in_place
+
+"""
+$(TYPEDSIGNATURES)
+
 Return `true` if the model has been defined with free initial time.
 """
 has_free_initial_time(ocp::OptimalControlModel) = (typeof(initial_time(ocp)) == Index)
diff --git a/src/optimal_control_model-setters.jl b/src/optimal_control_model-setters.jl
index 50849ea8..29a3991d 100644
--- a/src/optimal_control_model-setters.jl
+++ b/src/optimal_control_model-setters.jl
@@ -3,8 +3,6 @@
 # Interaction with the Model that affect it. Setters / Constructors.
 #
 
-# todo: use design pattern to generate functions in nlp_constraints!
-
 """
 $(TYPEDSIGNATURES)
 
@@ -37,11 +35,13 @@ julia> ocp = Model(autonomous=false, variable=true)
     - If the time dependence of the model is defined as nonautonomous, then, the dynamics function, the lagrange cost and the path constraints must be defined as functions of time and state, and possibly control. If the model is defined as autonomous, then, the dynamics function, the lagrange cost and the path constraints must be defined as functions of state, and possibly control.
 
 """
-function Model(; autonomous::Bool = true, variable::Bool = false)
+function Model(; autonomous::Bool = true, variable::Bool = false, in_place::Bool = false)
     time_dependence = autonomous ? Autonomous : NonAutonomous
     variable_dependence = variable ? NonFixed : Fixed
+    ocp = OptimalControlModel{time_dependence, variable_dependence}()
+    ocp.in_place = in_place
 
-    return OptimalControlModel{time_dependence, variable_dependence}()
+    return ocp
 end
 
 """
@@ -67,13 +67,13 @@ julia> ocp = Model(Autonomous, NonFixed)
 
 """
 function Model(
-    dependencies::DataType...,
+    dependencies::DataType...; in_place::Bool = false,
 )::OptimalControlModel{<:TimeDependence, <:VariableDependence}
     # some checkings: 
     __check_dependencies(dependencies)
     time_dependence = NonAutonomous ∈ dependencies ? NonAutonomous : Autonomous
     variable_dependence = NonFixed ∈ dependencies ? NonFixed : Fixed
-    return OptimalControlModel{time_dependence, variable_dependence}()
+    return OptimalControlModel{time_dependence, variable_dependence}(; in_place = in_place)
 end
 
 """
@@ -567,13 +567,19 @@ function constraint!(
     q = variable_dimension(ocp)
 
     # range
-    (typeof(rg) <: Int) && (rg = Index(rg))
+    (typeof(rg) <: Int) && (rg = Index(rg)) # todo: scalar range
 
     # core
+    BoundaryConstraint_ = is_in_place(ocp) ? BoundaryConstraint! : BoundaryConstraint
+    ControlConstraint_ = is_in_place(ocp) ? ControlConstraint! : ControlConstraint
+    StateConstraint_ = is_in_place(ocp) ? StateConstraint! : StateConstraint
+    MixedConstraint_ = is_in_place(ocp) ? MixedConstraint! : MixedConstraint
+    VariableConstraint_ = is_in_place(ocp) ? VariableConstraint! : VariableConstraint
+
     @match (rg, f, lb, ub) begin
         (::Nothing, ::Nothing, ::ctVector, ::ctVector) => begin
             if type ∈ [:initial, :final, :state]
-                rg = n == 1 ? Index(1) : 1:n
+                rg = n == 1 ? Index(1) : 1:n # todo: scalar range
                 txt = "the lower bound `lb`, the upper bound `ub` and the value `val` must be of dimension $n"
             elseif type == :control
                 rg = m == 1 ? Index(1) : 1:m
@@ -631,14 +637,25 @@ function constraint!(
                     ),
                 )
             end
+
             # set the constraint
             fun_rg = @match type begin
                 :initial =>
-                    V == Fixed ? BoundaryConstraint((x0, xf) -> x0[rg], V) :
-                    BoundaryConstraint((x0, xf, v) -> x0[rg], V)
+                    if is_in_place(ocp)
+                        V == Fixed ? BoundaryConstraint!((r, x0, xf) -> (@views r[:] .= x0[rg]; nothing), V) :
+                        BoundaryConstraint!((r, x0, xf, v) -> (@views r[:] .= x0[rg]; nothing), V)
+                    else
+                        V == Fixed ? BoundaryConstraint((x0, xf) -> x0[rg], V) :
+                        BoundaryConstraint((x0, xf, v) -> x0[rg], V)
+                    end
                 :final =>
-                    V == Fixed ? BoundaryConstraint((x0, xf) -> xf[rg], V) :
-                    BoundaryConstraint((x0, xf, v) -> xf[rg], V)
+                    if is_in_place(ocp)
+                        V == Fixed ? BoundaryConstraint!((r, x0, xf) -> (@views r[:] .= xf[rg]; nothing), V) :
+                        BoundaryConstraint!((r, x0, xf, v) -> (@views r[:] .= xf[rg]; nothing), V)
+                    else
+                        V == Fixed ? BoundaryConstraint((x0, xf) -> xf[rg], V) :
+                        BoundaryConstraint((x0, xf, v) -> xf[rg], V)
+                    end
                 :control || :state || :variable => rg
                 _ => throw(
                     IncorrectArgument(
@@ -651,18 +668,18 @@ function constraint!(
             ocp.constraints[label] = (type, fun_rg, lb, ub)
         end
 
-        (::Nothing, ::Function, ::ctVector, ::ctVector) => begin
-            # set the constraint
+        (::Nothing, ::Function, ::ctVector, ::ctVector) => begin 
+                       # set the constraint
             if type == :boundary
-                ocp.constraints[label] = (type, BoundaryConstraint(f, V), lb, ub)
+                ocp.constraints[label] = (type, BoundaryConstraint_(f, V), lb, ub)
             elseif type == :control
-                ocp.constraints[label] = (type, ControlConstraint(f, T, V), lb, ub)
+                  ocp.constraints[label] = (type, ControlConstraint_(f, T, V), lb, ub)
             elseif type == :state
-                ocp.constraints[label] = (type, StateConstraint(f, T, V), lb, ub)
+                ocp.constraints[label] = (type, StateConstraint_(f, T, V), lb, ub)
             elseif type == :mixed
-                ocp.constraints[label] = (type, MixedConstraint(f, T, V), lb, ub)
+                ocp.constraints[label] = (type, MixedConstraint_(f, T, V), lb, ub)
             elseif type == :variable
-                ocp.constraints[label] = (type, VariableConstraint(f), lb, ub)
+                ocp.constraints[label] = (type, VariableConstraint_(f), lb, ub)
             else
                 throw(
                     IncorrectArgument(
@@ -707,7 +724,8 @@ function dynamics!(
     __check_all_set(ocp)
     __is_dynamics_set(ocp) && throw(UnauthorizedCall("the dynamics has already been set."))
 
-    ocp.dynamics = Dynamics(f, T, V)
+    Dynamics_ = is_in_place(ocp) ? Dynamics! : Dynamics
+    ocp.dynamics = Dynamics_(f, T, V)
 
     return nothing
 end
@@ -758,10 +776,12 @@ function objective!(
     ocp.criterion = criterion
 
     # set the objective
+    Mayer_ = is_in_place(ocp) ? Mayer! : Mayer
+    Lagrange_ = is_in_place(ocp) ? Lagrange! : Lagrange
     if type == :mayer
-        ocp.mayer = Mayer(f, V)
+        ocp.mayer = Mayer_(f, V)
     elseif type == :lagrange
-        ocp.lagrange = Lagrange(f, T, V)
+        ocp.lagrange = Lagrange_(f, T, V)
     else
         throw(
             IncorrectArgument(
@@ -817,9 +837,11 @@ function objective!(
     ocp.criterion = criterion
 
     # set the objective
+    Mayer_ = is_in_place(ocp) ? Mayer! : Mayer
+    Lagrange_ = is_in_place(ocp) ? Lagrange! : Lagrange
     if type == :bolza
-        ocp.mayer = Mayer(g, V)
-        ocp.lagrange = Lagrange(f⁰, T, V)
+        ocp.mayer = Mayer_(g, V) 
+        ocp.lagrange = Lagrange_(f⁰, T, V)
     else
         throw(
             IncorrectArgument(
diff --git a/src/optimal_control_model-type.jl b/src/optimal_control_model-type.jl
index 65aff302..230cfffe 100644
--- a/src/optimal_control_model-type.jl
+++ b/src/optimal_control_model-type.jl
@@ -28,10 +28,10 @@ $(TYPEDFIELDS)
     variable_dimension::Union{Dimension, Nothing} = nothing
     variable_components_names::Union{Vector{String}, Nothing} = nothing
     variable_name::Union{String, Nothing} = nothing
-    lagrange::Union{Lagrange, Nothing} = nothing
-    mayer::Union{Mayer, Nothing} = nothing
+    lagrange::Union{Lagrange, Lagrange!, Nothing} = nothing
+    mayer::Union{Mayer, Mayer!, Nothing} = nothing
     criterion::Union{Symbol, Nothing} = nothing
-    dynamics::Union{Dynamics, Nothing} = nothing
+    dynamics::Union{Dynamics, Dynamics!, Nothing} = nothing
     constraints::Dict{Symbol, Tuple{Vararg{Any}}} = Dict{Symbol, Tuple{Vararg{Any}}}()
     dim_control_constraints::Union{Dimension, Nothing} = nothing
     dim_state_constraints::Union{Dimension, Nothing} = nothing
@@ -41,6 +41,7 @@ $(TYPEDFIELDS)
     dim_control_range::Union{Dimension, Nothing} = nothing
     dim_state_range::Union{Dimension, Nothing} = nothing
     dim_variable_range::Union{Dimension, Nothing} = nothing
+    in_place::Union{Bool, Nothing} = nothing 
 end
 
 # ----------------------------------------------------------------------
diff --git a/src/types.jl b/src/types.jl
index 82f76797..2c88d640 100644
--- a/src/types.jl
+++ b/src/types.jl
@@ -57,6 +57,10 @@ struct BoundaryConstraint{variable_dependence}
     f::Function
 end
 
+struct BoundaryConstraint!{variable_dependence}
+    f!::Function
+end
+
 """
 $(TYPEDEF)
 
@@ -107,6 +111,10 @@ struct Mayer{variable_dependence}
     f::Function
 end
 
+struct Mayer!{variable_dependence}
+    f!::Function
+end
+
 """
 $(TYPEDEF)
 
@@ -117,7 +125,7 @@ abstract type AbstractHamiltonian{time_dependence, variable_dependence} end
 """
 $(TYPEDEF)
 
-Abstract type for vectorfields.
+Abstract type for vector fields.
 """
 abstract type AbstractVectorField{time_dependence, variable_dependence} end
 
@@ -508,6 +516,10 @@ struct Lagrange{time_dependence, variable_dependence}
     f::Function
 end
 
+struct Lagrange!{time_dependence, variable_dependence}
+    f!::Function
+end
+
 """
 $(TYPEDEF)
 
@@ -580,6 +592,10 @@ struct Dynamics{time_dependence, variable_dependence}
     f::Function
 end
 
+struct Dynamics!{time_dependence, variable_dependence}
+    f!::Function
+end
+
 """
 $(TYPEDEF)
 
@@ -658,6 +674,10 @@ struct StateConstraint{time_dependence, variable_dependence}
     f::Function
 end
 
+struct StateConstraint!{time_dependence, variable_dependence}
+    f!::Function
+end
+
 """
 $(TYPEDEF)
 
@@ -730,6 +750,10 @@ struct ControlConstraint{time_dependence, variable_dependence}
     f::Function
 end
 
+struct ControlConstraint!{time_dependence, variable_dependence}
+    f!::Function
+end
+
 """
 $(TYPEDEF)
 
@@ -810,6 +834,10 @@ struct MixedConstraint{time_dependence, variable_dependence}
     f::Function
 end
 
+struct MixedConstraint!{time_dependence, variable_dependence}
+    f!::Function
+end
+
 """
 $(TYPEDEF)
 
@@ -851,6 +879,10 @@ struct VariableConstraint
     f::Function
 end
 
+struct VariableConstraint!
+    f!::Function
+end
+
 """
 $(TYPEDEF)
 
diff --git a/src/utils.jl b/src/utils.jl
index 2af577ad..6a4ac54d 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -198,7 +198,6 @@ Return `expand(matrix2vec(x, 1))`
 """
 expand(x::Matrix{<:ctNumber}) = expand(matrix2vec(x, 1))
 
-# transform a Matrix{<:ctNumber} to a Vector{<:Vector{<:ctNumber}}
 """
 $(TYPEDSIGNATURES)
 
@@ -225,3 +224,32 @@ function matrix2vec(
     end
     return y
 end
+
+"""
+$(TYPEDSIGNATURES)
+
+Tranform in place function to out of place. Pass the result size and type (default = `Float64`).
+Return a scalar when the result has size one. If `f!` is `nothing`, return `nothing`.
+"""
+function to_out_of_place(f!, n; T = Float64)
+    function f(args...; kwargs...)
+        r = zeros(T, n)
+        f!(r, args...; kwargs...)
+        return n == 1 ? r[1] : r
+    end
+    return isnothing(f!) ? nothing : f
+end
+
+# Adapt getters to test in place
+function __constraint(ocp, label)
+    if is_in_place(ocp)
+        n = length(constraints(ocp)[label][3]) # Size of lb 
+        return to_out_of_place(constraint(ocp, label), n)    
+    else
+        return constraint(ocp, label)
+    end
+end
+
+__dynamics(ocp) = is_in_place(ocp) ? to_out_of_place(dynamics(ocp), state_dimension(ocp)) : dynamics(ocp)
+__lagrange(ocp) = is_in_place(ocp) ? to_out_of_place(lagrange(ocp), 1) : lagrange(ocp)
+__mayer(ocp) = is_in_place(ocp) ? to_out_of_place(mayer(ocp), 1) : mayer(ocp)
\ No newline at end of file
diff --git a/test/test_goddard.jl b/test/test_goddard.jl
index 3bac4ad4..65cb1ed8 100644
--- a/test/test_goddard.jl
+++ b/test/test_goddard.jl
@@ -53,5 +53,5 @@ test_goddard() = begin
     x = x0
     u = 2
     tf = 1
-    @test dynamics(ocp)(x, u, tf) == F0(x) + u * F1(x)
+    @test __dynamics(ocp)(x, u, tf) == F0(x) + u * F1(x)
 end
diff --git a/test/test_model.jl b/test/test_model.jl
index ae025f6f..2be8b97b 100644
--- a/test/test_model.jl
+++ b/test/test_model.jl
@@ -1,6 +1,43 @@
 function test_model() # 30 55 185
     ∅ = Vector{Real}()
 
+    @testset "is_in_place" begin
+        ocp = Model()
+        @test !is_in_place(ocp)
+        ocp = Model(; in_place = false)
+        @test !is_in_place(ocp)
+        ocp = Model(; in_place = true)
+        @test is_in_place(ocp)
+
+        ocp = Model(; autonomous = true, variable = true)
+        @test !is_in_place(ocp)
+        ocp = Model(; autonomous = true, variable = true, in_place = false)
+        @test !is_in_place(ocp)
+        ocp = Model(; autonomous = true, variable = true, in_place = true)
+        @test is_in_place(ocp)
+
+        ocp = Model(; autonomous = false, variable = true)
+        @test !is_in_place(ocp)
+        ocp = Model(; autonomous = false, variable = true, in_place = false)
+        @test !is_in_place(ocp)
+        ocp = Model(; autonomous = false, variable = true, in_place = true)
+        @test is_in_place(ocp)
+
+        ocp = Model(; autonomous = true, variable = false)
+        @test !is_in_place(ocp)
+        ocp = Model(; autonomous = true, variable = false, in_place = false)
+        @test !is_in_place(ocp)
+        ocp = Model(; autonomous = true, variable = false, in_place = true)
+        @test is_in_place(ocp)
+
+        ocp = Model(; autonomous = false, variable = false)
+        @test !is_in_place(ocp)
+        ocp = Model(; autonomous = false, variable = false, in_place = false)
+        @test !is_in_place(ocp)
+        ocp = Model(; autonomous = false, variable = false, in_place = true)
+        @test is_in_place(ocp)
+    end
+
     @testset "variable!" begin
         ocp = Model(variable = false)
 
@@ -1445,6 +1482,207 @@ function test_model() # 30 55 185
         @test dim_variable_range(ocp) == 7
     end
 
+
+    @testset "nlp_constraints! without variable (in place)" begin
+        ocp = Model(; in_place = true)
+        time!(ocp; t0 = 0, tf = 1)
+        state!(ocp, 2)
+        control!(ocp, 1)
+        constraint!(ocp, :initial, rg = 2, lb = 10, ub = 10, label = :ci)
+        constraint!(ocp, :final, rg = 1, lb = 1, ub = 1, label = :cf)
+        constraint!(ocp, :control, lb = 0, ub = 1, label = :cu)
+        constraint!(ocp, :state, lb = [0, 1], ub = [1, 2], label = :cs)
+        constraint!(ocp, :boundary, f = (r, x0, xf) -> (r[:] .= x0[2] + xf[2]; nothing), lb = 0, ub = 1, label = :cb)
+        constraint!(ocp, :control, f = (r, u) -> (r[:] .= u; nothing), lb = 0, ub = 1, label = :cuu)
+        constraint!(ocp, :state, f = (r, x) -> (r[:] .= x; nothing), lb = [0, 1], ub = [1, 2], label = :css)
+        constraint!(ocp, :mixed, f = (r, x, u) -> (r[:] .= x[1] + u; nothing), lb = 1, ub = 1, label = :cm)
+
+        # dimensions (not set yet)
+        @test dim_control_constraints(ocp) === nothing
+        @test dim_state_constraints(ocp) === nothing
+        @test dim_mixed_constraints(ocp) === nothing
+        @test dim_path_constraints(ocp) === nothing
+        @test dim_boundary_constraints(ocp) === nothing
+        @test dim_variable_constraints(ocp) === nothing
+        @test dim_control_range(ocp) === nothing
+        @test dim_state_range(ocp) === nothing
+        @test dim_variable_range(ocp) === nothing
+
+        (ξl, ξ!, ξu),
+        (ηl, η!, ηu),
+        (ψl, ψ!, ψu),
+        (ϕl, ϕ!, ϕu),
+        (θl, θ!, θu),
+        (ul, uind, uu),
+        (xl, xind, xu),
+        (vl, vind, vu) = nlp_constraints!(ocp)
+
+        v = Real[]
+
+        # control
+        @test sort(ξl) == sort([0])
+        @test sort(ξu) == sort([1])
+        r = [0.]
+        ξ!(r, -1, 1, v)
+        @test sort(r) == sort([1])
+
+        # state
+        @test sort(ηl) == sort([0, 1])
+        @test sort(ηu) == sort([1, 2])
+        r = [0., 0.]
+        η!(r, -1, [1, 1], v)
+        @test sort(r) == sort([1, 1])
+
+        # mixed
+        @test sort(ψl) == sort([1])
+        @test sort(ψu) == sort([1])
+        r = [0.]
+        ψ!(r, -1, [1, 1], 2, v)
+        @test sort(r) == sort([3])
+
+        # boundary
+        @test sort(ϕl) == sort([10, 1, 0])
+        @test sort(ϕu) == sort([10, 1, 1])
+        r = [0., 0., 0.]
+        ϕ!(r, [1, 3], [4, 100], v)
+        @test sort(r) == sort([3, 4, 103])
+
+        # box constraint
+        @test sort(ul) == sort([0])
+        @test sort(uind) == sort([1])
+        @test sort(uu) == sort([1])
+        @test sort(xl) == sort([0, 1])
+        @test sort(xind) == sort([1, 2])
+        @test sort(xu) == sort([1, 2])
+
+        # variable
+        @test sort(vl) == sort([])
+        @test sort(vind) == sort([])
+        @test sort(vu) == sort([])
+        @test sort(θl) == sort([])
+        @test sort(θu) == sort([])
+        r = Real[]
+        θ!(r, v)
+        @test sort(r) == sort([])
+
+        # dimensions (set)
+        @test dim_control_constraints(ocp) == 1
+        @test dim_state_constraints(ocp) == 2
+        @test dim_mixed_constraints(ocp) == 1
+        @test dim_path_constraints(ocp) == 4
+        @test dim_boundary_constraints(ocp) == 3
+        @test dim_variable_constraints(ocp) == 0
+        @test dim_control_range(ocp) == 1
+        @test dim_state_range(ocp) == 2
+        @test dim_variable_range(ocp) == 0
+    end
+
+    @testset "nlp_constraints! with variable (in place)" begin
+        ocp = Model(; variable = true, in_place = true)
+        time!(ocp; t0 = 0, tf = 1)
+        variable!(ocp, 4)
+        state!(ocp, 2)
+        control!(ocp, 1)
+        constraint!(ocp, :initial, rg = 2, lb = 10, ub = 10, label = :ci)
+        constraint!(ocp, :final, rg = 1, lb = 1, ub = 1, label = :cf)
+        constraint!(ocp, :control, lb = 0, ub = 1, label = :cu)
+        constraint!(ocp, :state, lb = [0, 1], ub = [1, 2], label = :cs)
+        constraint!(
+            ocp,
+            :boundary,
+            f = (r, x0, xf, v) -> (r[:] .= x0[2] + xf[2] + v[1]; nothing),
+            lb = 0,
+            ub = 1,
+            label = :cb,
+        )
+        constraint!(ocp, :control, f = (r, u, v) -> (r[:] .= u + v[2]; nothing), lb = 0, ub = 1, label = :cuu)
+        constraint!(ocp, :state, f = (r, x, v) -> (r[:] .= x + v[1:2]; nothing), lb = [0, 1], ub = [1, 2], label = :css)
+        constraint!(ocp, :mixed, f = (r, x, u, v) -> (r[:] .= x[1] + u + v[2]; nothing), lb = 1, ub = 1, label = :cm)
+        constraint!(ocp, :variable, lb = [0, 0, 0, 0], ub = [5, 5, 5, 5], label = :cv1)
+        constraint!(ocp, :variable, rg = 1:2, lb = [1, 2], ub = [3, 4], label = :cv2)
+        constraint!(ocp, :variable, rg = 3, lb = 2, ub = 3, label = :cv3)
+        constraint!(ocp, :variable, f = (r, v) -> (r[:] .= v[3]^2; nothing), lb = 0, ub = 1, label = :cv4)
+
+        # dimensions (not set yet)
+        @test dim_control_constraints(ocp) === nothing
+        @test dim_state_constraints(ocp) === nothing
+        @test dim_mixed_constraints(ocp) === nothing
+        @test dim_path_constraints(ocp) === nothing
+        @test dim_boundary_constraints(ocp) === nothing
+        @test dim_variable_constraints(ocp) === nothing
+        @test dim_control_range(ocp) === nothing
+        @test dim_state_range(ocp) === nothing
+        @test dim_variable_range(ocp) === nothing
+
+        (ξl, ξ!, ξu),
+        (ηl, η!, ηu),
+        (ψl, ψ!, ψu),
+        (ϕl, ϕ!, ϕu),
+        (θl, θ!, θu),
+        (ul, uind, uu),
+        (xl, xind, xu),
+        (vl, vind, vu) = nlp_constraints!(ocp)
+
+        v = [1, 2, 3, 4]
+
+        # control
+        @test sort(ξl) == sort([0])
+        @test sort(ξu) == sort([1])
+        r = [0.]
+        ξ!(r, -1, 1, v)
+        @test sort(r) == sort([1 + v[2]])
+
+        # state
+        @test sort(ηl) == sort([0, 1])
+        @test sort(ηu) == sort([1, 2])
+        r = [0., 0.]
+        η!(r, -1, [1, 1], v)
+        @test sort(r) == sort([1, 1] + v[1:2])
+
+        # mixed
+        @test sort(ψl) == sort([1])
+        @test sort(ψu) == sort([1])
+        r = [0.]
+        ψ!(r, -1, [1, 1], 2, v)
+        @test sort(r) == sort([3 + v[2]])
+
+        # boundary
+        @test sort(ϕl) == sort([10, 1, 0])
+        @test sort(ϕu) == sort([10, 1, 1])
+        r = [0., 0., 0.]
+        ϕ!(r, [1, 3], [4, 100], v)
+        @test sort(r) == sort([3, 4, 103 + v[1]])
+
+        # variable
+        @test sort(θl) == sort([0])
+        @test sort(θu) == sort([1])
+        r = [0.]
+        θ!(r, v)
+        @test sort(r) == sort([v[3]^2])
+
+        # box constraint
+        @test sort(ul) == sort([0])
+        @test sort(uind) == sort([1])
+        @test sort(uu) == sort([1])
+        @test sort(xl) == sort([0, 1])
+        @test sort(xind) == sort([1, 2])
+        @test sort(xu) == sort([1, 2])
+        @test sort(vl) == sort([0, 0, 0, 0, 1, 2, 2])
+        @test sort(vind) == sort([1, 2, 3, 4, 1, 2, 3])
+        @test sort(vu) == sort([5, 5, 5, 5, 3, 4, 3])
+
+        # dimensions
+        @test dim_control_constraints(ocp) == 1
+        @test dim_state_constraints(ocp) == 2
+        @test dim_mixed_constraints(ocp) == 1
+        @test dim_path_constraints(ocp) == 4
+        @test dim_boundary_constraints(ocp) == 3
+        @test dim_variable_constraints(ocp) == 1
+        @test dim_control_range(ocp) == 1
+        @test dim_state_range(ocp) == 2
+        @test dim_variable_range(ocp) == 7
+    end
+
     @testset "val vs lb and ub, errors" begin
         ocp = Model(variable = true)
 
@@ -1835,6 +2073,223 @@ function test_model() # 30 55 185
         @test dim_variable_range(ocp) == 7
     end
 
+
+    @testset "val vs lb and ub, 1/2 (in place)" begin
+        ocp = Model(variable = true, in_place = true)
+
+        time!(ocp; t0 = 0, tf = 1)
+        variable!(ocp, 4)
+        state!(ocp, 2)
+        control!(ocp, 1)
+
+        constraint!(ocp, :initial; rg = Index(2), lb = 10, ub = 10, label = :ci)
+        constraint!(ocp, :final; rg = Index(1), lb = 1, ub = 1, label = :cf)
+        constraint!(ocp, :control; lb = 0, ub = 0, label = :cu)
+        constraint!(ocp, :state; lb = [0, 1], ub = [0, 1], label = :cs)
+        constraint!(
+            ocp,
+            :boundary;
+            f = (r, x0, xf, v) -> (r[:] .= x0[2] + xf[2] + v[1]; nothing),
+            lb = 2,
+            ub = 2,
+            label = :cb,
+        )
+        constraint!(ocp, :control; f = (r, u, v) -> (r[:] .= u + v[2]; nothing), lb = 20, ub = 20, label = :cuu)
+        constraint!(
+            ocp,
+            :state;
+            f = (r, x, v) -> (r[:] .= x + v[1:2]; nothing),
+            lb = [100, 101],
+            ub = [100, 101],
+            label = :css,
+        )
+        constraint!(ocp, :mixed; f = (r, x, u, v) -> (r[:] .= x[1] + u + v[2]; nothing), lb = -1, ub = -1, label = :cm)
+        constraint!(ocp, :variable; lb = [5, 5, 5, 5], ub = [5, 5, 5, 5], label = :cv1)
+        constraint!(ocp, :variable; rg = 1:2, lb = [10, 20], ub = [10, 20], label = :cv2)
+        constraint!(ocp, :variable; rg = Index(3), lb = 1000, ub = 1000, label = :cv3)
+        constraint!(ocp, :variable; f = (r, v) -> (r[:] .= v[3]^2; nothing), lb = -10, ub = -10, label = :cv4)
+
+        # dimensions (not set yet)
+        @test dim_control_constraints(ocp) === nothing
+        @test dim_state_constraints(ocp) === nothing
+        @test dim_mixed_constraints(ocp) === nothing
+        @test dim_path_constraints(ocp) === nothing
+        @test dim_boundary_constraints(ocp) === nothing
+        @test dim_variable_constraints(ocp) === nothing
+        @test dim_control_range(ocp) === nothing
+        @test dim_state_range(ocp) === nothing
+        @test dim_variable_range(ocp) === nothing
+
+        (ξl, ξ!, ξu),
+        (ηl, η!, ηu),
+        (ψl, ψ!, ψu),
+        (ϕl, ϕ!, ϕu),
+        (θl, θ!, θu),
+        (ul, uind, uu),
+        (xl, xind, xu),
+        (vl, vind, vu) = nlp_constraints!(ocp)
+
+        v = [1, 2, 3, 4]
+
+        # control
+        @test sort(ξl) == sort([20])
+        @test sort(ξu) == sort([20])
+        r = [0.]
+        ξ!(r, -1, 1, v)
+        @test sort(r) == sort([1 + v[2]])
+
+        # state
+        @test sort(ηl) == sort([100, 101])
+        @test sort(ηu) == sort([100, 101])
+        r = [0., 0.]
+        η!(r, -1, [1, 1], v)
+        @test sort(r) == sort([1, 1] + v[1:2])
+
+        # mixed
+        @test sort(ψl) == sort([-1])
+        @test sort(ψu) == sort([-1])
+        r = [0.]
+        ψ!(r, -1, [1, 1], 2, v)
+        @test sort(r) == sort([3 + v[2]])
+
+        # boundary
+        @test sort(ϕl) == sort([10, 1, 2])
+        @test sort(ϕu) == sort([10, 1, 2])
+        r = [0., 0., 0.]
+        ϕ!(r, [1, 3], [4, 100], v)
+        @test sort(r) == sort([3, 4, 103 + v[1]])
+
+        # variable
+        @test sort(θl) == sort([-10])
+        @test sort(θu) == sort([-10])
+        r = [0.]
+        θ!(r, v)
+        @test sort(r) == sort([v[3]^2])
+
+        # box constraint
+        @test sort(ul) == sort([0])
+        @test sort(uind) == sort([1])
+        @test sort(uu) == sort([0])
+        @test sort(xl) == sort([0, 1])
+        @test sort(xind) == sort([1, 2])
+        @test sort(xu) == sort([0, 1])
+        @test sort(vl) == sort([5, 5, 5, 5, 10, 20, 1000])
+        @test sort(vind) == sort([1, 2, 3, 4, 1, 2, 3])
+        @test sort(vu) == sort([5, 5, 5, 5, 10, 20, 1000])
+
+        # dimensions
+        @test dim_control_constraints(ocp) == 1
+        @test dim_state_constraints(ocp) == 2
+        @test dim_mixed_constraints(ocp) == 1
+        @test dim_path_constraints(ocp) == 4
+        @test dim_boundary_constraints(ocp) == 3
+        @test dim_variable_constraints(ocp) == 1
+        @test dim_control_range(ocp) == 1
+        @test dim_state_range(ocp) == 2
+        @test dim_variable_range(ocp) == 7
+    end
+
+    @testset "val vs lb and ub, 2/2 (in place)" begin
+        ocp = Model(variable = true, in_place = true)
+
+        time!(ocp; t0 = 0, tf = 1)
+        variable!(ocp, 4)
+        state!(ocp, 2)
+        control!(ocp, 1)
+
+        constraint!(ocp, :initial; rg = Index(2), val = 10, label = :ci)
+        constraint!(ocp, :final; rg = Index(1), val = 1, label = :cf)
+        constraint!(ocp, :control; val = 0, label = :cu)
+        constraint!(ocp, :state; val = [0, 1], label = :cs)
+        constraint!(ocp, :boundary; f = (r, x0, xf, v) -> (r[:] .= x0[2] + xf[2] + v[1]; nothing), val = 2, label = :cb)
+        constraint!(ocp, :control; f = (r, u, v) -> (r[:] .= u + v[2]; nothing), val = 20, label = :cuu)
+        constraint!(ocp, :state; f = (r, x, v) -> (r[:] .= x + v[1:2]; nothing), val = [100, 101], label = :css)
+        constraint!(ocp, :mixed; f = (r, x, u, v) -> (r[:] .= x[1] + u + v[2]; nothing), val = -1, label = :cm)
+        constraint!(ocp, :variable; val = [5, 5, 5, 5], label = :cv1)
+        constraint!(ocp, :variable; rg = 1:2, val = [10, 20], label = :cv2)
+        constraint!(ocp, :variable; rg = Index(3), val = 1000, label = :cv3)
+        constraint!(ocp, :variable; f = (r, v) -> (r[:] .= v[3]^2; nothing), val = -10, label = :cv4)
+
+        # dimensions (not set yet)
+        @test dim_control_constraints(ocp) === nothing
+        @test dim_state_constraints(ocp) === nothing
+        @test dim_mixed_constraints(ocp) === nothing
+        @test dim_path_constraints(ocp) === nothing
+        @test dim_boundary_constraints(ocp) === nothing
+        @test dim_variable_constraints(ocp) === nothing
+        @test dim_control_range(ocp) === nothing
+        @test dim_state_range(ocp) === nothing
+        @test dim_variable_range(ocp) === nothing
+
+        (ξl, ξ!, ξu),
+        (ηl, η!, ηu),
+        (ψl, ψ!, ψu),
+        (ϕl, ϕ!, ϕu),
+        (θl, θ!, θu),
+        (ul, uind, uu),
+        (xl, xind, xu),
+        (vl, vind, vu) = nlp_constraints!(ocp)
+
+        v = [1, 2, 3, 4]
+
+        # control
+        @test sort(ξl) == sort([20])
+        @test sort(ξu) == sort([20])
+        r = [0.]
+        ξ!(r, -1, 1, v)
+        @test sort(r) == sort([1 + v[2]])
+
+        # state
+        @test sort(ηl) == sort([100, 101])
+        @test sort(ηu) == sort([100, 101])
+        r = [0., 0.]
+        η!(r, -1, [1, 1], v)
+        @test sort(r) == sort([1, 1] + v[1:2])
+
+        # mixed
+        @test sort(ψl) == sort([-1])
+        @test sort(ψu) == sort([-1])
+        r = [0.]
+        ψ!(r, -1, [1, 1], 2, v)
+        @test sort(r) == sort([3 + v[2]])
+
+        # boundary
+        @test sort(ϕl) == sort([10, 1, 2])
+        @test sort(ϕu) == sort([10, 1, 2])
+        r = [0., 0., 0.]
+        ϕ!(r, [1, 3], [4, 100], v)
+        @test sort(r) == sort([3, 4, 103 + v[1]])
+
+        # variable
+        @test sort(θl) == sort([-10])
+        @test sort(θu) == sort([-10])
+        r = [0.]
+        θ!(r, v)
+        @test sort(r) == sort([v[3]^2])
+
+        # box constraint
+        @test sort(ul) == sort([0])
+        @test sort(uind) == sort([1])
+        @test sort(uu) == sort([0])
+        @test sort(xl) == sort([0, 1])
+        @test sort(xind) == sort([1, 2])
+        @test sort(xu) == sort([0, 1])
+        @test sort(vl) == sort([5, 5, 5, 5, 10, 20, 1000])
+        @test sort(vind) == sort([1, 2, 3, 4, 1, 2, 3])
+        @test sort(vu) == sort([5, 5, 5, 5, 10, 20, 1000])
+
+        # dimensions
+        @test dim_control_constraints(ocp) == 1
+        @test dim_state_constraints(ocp) == 2
+        @test dim_mixed_constraints(ocp) == 1
+        @test dim_path_constraints(ocp) == 4
+        @test dim_boundary_constraints(ocp) == 3
+        @test dim_variable_constraints(ocp) == 1
+        @test dim_control_range(ocp) == 1
+        @test dim_state_range(ocp) == 2
+        @test dim_variable_range(ocp) == 7
+    end
+
     @testset "objective!" begin
         ocp = Model()
         @test_throws UnauthorizedCall objective!(ocp, :lagrange, (x, u) -> 0.5u^2)
diff --git a/test/test_onepass.jl b/test/test_onepass.jl
index f90bcbb4..5060a1b6 100644
--- a/test/test_onepass.jl
+++ b/test/test_onepass.jl
@@ -31,7 +31,7 @@ function test_onepass()
         d = 4
         x = 10
         u = 20
-        @test dynamics(o)(x, u) == x + u + b + 3 + d
+        @test __dynamics(o)(x, u) == x + u + b + 3 + d
     end
 
     @testset "log" begin
@@ -116,10 +116,10 @@ function test_onepass()
         x0 = 2 * x
         xf = 3 * x
         u = 3
-        @test constraint(o, :eq1)(x0, xf) == x0[1]
-        @test constraint(o, Symbol("♡"))(x0, xf) == x0[2]
-        @test dynamics(o)(x, u) == [x[2], (x[1] + 2x[2])^2]
-        @test lagrange(o)(x, u) == u^2 + x[1]
+        @test __constraint(o, :eq1)(x0, xf) == x0[1]
+        @test __constraint(o, Symbol("♡"))(x0, xf) == x0[2]
+        @test __dynamics(o)(x, u) == [x[2], (x[1] + 2x[2])^2]
+        @test __lagrange(o)(x, u) == u^2 + x[1]
 
         @def o begin
             t ∈ [0, 1], time
@@ -135,10 +135,10 @@ function test_onepass()
         x0 = 2 * x
         xf = 3 * x
         u = 3
-        @test constraint(o, :eq1)(x0, xf) == x0[1]
-        @test constraint(o, Symbol("♡"))(x0, xf) == x0[2]
-        @test dynamics(o)(x, u) == [x[2], (x[1] + 2x[2])^2]
-        @test lagrange(o)(x, u) == u^2 + x[1]
+        @test __constraint(o, :eq1)(x0, xf) == x0[1]
+        @test __constraint(o, Symbol("♡"))(x0, xf) == x0[2]
+        @test __dynamics(o)(x, u) == [x[2], (x[1] + 2x[2])^2]
+        @test __lagrange(o)(x, u) == u^2 + x[1]
 
         @def o begin
             t ∈ [0, 1], time
@@ -156,10 +156,10 @@ function test_onepass()
         xf = 3 * x
         u = 3
         c = [u, 0]
-        @test constraint(o, :eq1)(x0, xf) == x0[1]
-        @test constraint(o, Symbol("♡"))(x0, xf) == x0[2]
-        @test dynamics(o)(x, c) == [x[2], (x[1] + 2x[2])^2]
-        @test lagrange(o)(x, c) == u^2 + x[1]
+        @test __constraint(o, :eq1)(x0, xf) == x0[1]
+        @test __constraint(o, Symbol("♡"))(x0, xf) == x0[2]
+        @test __dynamics(o)(x, c) == [x[2], (x[1] + 2x[2])^2]
+        @test __lagrange(o)(x, c) == u^2 + x[1]
 
         @def o begin
             t ∈ [0, 1], time
@@ -171,7 +171,7 @@ function test_onepass()
         @test control_dimension(o) == 2
         x = [1, 2, 3]
         u = [-1, 2]
-        @test dynamics(o)(x, u) == [x[1] + 2u[2], 2x[3], x[1] + u[2]]
+        @test __dynamics(o)(x, u) == [x[1] + 2u[2], 2x[3], x[1] + u[2]]
 
         t0 = 0.0
         tf = 0.1
@@ -542,7 +542,7 @@ function test_onepass()
         @test control_dimension(o) == 2
         x = [1, 2, 3]
         u = [-1, 2]
-        @test dynamics(o)(x, u) == [x[1] + 2u[2], 2x[3], x[1] + u[2]]
+        @test __dynamics(o)(x, u) == [x[1] + 2u[2], 2x[3], x[1] + u[2]]
         @def o begin
             z ∈ R², variable
             s ∈ [0, z₁], time
@@ -556,7 +556,7 @@ function test_onepass()
         z = [5, 6]
         y = [1, 2, 3, 4]
         w = 9
-        @test dynamics(o)(y, w, z) == [y[1], y[3]^2 + w + z[1], 0, 0]
+        @test __dynamics(o)(y, w, z) == [y[1], y[3]^2 + w + z[1], 0, 0]
 
         @def o begin
             z ∈ R², variable
@@ -571,7 +571,7 @@ function test_onepass()
         z = [5, 6]
         y = [1, 2, 3, 4]
         w = 9
-        @test_throws MethodError dynamics(o)(y, w, z)
+        @test_throws MethodError __dynamics(o)(y, w, z)
 
         @def o begin
             z ∈ R², variable
@@ -588,7 +588,7 @@ function test_onepass()
         y0 = y
         yf = 3y0
         ww = 19
-        @test dynamics(o)(y, ww, z) == [y[1] + ww^2 + y[4]^3 + z[2], y[3]^2, 0, 0]
+        @test __dynamics(o)(y, ww, z) == [y[1] + ww^2 + y[4]^3 + z[2], y[3]^2, 0, 0]
 
         @def o begin
             z ∈ R², variable
@@ -606,8 +606,8 @@ function test_onepass()
         y0 = y
         yf = 3y0
         w = 11
-        @test_throws MethodError dynamics(o)(y, w, z)
-        @test mayer(o)(y0, yf, z) == y0[1] + y0[4]^3 + z[2] + yf[2]
+        @test_throws MethodError __dynamics(o)(y, w, z)
+        @test __mayer(o)(y0, yf, z) == y0[1] + y0[4]^3 + z[2] + yf[2]
     end
 
     # ---------------------------------------------------------------
@@ -628,7 +628,7 @@ function test_onepass()
         tf = 2
         x0 = [1, 2]
         xf = [3, 4]
-        @test constraint(o, :eq1)(x0, xf, tf) == x0[1] + (xf[1] + 2xf[2]^3) - tf^2
+        @test __constraint(o, :eq1)(x0, xf, tf) == x0[1] + (xf[1] + 2xf[2]^3) - tf^2
 
         n = 11
         m = 6
@@ -647,22 +647,22 @@ function test_onepass()
             [0, 0] ≤ u[1:2](t) ≤ [1, 1], (7)
             [0, 0] ≤ u[1:2:4](t) ≤ [1, 1], (8)
             0 ≤ u₂(t)^2 ≤ 1, (9)
-            u₁(t) * x[1:2](t) == 1, (10)
-            0 ≤ u₁(t) * x[1:2](t) .^ 3 ≤ 1, (11)
+            u₁(t) * x[1:2](t) == [1, 1], (10)
+            [0, 0] ≤ u₁(t) * x[1:2](t) .^ 3 ≤ [1, 1], (11)
         end
         x = Vector{Float64}(1:n)
         u = 2 * Vector{Float64}(1:m)
-        @test constraint(o, :eq1)(x) == x[1]
-        @test constraint(o, :eq2)(x) == x
-        @test constraint(o, :eq3)(x) == x[1:2]
-        @test constraint(o, :eq4)(x) == x[1:2:4]
-        @test constraint(o, :eq5)(x) == x[2]^2
-        @test constraint(o, :eq6)(u) == u
-        @test constraint(o, :eq7)(u) == u[1:2]
-        @test constraint(o, :eq8)(u) == u[1:2:4]
-        @test constraint(o, :eq9)(u) == u[2]^2
-        @test constraint(o, :eq10)(x, u) == u[1] * x[1:2]
-        @test constraint(o, :eq11)(x, u) == u[1] * x[1:2] .^ 3
+        @test __constraint(o, :eq1)(x) == x[1]
+        @test __constraint(o, :eq2)(x) == x
+        @test __constraint(o, :eq3)(x) == x[1:2]
+        @test __constraint(o, :eq4)(x) == x[1:2:4]
+        @test __constraint(o, :eq5)(x) == x[2]^2
+        @test __constraint(o, :eq6)(u) == u
+        @test __constraint(o, :eq7)(u) == u[1:2]
+        @test __constraint(o, :eq8)(u) == u[1:2:4]
+        @test __constraint(o, :eq9)(u) == u[2]^2
+        @test __constraint(o, :eq10)(x, u) == u[1] * x[1:2]
+        @test __constraint(o, :eq11)(x, u) == u[1] * x[1:2] .^ 3
 
         n = 11
         m = 6
@@ -682,24 +682,24 @@ function test_onepass()
             [0, 0] ≤ u[1:2](t) ≤ [1, 1], (7)
             [0, 0] ≤ u[1:2:4](t) ≤ [1, 1], (8)
             0 ≤ u₂(t)^2 ≤ 1, (9)
-            u₁(t) * x[1:2](t) + z + f() == 1, (10)
-            0 ≤ u₁(t) * x[1:2](t) .^ 3 + z ≤ 1, (11)
+            u₁(t) * x[1:2](t) + z + f() == [1, 1], (10)
+            [0, 0] ≤ u₁(t) * x[1:2](t) .^ 3 + z ≤ [1, 1], (11)
         end
         f() = [1, 1]
         z = 3 * Vector{Float64}(1:2)
         x = Vector{Float64}(1:n)
         u = 2 * Vector{Float64}(1:m)
-        @test constraint(o, :eq1)(x, z) == x[1]
-        @test constraint(o, :eq2)(x, z) == x
-        @test constraint(o, :eq3)(x, z) == x[1:2] - [z[1], 1]
-        @test constraint(o, :eq4)(x, z) == x[1:2:4]
-        @test constraint(o, :eq5)(x, z) == x[2]^2
-        @test constraint(o, :eq6)(u, z) == u
-        @test constraint(o, :eq7)(u, z) == u[1:2]
-        @test constraint(o, :eq8)(u, z) == u[1:2:4]
-        @test constraint(o, :eq9)(u, z) == u[2]^2
-        @test constraint(o, :eq10)(x, u, z) == u[1] * x[1:2] + z + f()
-        @test constraint(o, :eq11)(x, u, z) == u[1] * x[1:2] .^ 3 + z
+        @test __constraint(o, :eq1)(x, z) == x[1]
+        @test __constraint(o, :eq2)(x, z) == x
+        @test __constraint(o, :eq3)(x, z) == x[1:2] - [z[1], 1]
+        @test __constraint(o, :eq4)(x, z) == x[1:2:4]
+        @test __constraint(o, :eq5)(x, z) == x[2]^2
+        @test __constraint(o, :eq6)(u, z) == u
+        @test __constraint(o, :eq7)(u, z) == u[1:2]
+        @test __constraint(o, :eq8)(u, z) == u[1:2:4]
+        @test __constraint(o, :eq9)(u, z) == u[2]^2
+        @test __constraint(o, :eq10)(x, u, z) == u[1] * x[1:2] + z + f()
+        @test __constraint(o, :eq11)(x, u, z) == u[1] * x[1:2] .^ 3 + z
 
         @def o begin
             t ∈ [0, 1], time
@@ -719,10 +719,10 @@ function test_onepass()
         x0 = 2 * x
         xf = 3 * x
         u = 3
-        @test constraint(o, :eq1)(x0, xf) == x0[1]
-        @test constraint(o, Symbol("♡"))(x0, xf) == x0[2]
-        @test dynamics(o)(x, u) == [x[2], (x[1] + 2x[2])^2]
-        @test lagrange(o)(x, u) == u^2 + x[1]
+        @test __constraint(o, :eq1)(x0, xf) == x0[1]
+        @test __constraint(o, Symbol("♡"))(x0, xf) == x0[2]
+        @test __dynamics(o)(x, u) == [x[2], (x[1] + 2x[2])^2]
+        @test __lagrange(o)(x, u) == u^2 + x[1]
 
         @def o begin
             t ∈ [0, 1], time
@@ -740,10 +740,10 @@ function test_onepass()
         x0 = 2 * x
         xf = 3 * x
         u = 3
-        @test constraint(o, :eq1)(x0, xf) == x0[1]
-        @test constraint(o, Symbol("♡"))(x0, xf) == x0[2]
-        @test dynamics(o)(x, u) == [x[2], (x[1] + 2x[2])^2]
-        @test lagrange(o)(x, u) == u^2 + x[1]
+        @test __constraint(o, :eq1)(x0, xf) == x0[1]
+        @test __constraint(o, Symbol("♡"))(x0, xf) == x0[2]
+        @test __dynamics(o)(x, u) == [x[2], (x[1] + 2x[2])^2]
+        @test __lagrange(o)(x, u) == u^2 + x[1]
 
         @def o begin
             z ∈ R², variable
@@ -763,10 +763,10 @@ function test_onepass()
         xf = 3 * [1, 2]
         u = 3
         z = [4, 5]
-        @test constraint(o, :eq1)(x0, xf, z) == x0[1]
-        @test constraint(o, Symbol("♡"))(x0, xf, z) == x0[2]
-        @test dynamics(o)(x, u, z) == [x[2], (x[1] + 2x[2])^2 + z[1]]
-        @test lagrange(o)(x, u, z) == u^2 + z[2] * x[1]
+        @test __constraint(o, :eq1)(x0, xf, z) == x0[1]
+        @test __constraint(o, Symbol("♡"))(x0, xf, z) == x0[2]
+        @test __dynamics(o)(x, u, z) == [x[2], (x[1] + 2x[2])^2 + z[1]]
+        @test __lagrange(o)(x, u, z) == u^2 + z[2] * x[1]
 
         @def o begin
             t ∈ [0, 1], time
@@ -783,10 +783,10 @@ function test_onepass()
         xf = [4, 5]
         x = [1, 2]
         u = 3
-        @test constraint(o, :eq1)(x0, xf) == x0[1]^2 + xf[2]
-        @test constraint(o, Symbol("♡"))(x0, xf) == x0[2]
-        @test dynamics(o)(x, u) == [x[2], x[1]^2]
-        @test lagrange(o)(x, u) == u^2 + x[1]
+        @test __constraint(o, :eq1)(x0, xf) == x0[1]^2 + xf[2]
+        @test __constraint(o, Symbol("♡"))(x0, xf) == x0[2]
+        @test __dynamics(o)(x, u) == [x[2], x[1]^2]
+        @test __lagrange(o)(x, u) == u^2 + x[1]
 
         @def o begin
             z ∈ R, variable
@@ -805,10 +805,10 @@ function test_onepass()
         x = [1, 2]
         u = 3
         z = 4
-        @test constraint(o, :eq1)(x0, xf, z) == x0[1] - z
-        @test constraint(o, Symbol("♡"))(x0, xf, z) == x0[2]
-        @test dynamics(o)(x, u, z) == [x[2], x[1]^2 + z]
-        @test lagrange(o)(x, u, z) == u^2 + z * x[1]
+        @test __constraint(o, :eq1)(x0, xf, z) == x0[1] - z
+        @test __constraint(o, Symbol("♡"))(x0, xf, z) == x0[2]
+        @test __dynamics(o)(x, u, z) == [x[2], x[1]^2 + z]
+        @test __lagrange(o)(x, u, z) == u^2 + z * x[1]
 
         @def o begin
             z ∈ R, variable
@@ -828,11 +828,11 @@ function test_onepass()
         x = [1, 2]
         u = 3
         z = 4
-        @test constraint(o, :eq1)(x0, xf, z) == x0[1] - z
-        @test constraint(o, :eq2)(x0, xf, z) == xf[2]^2
-        @test constraint(o, Symbol("♡"))(x0, xf, z) == x0
-        @test dynamics(o)(x, u, z) == [x[2], x[1]^2 + z]
-        @test lagrange(o)(x, u, z) == u^2 + z * x[1]
+        @test __constraint(o, :eq1)(x0, xf, z) == x0[1] - z
+        @test __constraint(o, :eq2)(x0, xf, z) == xf[2]^2
+        @test __constraint(o, Symbol("♡"))(x0, xf, z) == x0
+        @test __dynamics(o)(x, u, z) == [x[2], x[1]^2 + z]
+        @test __lagrange(o)(x, u, z) == u^2 + z * x[1]
 
         @def o begin
             z ∈ R, variable
@@ -852,11 +852,11 @@ function test_onepass()
         x = [1, 2]
         u = 3
         z = 4
-        @test constraint(o, :eq1)(x0, xf, z) == x0[1] - z
-        @test constraint(o, :eq2)(x0, xf, z) == xf[2]^2
-        @test constraint(o, :eq3)(x0, xf, z) == x0
-        @test dynamics(o)(x, u, z) == [x[2], x[1]^2 + z]
-        @test lagrange(o)(x, u, z) == u^2 + z * x[1]
+        @test __constraint(o, :eq1)(x0, xf, z) == x0[1] - z
+        @test __constraint(o, :eq2)(x0, xf, z) == xf[2]^2
+        @test __constraint(o, :eq3)(x0, xf, z) == x0
+        @test __dynamics(o)(x, u, z) == [x[2], x[1]^2 + z]
+        @test __lagrange(o)(x, u, z) == u^2 + z * x[1]
         @test constraints(o)[:eq1][3] == 0
         @test constraints(o)[:eq1][4] == 1
         @test constraints(o)[:eq2][3] == 0
@@ -894,20 +894,20 @@ function test_onepass()
         u = 4
         v = [5, 6]
         z = v[1] + 2v[2]
-        @test constraint(o, :eq1)(x0, xf, v) == x0 - v[1]
-        @test constraint(o, :eq2)(x0, xf, v) == xf - v[1]
-        @test constraint(o, :eq3)(x0, xf, v) == x0 - v[1]
-        @test constraint(o, :eq4)(x0, xf, v) == xf - v[1]
-        @test constraint(o, :eq5)(x0, xf, v) == x0 + xf - v[2]
-        @test constraint(o, :eq6)(x0, xf, v) == x0 + xf - v[2]
-        @test constraint(o, :eq7)(x, v) == x - v[1]
-        @test constraint(o, :eq9)(x, v) == x - z
-        @test constraint(o, :eq10)(u, v) == u - z
-        @test constraint(o, :eq11)(x, u, v) == x + u - z
-        @test constraint(o, :eq12)(v) == v[1]
-        @test constraint(o, :eq13)(v) == v[1]
-        @test constraint(o, :eq14)(v) == v[1] + 2v[2]
-        @test constraint(o, :eq15)(v) == v[1] + 2v[2]
+        @test __constraint(o, :eq1)(x0, xf, v) == x0 - v[1]
+        @test __constraint(o, :eq2)(x0, xf, v) == xf - v[1]
+        @test __constraint(o, :eq3)(x0, xf, v) == x0 - v[1]
+        @test __constraint(o, :eq4)(x0, xf, v) == xf - v[1]
+        @test __constraint(o, :eq5)(x0, xf, v) == x0 + xf - v[2]
+        @test __constraint(o, :eq6)(x0, xf, v) == x0 + xf - v[2]
+        @test __constraint(o, :eq7)(x, v) == x - v[1]
+        @test __constraint(o, :eq9)(x, v) == x - z
+        @test __constraint(o, :eq10)(u, v) == u - z
+        @test __constraint(o, :eq11)(x, u, v) == x + u - z
+        @test __constraint(o, :eq12)(v) == v[1]
+        @test __constraint(o, :eq13)(v) == v[1]
+        @test __constraint(o, :eq14)(v) == v[1] + 2v[2]
+        @test __constraint(o, :eq15)(v) == v[1] + 2v[2]
 
         @def o begin
             v ∈ R, variable
@@ -1573,9 +1573,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == 0.5u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == 0.5u^2
         @test criterion(o) == :min
 
         t0 = 0
@@ -1601,9 +1601,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == -0.5u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == -0.5u^2
         @test criterion(o) == :min
 
         t0 = 0
@@ -1629,9 +1629,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == 0.5u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == 0.5u^2
         @test criterion(o) == :min
 
         t0 = 0
@@ -1657,9 +1657,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == 0.5u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == 0.5u^2
         @test criterion(o) == :min
 
         t0 = 0
@@ -1685,9 +1685,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == -0.5u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == -0.5u^2
         @test criterion(o) == :min
 
         t0 = 0
@@ -1713,9 +1713,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == (-0.5 + tf) * u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == (-0.5 + tf) * u^2
         @test criterion(o) == :min
 
         t0 = 0
@@ -1767,9 +1767,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == 0.5u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == 0.5u^2
         @test criterion(o) == :max
 
         t0 = 0
@@ -1795,9 +1795,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == -0.5u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == -0.5u^2
         @test criterion(o) == :max
 
         t0 = 0
@@ -1823,9 +1823,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == 0.5u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == 0.5u^2
         @test criterion(o) == :max
 
         t0 = 0
@@ -1851,9 +1851,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == 0.5u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == 0.5u^2
         @test criterion(o) == :max
 
         t0 = 0
@@ -1879,9 +1879,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == -0.5u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == -0.5u^2
         @test criterion(o) == :max
 
         # -----------------------------------
@@ -1948,8 +1948,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 3xf
-        @test lagrange(o)(x, u) == x + u
+        @test __mayer(o)(x0, xf) == x0 + 3xf
+        @test __lagrange(o)(x, u) == x + u
         @test criterion(o) == :min
 
         @def o begin
@@ -1962,8 +1962,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 2xf
-        @test lagrange(o)(x, u) == x + u
+        @test __mayer(o)(x0, xf) == x0 + 2xf
+        @test __lagrange(o)(x, u) == x + u
         @test criterion(o) == :min
 
         @def o begin
@@ -1976,8 +1976,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 2xf
-        @test lagrange(o)(x, u) == 2(x + u)
+        @test __mayer(o)(x0, xf) == x0 + 2xf
+        @test __lagrange(o)(x, u) == 2(x + u)
         @test criterion(o) == :min
 
         @def o begin
@@ -1990,8 +1990,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 2xf
-        @test lagrange(o)(x, u) == -(x + u)
+        @test __mayer(o)(x0, xf) == x0 + 2xf
+        @test __lagrange(o)(x, u) == -(x + u)
         @test criterion(o) == :min
 
         @def o begin
@@ -2004,8 +2004,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 2xf
-        @test lagrange(o)(x, u) == -2(x + u)
+        @test __mayer(o)(x0, xf) == x0 + 2xf
+        @test __lagrange(o)(x, u) == -2(x + u)
         @test criterion(o) == :min
 
         @test_throws ParsingError @def o begin
@@ -2035,8 +2035,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 5xf
-        @test lagrange(o)(x, u) == x + u
+        @test __mayer(o)(x0, xf) == x0 + 5xf
+        @test __lagrange(o)(x, u) == x + u
         @test criterion(o) == :max
 
         @def o begin
@@ -2049,8 +2049,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 2xf
-        @test lagrange(o)(x, u) == 2(x + u)
+        @test __mayer(o)(x0, xf) == x0 + 2xf
+        @test __lagrange(o)(x, u) == 2(x + u)
         @test criterion(o) == :max
 
         @def o begin
@@ -2063,8 +2063,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 2xf
-        @test lagrange(o)(x, u) == -(x + u)
+        @test __mayer(o)(x0, xf) == x0 + 2xf
+        @test __lagrange(o)(x, u) == -(x + u)
         @test criterion(o) == :max
 
         @def o begin
@@ -2077,8 +2077,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 2xf
-        @test lagrange(o)(x, u) == -2(x + u)
+        @test __mayer(o)(x0, xf) == x0 + 2xf
+        @test __lagrange(o)(x, u) == -2(x + u)
         @test criterion(o) == :max
 
         @test_throws ParsingError @def o begin
@@ -2108,8 +2108,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 2xf
-        @test lagrange(o)(x, u) == x + u
+        @test __mayer(o)(x0, xf) == x0 + 2xf
+        @test __lagrange(o)(x, u) == x + u
         @test criterion(o) == :min
 
         @def o begin
@@ -2122,8 +2122,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 2xf
-        @test lagrange(o)(x, u) == 2(x + u)
+        @test __mayer(o)(x0, xf) == x0 + 2xf
+        @test __lagrange(o)(x, u) == 2(x + u)
         @test criterion(o) == :min
 
         @def o begin
@@ -2136,8 +2136,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == -(x0 + 2xf)
-        @test lagrange(o)(x, u) == x + u
+        @test __mayer(o)(x0, xf) == -(x0 + 2xf)
+        @test __lagrange(o)(x, u) == x + u
         @test criterion(o) == :min
 
         @def o begin
@@ -2150,8 +2150,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == -(x0 + 2xf)
-        @test lagrange(o)(x, u) == 2(x + u)
+        @test __mayer(o)(x0, xf) == -(x0 + 2xf)
+        @test __lagrange(o)(x, u) == 2(x + u)
         @test criterion(o) == :min
 
         @test_throws ParsingError @def o begin
@@ -2181,8 +2181,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 2xf
-        @test lagrange(o)(x, u) == x + u
+        @test __mayer(o)(x0, xf) == x0 + 2xf
+        @test __lagrange(o)(x, u) == x + u
         @test criterion(o) == :max
 
         @def o begin
@@ -2195,8 +2195,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == x0 + 2xf
-        @test lagrange(o)(x, u) == 2(x + u)
+        @test __mayer(o)(x0, xf) == x0 + 2xf
+        @test __lagrange(o)(x, u) == 2(x + u)
         @test criterion(o) == :max
 
         @def o begin
@@ -2209,8 +2209,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == -(x0 + 2xf)
-        @test lagrange(o)(x, u) == x + u
+        @test __mayer(o)(x0, xf) == -(x0 + 2xf)
+        @test __lagrange(o)(x, u) == x + u
         @test criterion(o) == :max
 
         @def o begin
@@ -2223,8 +2223,8 @@ function test_onepass()
         u = 2
         x0 = 3
         xf = 4
-        @test mayer(o)(x0, xf) == -(x0 + 2xf)
-        @test lagrange(o)(x, u) == 2(x + u)
+        @test __mayer(o)(x0, xf) == -(x0 + 2xf)
+        @test __lagrange(o)(x, u) == 2(x + u)
         @test criterion(o) == :max
 
         @test_throws ParsingError @def o begin
@@ -2258,7 +2258,7 @@ function test_onepass()
         y0 = [1, 2, 3, 4]
         yf = 2 * [1, 2, 3, 4]
         @test is_min(o)
-        @test mayer(o)(y0, yf) == y0[3] + yf[4]
+        @test __mayer(o)(y0, yf) == y0[3] + yf[4]
 
         @def o begin
             s ∈ [0, 1], time
@@ -2271,7 +2271,7 @@ function test_onepass()
         y0 = [1, 2, 3, 4]
         yf = 2 * [1, 2, 3, 4]
         @test is_max(o)
-        @test mayer(o)(y0, yf) == y0[3] + yf[4]
+        @test __mayer(o)(y0, yf) == y0[3] + yf[4]
 
         @def o begin
             z ∈ R^2, variable
@@ -2286,7 +2286,7 @@ function test_onepass()
         y0 = [1, 2, 3, 4]
         yf = 2 * [1, 2, 3, 4]
         @test is_min(o)
-        @test mayer(o)(y0, yf, z) == y0[3] + yf[4] + z[2]
+        @test __mayer(o)(y0, yf, z) == y0[3] + yf[4] + z[2]
 
         @def o begin
             z ∈ R², variable
@@ -2304,8 +2304,8 @@ function test_onepass()
         y0 = y
         yf = 3y0
         w = 7
-        @test dynamics(o)(y, w, z) == [y[1] + w^2 + y[4]^3 + z[2], y[3]^2, 0, 0]
-        @test mayer(o)(y0, yf, z) == y0[3] + yf[4] + z[2]
+        @test __dynamics(o)(y, w, z) == [y[1] + w^2 + y[4]^3 + z[2], y[3]^2, 0, 0]
+        @test __mayer(o)(y0, yf, z) == y0[3] + yf[4] + z[2]
 
         @def o begin
             z ∈ R², variable
@@ -2323,8 +2323,8 @@ function test_onepass()
         y0 = y
         yf = 3y0
         w = 7
-        @test dynamics(o)(y, w, z) == [y[1] + w^2 + y[4]^3 + z[2], y[3]^2, 0, 0]
-        @test mayer(o)(y0, yf, z) == y0[3] + yf[4] + z[2]
+        @test __dynamics(o)(y, w, z) == [y[1] + w^2 + y[4]^3 + z[2], y[3]^2, 0, 0]
+        @test __mayer(o)(y0, yf, z) == y0[3] + yf[4] + z[2]
 
         @def o begin
             z ∈ R², variable
@@ -2339,7 +2339,7 @@ function test_onepass()
         z = [5, 6]
         y0 = y
         yf = 3y0
-        @test mayer(o)(y0, yf, z) == y0[1] + y0[4]^3 + z[2] + yf[2]
+        @test __mayer(o)(y0, yf, z) == y0[1] + y0[4]^3 + z[2] + yf[2]
 
         @def o begin
             z ∈ R², variable
@@ -2357,8 +2357,8 @@ function test_onepass()
         y0 = y
         yf = 3y0
         w = 11
-        @test dynamics(o)(y, w, z) == [y[1] + w^2 + y[4]^3 + z[2], y[3]^2, 0, 0]
-        @test_throws UndefVarError mayer(o)(y0, yf, z)
+        @test __dynamics(o)(y, w, z) == [y[1] + w^2 + y[4]^3 + z[2], y[3]^2, 0, 0]
+        @test_throws UndefVarError __mayer(o)(y0, yf, z)
     end
 
     # ---------------------------------------------------------------
@@ -2382,7 +2382,7 @@ function test_onepass()
         d = 4
         x = 10
         u = 20
-        @test dynamics(o)(x, u) == x + u + b + 3 + d
+        @test __dynamics(o)(x, u) == x + u + b + 3 + d
     end
 
     # ---------------------------------------------------------------
@@ -2486,9 +2486,9 @@ function test_onepass()
             0
             1
         ]
-        @test constraint(o, :eq1)(x0, xf) == x0
-        @test dynamics(o)(x, u) == A * x + B * u
-        @test lagrange(o)(x, u) == 0.5u^2
+        @test __constraint(o, :eq1)(x0, xf) == x0
+        @test __dynamics(o)(x, u) == A * x + B * u
+        @test __lagrange(o)(x, u) == 0.5u^2
         @test criterion(o) == :min
 
         @def o begin
@@ -2510,12 +2510,12 @@ function test_onepass()
         x = [1, 2]
         u = 3
         z = 4
-        @test constraint(o, :eq1)(x0, xf, z) == x0[1] - z
-        @test constraint(o, :eq2)(x0, xf, z) == xf[2]^2
-        @test constraint(o, Symbol("♡"))(x0, xf, z) == x0
-        @test constraint(o, :eq3)(z) == z
-        @test dynamics(o)(x, u, z) == [x[2], x[1]^2 + z]
-        @test lagrange(o)(x, u, z) == u^2 + z * x[1]
+        @test __constraint(o, :eq1)(x0, xf, z) == x0[1] - z
+        @test __constraint(o, :eq2)(x0, xf, z) == xf[2]^2
+        @test __constraint(o, Symbol("♡"))(x0, xf, z) == x0
+        @test __constraint(o, :eq3)(z) == z
+        @test __dynamics(o)(x, u, z) == [x[2], x[1]^2 + z]
+        @test __lagrange(o)(x, u, z) == u^2 + z * x[1]
 
         @def o begin
             z in R, variable
@@ -2536,12 +2536,12 @@ function test_onepass()
         x = [1, 2]
         u = 3
         z = 4
-        @test constraint(o, :eq1)(x0, xf, z) == x0[1] - z
-        @test constraint(o, :eq2)(x0, xf, z) == xf[2]^2
-        @test constraint(o, Symbol("♡"))(x0, xf, z) == x0
-        @test constraint(o, :eq3)(z) == z
-        @test dynamics(o)(x, u, z) == [x[2], x[1]^2 + z]
-        @test lagrange(o)(x, u, z) == u^2 + z * x[1]
+        @test __constraint(o, :eq1)(x0, xf, z) == x0[1] - z
+        @test __constraint(o, :eq2)(x0, xf, z) == xf[2]^2
+        @test __constraint(o, Symbol("♡"))(x0, xf, z) == x0
+        @test __constraint(o, :eq3)(z) == z
+        @test __dynamics(o)(x, u, z) == [x[2], x[1]^2 + z]
+        @test __lagrange(o)(x, u, z) == u^2 + z * x[1]
 
         @def o begin
             z in R^2, variable
@@ -2564,11 +2564,11 @@ function test_onepass()
         x = [1, 2]
         u = [3, 0]
         z = [4, 1]
-        @test constraint(o, :eq1)(x0, xf, z) == x0[1] - z[1]
-        @test constraint(o, :eq2)(x0, xf, z) == xf[2]^2
-        @test constraint(o, Symbol("♡"))(x0, xf, z) == x0
-        @test constraint(o, :eq3)(z) == z[1]
-        @test dynamics(o)(x, u, z) == [x[2], x[1]^2 + z[1]]
-        @test lagrange(o)(x, u, z) == u[1]^2 + z[1] * x[1]
+        @test __constraint(o, :eq1)(x0, xf, z) == x0[1] - z[1]
+        @test __constraint(o, :eq2)(x0, xf, z) == xf[2]^2
+        @test __constraint(o, Symbol("♡"))(x0, xf, z) == x0
+        @test __constraint(o, :eq3)(z) == z[1]
+        @test __dynamics(o)(x, u, z) == [x[2], x[1]^2 + z[1]]
+        @test __lagrange(o)(x, u, z) == u[1]^2 + z[1] * x[1]
     end
 end
diff --git a/test/test_plot.jl b/test/test_plot.jl
index 7e6fc3ac..2652dd18 100644
--- a/test/test_plot.jl
+++ b/test/test_plot.jl
@@ -51,4 +51,4 @@ function test_plot()
     @test plot(sol, layout = :split) isa Plots.Plot
     @test plot(sol, layout = :group) isa Plots.Plot
     @test display(sol) isa Nothing
-end
+end
\ No newline at end of file
diff --git a/test/test_plot_joseph.jl b/test/test_plot_joseph.jl
deleted file mode 100644
index 0f401559..00000000
--- a/test/test_plot_joseph.jl
+++ /dev/null
@@ -1,54 +0,0 @@
-# using OptimalControl
-
-# @def ocp begin
-
-#     t ∈ [ 0, 1 ], time
-#     x ∈ R, state
-#     u ∈ R, control
-#     x(0) == 0
-#     -1 ≤ u(t) ≤ 1
-#     ẋ(t) == -x(t) + u(t)
-#     x(1) → min
-
-# end
-
-# sol = solve(ocp, grid_size=20, print_level=5)
-# plt = plot(sol)
-
-# a la main
-using CTBase
-using PlotUtils
-
-n = 1
-m = 1
-t0 = 0
-tf = 1
-x0 = 0
-x = t -> -1 + 1e-8 * rand()
-p = t -> 0 + 1e-8 * rand()
-u = t -> 0
-objective = 1
-#
-N = 201
-times = range(t0, tf, N)
-#
-
-sol = OptimalControlSolution()
-sol.state_dimension = n
-sol.control_dimension = m
-sol.time_grid = times
-sol.time_name = "t"
-sol.state = x
-sol.state_name = "x"
-sol.state_components_names = ["x"]
-sol.costate = p
-sol.control = u
-sol.control_name = "u"
-sol.control_components_names = ["u"]
-sol.objective = objective
-sol.iterations = 0
-sol.stopping = :dummy
-sol.message = "ceci est un test"
-sol.success = true
-
-plt = plot(sol)
diff --git a/test/test_plot_manual.jl b/test/test_plot_manual.jl
deleted file mode 100644
index a5ce33a0..00000000
--- a/test/test_plot_manual.jl
+++ /dev/null
@@ -1,211 +0,0 @@
-using CTBase
-using Plots
-
-layout = :split
-size = (900, 600)
-control_plt = :all
-
-#
-do_plot_2 = true
-
-# ----------------------------------------
-# SOL 1
-n = 2
-m = 1
-t0 = 0.0
-tf = 1.0
-x0 = [-1.0, 0.0]
-xf = [0.0, 0.0]
-a = x0[1]
-b = x0[2]
-C = [
-    -(tf - t0)^3/6.0 (tf - t0)^2/2.0
-    -(tf - t0)^2/2.0 (tf-t0)
-]
-D = [-a - b * (tf - t0), -b] + xf
-p0 = C \ D
-α = p0[1]
-β = p0[2]
-x =
-    t -> [
-        a + b * (t - t0) + β * (t - t0)^2 / 2.0 - α * (t - t0)^3 / 6.0,
-        b + β * (t - t0) - α * (t - t0)^2 / 2.0,
-    ]
-p = t -> [α, -α * (t - t0) + β]
-u = t -> [p(t)[2]]
-objective = 0.5 * (α^2 * (tf - t0)^3 / 3 + β^2 * (tf - t0) - α * β * (tf - t0)^2)
-#
-N = 201
-times = range(t0, tf, N)
-#
-
-sol = OptimalControlSolution()
-sol.state_dimension = n
-sol.control_dimension = m
-sol.time_grid = times
-sol.time_name = "t"
-sol.state = x
-sol.state_name = "x"
-sol.state_components_names = ["x" * ctindices(i) for i ∈ range(1, n)]
-sol.costate = p
-sol.control = u
-sol.control_name = "u"
-sol.control_components_names = ["u"]
-sol.objective = objective
-sol.iterations = 0
-sol.stopping = :dummy
-sol.message = "ceci est un test"
-sol.success = true
-
-#
-plt = plot(
-    sol,
-    layout = layout,
-    control = control_plt,
-    size = size,
-    flip = true,
-    linewidth = 5,
-    solution_label = "sol1",
-)
-#plot(sol, layout=:group)
-#ps=plot(sol, :time, (:state, 1))
-#plot!(ps, sol, :time, (:control, 1))
-
-# ----------------------------------------
-# SOL 2
-n = 2
-m = 1
-t0 = 0.0
-tf = 1.0
-x0 = [-1.0, -1.0]
-xf = [0.0, 0.0]
-a = x0[1]
-b = x0[2]
-C = [
-    -(tf - t0)^3/6.0 (tf - t0)^2/2.0
-    -(tf - t0)^2/2.0 (tf-t0)
-]
-D = [-a - b * (tf - t0), -b] + xf
-p0 = C \ D
-α = p0[1]
-β = p0[2]
-x =
-    t -> [
-        a + b * (t - t0) + β * (t - t0)^2 / 2.0 - α * (t - t0)^3 / 6.0,
-        b + β * (t - t0) - α * (t - t0)^2 / 2.0,
-    ]
-p = t -> [α, -α * (t - t0) + β]
-u = t -> [p(t)[2]]
-objective = 0.5 * (α^2 * (tf - t0)^3 / 3 + β^2 * (tf - t0) - α * β * (tf - t0)^2)
-#
-N = 201
-times = range(t0, tf, N)
-#
-
-sol = OptimalControlSolution()
-sol.state_dimension = n
-sol.control_dimension = m
-sol.time_grid = times
-sol.time_name = "t"
-sol.state = x
-sol.state_name = "y"
-sol.state_components_names = ["y" * ctindices(i) for i ∈ range(1, n)]
-sol.costate = p
-sol.control = u
-sol.control_name = "v"
-sol.control_components_names = ["v"]
-sol.objective = objective
-sol.iterations = 0
-sol.stopping = :dummy
-sol.message = "ceci est un test"
-sol.success = true
-
-if do_plot_2
-    plot!(plt, sol, layout = layout, size = size, control = control_plt, solution_label = "sol2")
-else
-    plt
-end
-
-# ----------------------------------------
-# Orbital transfer consumption
-#=
-n=4
-m=2
-
-x0 = [-42272.67, 0, 0, -5796.72]
-μ      = 5.1658620912*1e12
-rf     = 42165.0 ;
-rf3    = rf^3  ;
-m0     = 2000.0
-F_max = 100.0
-γ_max  = F_max*3600.0^2/(m0*10^3)
-t0     = 0.0
-α      = sqrt(μ/rf3);
-β      = 0.0
-
-tol    = 1e-9;
-
-F_max_100  = 100.0
-
-tf_min = 13.40318195708344 # minimal time for Fmax = 100 N
-tf = 1.5*tf_min
-
-Th(F) = F*3600.0^2/(10^3)
-u_max = Th(F_max)
-
-# the solution
-x0 = [x0; 0]
-
-u0(x,p) = [0, 0]
-u1(x,p) = p[3:4]/norm(p[3:4])
-
-Hc(x,p) = p[1]*x[3] + p[2]*x[4] + p[3]*(-μ*x[1]/norm(x[1:2])^3) + p[4]*(-μ*x[2]/norm(x[1:2])^3)
-H(x,p,u) = -norm(u) + Hc(x,p) + u[1]*p[3]*γ_max + u[2]*p[4]*γ_max + p[5]*norm(u)
-H0(x,p) = H(x,p,u0(x,p)) 
-H1(x,p) = H(x,p,u1(x,p))
-
-# Flow
-f0 = Flow(Hamiltonian(H0));
-f1 = Flow(Hamiltonian(H1));
-
-# Initial guess
-p0_guess = [0.02698412111231433, 0.006910835140705538, 0.050397371862031096, -0.0032972040120747836, -1.0076835239866583e-23]
-ti_guess = [0.4556797711668658, 3.6289692721936913, 11.683607683450061, 12.505465498856514]
-ξ_guess  = [p0_guess;ti_guess]
-
-p0 = p0_guess
-t1, t2, t3, t4 = ti_guess
-
-# computing x, p, u
-f = f1 * (t1, f0) * (t2, f1) * (t3, f0) * (t4, f1)
-ode_sol  = f((t0, tf), x0, p0)
-
-x = t -> ode_sol(t)[1:4]
-p = t -> ode_sol(t)[6:9]
-u = t -> [0,0]*(t ∈ Interval(t1,t2)∪Interval(t3,t4)) +
-         p(t)[3:4]/norm(p(t)[3:4])*(t ∈ Interval(t0,t1)∪Interval(t2,t3)∪Interval(t4,tf))
-objective = ode_sol(tf)[5]
-
-#
-N=201
-times = range(t0, tf, N)
-#
-sol = OptimalControlSolution() #n, m, times, x, p, u)
-sol.state_dimension = n
-sol.control_dimension = m
-sol.time_grid = times
-sol.state = x
-sol.state_name = "y"
-sol.state_components_names = [ "x" * ctindices(1), "x" * ctindices(2), "v" * ctindices(1), "v" * ctindices(2)]
-sol.adjoint = p
-sol.control = u
-sol.control_name = "u"
-sol.control_components_names = [ "u" * ctindices(i) for i ∈ range(1, m)]
-sol.objective = objective
-sol.iterations = 0
-sol.message = "structure: B+B0B+B0B+"
-sol.success = true
-sol.infos[:resolution] = :numerical
-
-plt_transfert = plot(sol, layout=:split, size=(900, 600))
-=#
diff --git a/test/test_utils.jl b/test/test_utils.jl
index a2dc236f..2d34ffa8 100644
--- a/test/test_utils.jl
+++ b/test/test_utils.jl
@@ -78,4 +78,29 @@ function test_utils()
         @test CTBase.ctupperscripts(019) == "¹⁹"
         @test CTBase.ctupperscripts(109) == "¹⁰⁹"
     end
-end
+
+    @testset "to_out_of_place" begin
+
+    function f1!(r, x)
+        r[1] = x
+        r[2] = x + 1
+        return nothing
+    end
+    @test CTBase.to_out_of_place(f1!, 2)(1.0) == [1.0, 2.0]
+
+    function f2!(r, x, y)
+        r[:] .= x + y
+        return nothing
+    end
+    @test CTBase.to_out_of_place(f2!, 1; T = Int32)(1, 2) == 3 
+
+    function f3!(r, x; y = 1)
+        r[:] .= x + y
+        return nothing
+    end
+    @test CTBase.to_out_of_place(f3!, 1; T = Int32)(1; y = 2) == 3 
+
+    @test isnothing( CTBase.to_out_of_place(nothing, 1) )
+
+    end
+end
\ No newline at end of file