diff --git a/docs/make.jl b/docs/make.jl
index 3e732fc..e4f39eb 100644
--- a/docs/make.jl
+++ b/docs/make.jl
@@ -2,16 +2,16 @@ using Documenter
 using DocumenterMermaid
 using CTFlows
 
-makedocs(
-    sitename = "CTFlows.jl",
-    format = Documenter.HTML(
-        prettyurls = false,
-        assets = [
+makedocs(;
+    sitename="CTFlows.jl",
+    format=Documenter.HTML(;
+        prettyurls=false,
+        assets=[
             asset("https://control-toolbox.org/assets/css/documentation.css"),
             asset("https://control-toolbox.org/assets/js/documentation.js"),
         ],
     ),
-    pages = ["Introduction" => "index.md", "API" => "api.md", "Developers" => "dev.md"],
+    pages=["Introduction" => "index.md", "API" => "api.md", "Developers" => "dev.md"],
 )
 
-deploydocs(repo = "github.com/control-toolbox/CTFlows.jl.git", devbranch = "main")
+deploydocs(; repo="github.com/control-toolbox/CTFlows.jl.git", devbranch="main")
diff --git a/ext/CTFlowsODE.jl b/ext/CTFlowsODE.jl
index 03fd1e1..7ebbf1f 100644
--- a/ext/CTFlowsODE.jl
+++ b/ext/CTFlowsODE.jl
@@ -17,7 +17,7 @@ const DCoTangent = ctVector
 const ctgradient = CTBase.__ctgradient
 
 # types
-abstract type AbstractFlow{D, U} end
+abstract type AbstractFlow{D,U} end
 
 # from CTFlows
 const __create_hamiltonian = CTFlows.__create_hamiltonian
diff --git a/ext/concatenation.jl b/ext/concatenation.jl
index f56a308..058d60a 100644
--- a/ext/concatenation.jl
+++ b/ext/concatenation.jl
@@ -1,12 +1,15 @@
-function __concat_rhs(F::AbstractFlow{D, U}, G::AbstractFlow{D, U}, t_switch::Time) where {D, U}
+function __concat_rhs(
+    F::AbstractFlow{D,U}, G::AbstractFlow{D,U}, t_switch::Time
+) where {D,U}
     function rhs!(du::D, u::U, p::Variable, t::Time)
-        t < t_switch ? F.rhs!(du, u, p, t) : G.rhs!(du, u, p, t)
+        return t < t_switch ? F.rhs!(du, u, p, t) : G.rhs!(du, u, p, t)
     end
     return rhs!
 end
 
 function __concat_rhs(F::VectorFieldFlow, G::VectorFieldFlow, t_switch::Time)
-    return (x::State, v::Variable, t::Time) -> (t < t_switch ? F.rhs(x, v, t) : G.rhs(x, v, t))
+    return (x::State, v::Variable, t::Time) ->
+        (t < t_switch ? F.rhs(x, v, t) : G.rhs(x, v, t))
 end
 
 function __concat_rhs(F::ODEFlow, G::ODEFlow, t_switch::Time)
@@ -23,16 +26,18 @@ end
 
 function __concat_feedback_control(F::AbstractFlow, G::AbstractFlow, t_switch::Time)
     function _feedback_control(t, x, u, v)
-        t < t_switch ? F.feedback_control(t, x, u, v) : G.feedback_control(t, x, u, v)
+        return if t < t_switch
+            F.feedback_control(t, x, u, v)
+        else
+            G.feedback_control(t, x, u, v)
+        end
     end
     feedback_control = ControlLaw(_feedback_control, NonAutonomous, NonFixed)
     return feedback_control
 end
 
 function __concat_jumps(
-    F::AbstractFlow,
-    G::AbstractFlow,
-    jump::Union{Nothing, Tuple{Time, Any}} = nothing,
+    F::AbstractFlow, G::AbstractFlow, jump::Union{Nothing,Tuple{Time,Any}}=nothing
 )
     jumps = F.jumps
     append!(jumps, G.jumps)
@@ -42,7 +47,7 @@ end
 
 # --------------------------------------------------------------------------------------------
 # concatenate two flows with a prescribed switching time
-function concatenate(F::TF, g::Tuple{ctNumber, TF})::TF where {TF <: AbstractFlow}
+function concatenate(F::TF, g::Tuple{ctNumber,TF})::TF where {TF<:AbstractFlow}
     t_switch, G = g
     rhs! = __concat_rhs(F, G, t_switch)       # concatenation of the right and sides
     tstops = __concat_tstops(F, G, t_switch)    # times to break integration
@@ -50,7 +55,7 @@ function concatenate(F::TF, g::Tuple{ctNumber, TF})::TF where {TF <: AbstractFlo
     return TF(F.f, rhs!, tstops, jumps)  # we choose default values and options of F
 end
 
-function concatenate(F::TF, g::Tuple{ctNumber, Any, TF})::TF where {TF <: AbstractFlow}
+function concatenate(F::TF, g::Tuple{ctNumber,Any,TF})::TF where {TF<:AbstractFlow}
     t_switch, η_switch, G = g
     rhs! = __concat_rhs(F, G, t_switch)       # concatenation of the right and sides
     tstops = __concat_tstops(F, G, t_switch)    # times to break integration
@@ -60,7 +65,7 @@ end
 
 # --------------------------------------------------------------------------------------------
 # concatenate two flows with a prescribed switching time
-function concatenate(F::TF, g::Tuple{ctNumber, TF})::TF where {TF <: OptimalControlFlow}
+function concatenate(F::TF, g::Tuple{ctNumber,TF})::TF where {TF<:OptimalControlFlow}
     t_switch, G = g
     rhs! = __concat_rhs(F, G, t_switch)               # concatenation of the right and sides
     tstops = __concat_tstops(F, G, t_switch)            # times to break integration
@@ -69,7 +74,7 @@ function concatenate(F::TF, g::Tuple{ctNumber, TF})::TF where {TF <: OptimalCont
     return OptimalControlFlow(F.f, rhs!, feedback_u, F.ocp, F.kwargs_Flow, tstops, jumps) # we choose default values and options of F
 end
 
-function concatenate(F::TF, g::Tuple{ctNumber, Any, TF})::TF where {TF <: OptimalControlFlow}
+function concatenate(F::TF, g::Tuple{ctNumber,Any,TF})::TF where {TF<:OptimalControlFlow}
     t_switch, η_switch, G = g
     rhs! = __concat_rhs(F, G, t_switch)               # concatenation of the right and sides
     tstops = __concat_tstops(F, G, t_switch)            # times to break integration
@@ -79,5 +84,5 @@ function concatenate(F::TF, g::Tuple{ctNumber, Any, TF})::TF where {TF <: Optima
 end
 
 # --------------------------------------------------------------------------------------------
-*(F::TF, g::Tuple{ctNumber, TF}) where {TF <: AbstractFlow} = concatenate(F, g)
-*(F::TF, g::Tuple{ctNumber, Any, TF}) where {TF <: AbstractFlow} = concatenate(F, g)
+*(F::TF, g::Tuple{ctNumber,TF}) where {TF<:AbstractFlow} = concatenate(F, g)
+*(F::TF, g::Tuple{ctNumber,Any,TF}) where {TF<:AbstractFlow} = concatenate(F, g)
diff --git a/ext/function.jl b/ext/function.jl
index 7da5a26..48443ba 100644
--- a/ext/function.jl
+++ b/ext/function.jl
@@ -7,34 +7,36 @@ function ode_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
 
     # kwargs has priority wrt kwargs_flow
     function f(
-        tspan::Tuple{Time, Time},
+        tspan::Tuple{Time,Time},
         x0,
-        v = nothing;
+        v=nothing;
         jumps,
         _t_stops_interne,
         DiffEqRHS,
-        tstops = __tstops(),
-        callback = __callback(),
+        tstops=__tstops(),
+        callback=__callback(),
         kwargs...,
     )
 
         # ode
-        ode =
-            isnothing(v) ? OrdinaryDiffEq.ODEProblem(DiffEqRHS, x0, tspan) :
+        ode = if isnothing(v)
+            OrdinaryDiffEq.ODEProblem(DiffEqRHS, x0, tspan)
+        else
             OrdinaryDiffEq.ODEProblem(DiffEqRHS, x0, tspan, v)
+        end
 
         # jumps and callbacks
         cb, t_stops_all = __callbacks(callback, jumps, nothing, _t_stops_interne, tstops)
 
         # solve
         sol = OrdinaryDiffEq.solve(
-            ode,
-            alg = alg,
-            abstol = abstol,
-            reltol = reltol,
-            saveat = saveat,
-            tstops = t_stops_all,
-            callback = cb;
+            ode;
+            alg=alg,
+            abstol=abstol,
+            reltol=reltol,
+            saveat=saveat,
+            tstops=t_stops_all,
+            callback=cb,
             kwargs_Flow...,
             kwargs...,
         )
@@ -42,7 +44,7 @@ function ode_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
         return sol
     end
 
-    function f(t0::Time, x0, tf::Time, v = nothing; kwargs...)
+    function f(t0::Time, x0, tf::Time, v=nothing; kwargs...)
         sol = f((t0, tf), x0, v; kwargs...)
         return sol.u[end]
     end
@@ -53,12 +55,12 @@ end
 # --------------------------------------------------------------------------------------------
 function CTFlows.Flow(
     dyn::Function;
-    autonomous = true,
-    variable = false,
-    alg = __alg(),
-    abstol = __abstol(),
-    reltol = __reltol(),
-    saveat = __saveat(),
+    autonomous=true,
+    variable=false,
+    alg=__alg(),
+    abstol=__abstol(),
+    reltol=__reltol(),
+    saveat=__saveat(),
     kwargs_Flow...,
 )
     #
diff --git a/ext/hamiltonian.jl b/ext/hamiltonian.jl
index 8b35539..4e1837c 100644
--- a/ext/hamiltonian.jl
+++ b/ext/hamiltonian.jl
@@ -6,15 +6,15 @@ Returns a function that solves ODE problem associated to Hamiltonian vector fiel
 """
 function hamiltonian_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
     function f(
-        tspan::Tuple{Time, Time},
+        tspan::Tuple{Time,Time},
         x0::State,
         p0::Costate,
-        v::Variable = __variable(x0, p0);
+        v::Variable=__variable(x0, p0);
         jumps,
         _t_stops_interne,
         DiffEqRHS,
-        tstops = __tstops(),
-        callback = __callback(),
+        tstops=__tstops(),
+        callback=__callback(),
         kwargs...,
     )
 
@@ -23,17 +23,19 @@ function hamiltonian_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
 
         # jumps and callbacks
         n = size(x0, 1)
-        cb, t_stops_all = __callbacks(callback, jumps, rg(n + 1, 2n), _t_stops_interne, tstops)
+        cb, t_stops_all = __callbacks(
+            callback, jumps, rg(n + 1, 2n), _t_stops_interne, tstops
+        )
 
         # solve
         sol = OrdinaryDiffEq.solve(
-            ode,
-            alg = alg,
-            abstol = abstol,
-            reltol = reltol,
-            saveat = saveat,
-            tstops = t_stops_all,
-            callback = cb;
+            ode;
+            alg=alg,
+            abstol=abstol,
+            reltol=reltol,
+            saveat=saveat,
+            tstops=t_stops_all,
+            callback=cb,
             kwargs_Flow...,
             kwargs...,
         )
@@ -46,7 +48,7 @@ function hamiltonian_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
         x0::State,
         p0::Costate,
         tf::Time,
-        v::Variable = __variable(x0, p0);
+        v::Variable=__variable(x0, p0);
         kwargs...,
     )
         sol = f((t0, tf), x0, p0, v; kwargs...)
@@ -68,7 +70,7 @@ function rhs(h::AbstractHamiltonian)
         foo(z) = h(t, z[rg(1, n)], z[rg(n + 1, 2n)], v)
         dh = ctgradient(foo, z)
         dz[1:n] = dh[(n + 1):(2n)]
-        dz[(n + 1):(2n)] = -dh[1:n]
+        return dz[(n + 1):(2n)] = -dh[1:n]
     end
     return rhs!
 end
@@ -77,10 +79,10 @@ end
 # Flow from a Hamiltonian
 function CTFlows.Flow(
     h::AbstractHamiltonian;
-    alg = __alg(),
-    abstol = __abstol(),
-    reltol = __reltol(),
-    saveat = __saveat(),
+    alg=__alg(),
+    abstol=__abstol(),
+    reltol=__reltol(),
+    saveat=__saveat(),
     kwargs_Flow...,
 )
     #
@@ -93,17 +95,17 @@ end
 # Flow from a Hamiltonian Vector Field
 function CTFlows.Flow(
     hv::HamiltonianVectorField;
-    alg = __alg(),
-    abstol = __abstol(),
-    reltol = __reltol(),
-    saveat = __saveat(),
+    alg=__alg(),
+    abstol=__abstol(),
+    reltol=__reltol(),
+    saveat=__saveat(),
     kwargs_Flow...,
 )
     #
     f = hamiltonian_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
     function rhs!(dz::DCoTangent, z::CoTangent, v::Variable, t::Time)
         n = size(z, 1) ÷ 2
-        dz[rg(1, n)], dz[rg(n + 1, 2n)] = hv(t, z[rg(1, n)], z[rg(n + 1, 2n)], v)
+        return dz[rg(1, n)], dz[rg(n + 1, 2n)] = hv(t, z[rg(1, n)], z[rg(n + 1, 2n)], v)
     end
     return HamiltonianFlow(f, rhs!)
 end
diff --git a/ext/optimal_control_problem.jl b/ext/optimal_control_problem.jl
index ca2fad6..501e97e 100644
--- a/ext/optimal_control_problem.jl
+++ b/ext/optimal_control_problem.jl
@@ -15,14 +15,14 @@ julia> f = Flow(ocp, (x, p) -> p)
     The dimension of the output of the control function must be consistent with the dimension usage of the control of the optimal control problem.
 """
 function CTFlows.Flow(
-    ocp::OptimalControlModel{T, V},
-    u_::Union{Function, ControlLaw{T, V}};
-    alg = __alg(),
-    abstol = __abstol(),
-    reltol = __reltol(),
-    saveat = __saveat(),
+    ocp::OptimalControlModel{T,V},
+    u_::Union{Function,ControlLaw{T,V}};
+    alg=__alg(),
+    abstol=__abstol(),
+    reltol=__reltol(),
+    saveat=__saveat(),
     kwargs_Flow...,
-) where {T, V}
+) where {T,V}
     h, u = __create_hamiltonian(ocp, u_) # construction of the Hamiltonian
     return __ocp_Flow(ocp, h, u, alg, abstol, reltol, saveat; kwargs_Flow...)
 end
@@ -46,23 +46,23 @@ julia> f = Flow(ocp, (t, x, p) -> p[1], (t, x, u) -> x[1] - 1, (t, x, p) -> x[1]
     The dimension of the output of the control function must be consistent with the dimension usage of the control of the optimal control problem.
 """
 function CTFlows.Flow(
-    ocp::OptimalControlModel{T, V},
-    u_::Union{Function, ControlLaw{T, V}, FeedbackControl{T, V}},
-    g_::Union{Function, MixedConstraint{T, V}, StateConstraint{T, V}},
-    μ_::Union{Function, Multiplier{T, V}};
-    alg = __alg(),
-    abstol = __abstol(),
-    reltol = __reltol(),
-    saveat = __saveat(),
+    ocp::OptimalControlModel{T,V},
+    u_::Union{Function,ControlLaw{T,V},FeedbackControl{T,V}},
+    g_::Union{Function,MixedConstraint{T,V},StateConstraint{T,V}},
+    μ_::Union{Function,Multiplier{T,V}};
+    alg=__alg(),
+    abstol=__abstol(),
+    reltol=__reltol(),
+    saveat=__saveat(),
     kwargs_Flow...,
-) where {T, V}
+) where {T,V}
     h, u = __create_hamiltonian(ocp, u_, g_, μ_) # construction of the Hamiltonian
     return __ocp_Flow(ocp, h, u, alg, abstol, reltol, saveat; kwargs_Flow...)
 end
 
 # ---------------------------------------------------------------------------------------------------
 function __ocp_Flow(
-    ocp::OptimalControlModel{T, V},
+    ocp::OptimalControlModel{T,V},
     h::Hamiltonian,
     u::ControlLaw,
     alg,
@@ -70,9 +70,9 @@ function __ocp_Flow(
     reltol,
     saveat;
     kwargs_Flow...,
-) where {T, V}
+) where {T,V}
     rhs! = rhs(h) # right and side: same as for a flow from a Hamiltonian
     f = hamiltonian_usage(alg, abstol, reltol, saveat; kwargs_Flow...) # flow function
-    kwargs_Flow = (kwargs_Flow..., alg = alg, abstol = abstol, reltol = reltol, saveat = saveat)
+    kwargs_Flow = (kwargs_Flow..., alg=alg, abstol=abstol, reltol=reltol, saveat=saveat)
     return OptimalControlFlow(f, rhs!, u, ocp, kwargs_Flow)
 end
diff --git a/ext/types.jl b/ext/types.jl
index 11d92e4..8dee5e0 100644
--- a/ext/types.jl
+++ b/ext/types.jl
@@ -1,65 +1,65 @@
 # ---------------------------------------------------------------------------------------------------
 # This is the flow returned by the function Flow
 # The call to the flow is given after.
-struct HamiltonianFlow <: AbstractFlow{DCoTangent, CoTangent}
+struct HamiltonianFlow <: AbstractFlow{DCoTangent,CoTangent}
     f::Function      # f(args..., rhs): compute the flow
     rhs!::Function   # DifferentialEquations rhs
     tstops::Times    # stopping times
-    jumps::Vector{Tuple{Time, Costate}} # specific jumps the integrator must perform
+    jumps::Vector{Tuple{Time,Costate}} # specific jumps the integrator must perform
     function HamiltonianFlow(
         f,
         rhs!,
-        tstops::Times = Vector{Time}(),
-        jumps::Vector{Tuple{Time, Costate}} = Vector{Tuple{Time, Costate}}(),
+        tstops::Times=Vector{Time}(),
+        jumps::Vector{Tuple{Time,Costate}}=Vector{Tuple{Time,Costate}}(),
     )
         return new(f, rhs!, tstops, jumps)
     end
 end
 
 # call F.f
-(F::HamiltonianFlow)(args...; kwargs...) = begin
-    F.f(args...; jumps = F.jumps, _t_stops_interne = F.tstops, DiffEqRHS = F.rhs!, kwargs...)
+function (F::HamiltonianFlow)(args...; kwargs...)
+    F.f(args...; jumps=F.jumps, _t_stops_interne=F.tstops, DiffEqRHS=F.rhs!, kwargs...)
 end
 
 # ---------------------------------------------------------------------------------------------------
-struct VectorFieldFlow <: AbstractFlow{DState, State}
+struct VectorFieldFlow <: AbstractFlow{DState,State}
     f::Function     # f(args..., rhs): compute the flow
     rhs::Function  # DifferentialEquations rhs
     tstops::Times   # stopping times
-    jumps::Vector{Tuple{Time, State}} # specific jumps the integrator must perform
+    jumps::Vector{Tuple{Time,State}} # specific jumps the integrator must perform
     function VectorFieldFlow(
         f,
         rhs,
-        tstops::Times = Vector{Time}(),
-        jumps::Vector{Tuple{Time, State}} = Vector{Tuple{Time, State}}(),
+        tstops::Times=Vector{Time}(),
+        jumps::Vector{Tuple{Time,State}}=Vector{Tuple{Time,State}}(),
     )
         return new(f, rhs, tstops, jumps)
     end
 end
 
 # call F.f
-(F::VectorFieldFlow)(args...; kwargs...) = begin
-    F.f(args...; jumps = F.jumps, _t_stops_interne = F.tstops, DiffEqRHS = F.rhs, kwargs...)
+function (F::VectorFieldFlow)(args...; kwargs...)
+    F.f(args...; jumps=F.jumps, _t_stops_interne=F.tstops, DiffEqRHS=F.rhs, kwargs...)
 end
 
 # ---------------------------------------------------------------------------------------------------
-struct ODEFlow <: AbstractFlow{Any, Any}
+struct ODEFlow <: AbstractFlow{Any,Any}
     f::Function     # f(args..., rhs): compute the flow
     rhs::Function   # DifferentialEquations rhs
     tstops::Times   # stopping times
-    jumps::Vector{Tuple{Time, Any}} # specific jumps the integrator must perform
+    jumps::Vector{Tuple{Time,Any}} # specific jumps the integrator must perform
     function ODEFlow(
         f,
         rhs!,
-        tstops::Times = Vector{Time}(),
-        jumps::Vector{Tuple{Time, Any}} = Vector{Tuple{Time, Any}}(),
+        tstops::Times=Vector{Time}(),
+        jumps::Vector{Tuple{Time,Any}}=Vector{Tuple{Time,Any}}(),
     )
         return new(f, rhs!, tstops, jumps)
     end
 end
 
-(F::ODEFlow)(args...; kwargs...) = begin
-    F.f(args...; jumps = F.jumps, _t_stops_interne = F.tstops, DiffEqRHS = F.rhs, kwargs...)
+function (F::ODEFlow)(args...; kwargs...)
+    F.f(args...; jumps=F.jumps, _t_stops_interne=F.tstops, DiffEqRHS=F.rhs, kwargs...)
 end
 
 # ---------------------------------------------------------------------------------------------------
@@ -117,14 +117,14 @@ function CTBase.OptimalControlSolution(ocfs::OptimalControlFlowSolution; kwargs.
     end
 
     # we provide the variable only if the problem is NonFixed
-    kwargs_OCS = CTBase.is_fixed(ocp) ? () : (variable = v,)
+    kwargs_OCS = CTBase.is_fixed(ocp) ? () : (variable=v,)
     kwargs_OCS = (
         kwargs_OCS...,
-        time_grid = T,
-        state = t -> x(t),
-        costate = t -> p(t),
-        control = t -> u(t),
-        objective = obj,
+        time_grid=T,
+        state=t -> x(t),
+        costate=t -> p(t),
+        control=t -> u(t),
+        objective=obj,
     )
     sol = CTBase.OptimalControlSolution(ocp; kwargs_OCS...)
 
@@ -132,16 +132,16 @@ function CTBase.OptimalControlSolution(ocfs::OptimalControlFlowSolution; kwargs.
 end
 
 # ---------------------------------------------------------------------------------------------------
-struct OptimalControlFlow{VD} <: AbstractFlow{DCoTangent, CoTangent}
+struct OptimalControlFlow{VD} <: AbstractFlow{DCoTangent,CoTangent}
     # 
     f::Function      # the mere function which depends on the kind of flow (Hamiltonian or classical) 
     # this function takes a right and side as input
     rhs!::Function   # the right and side of the form: rhs!(du::D, u::U, p::V, t::T)
     tstops::Times    # specific times  the integrator must stop
     # useful when the rhs is not smooth at such times
-    jumps::Vector{Tuple{Time, Costate}} # specific jumps the integrator must perform
+    jumps::Vector{Tuple{Time,Costate}} # specific jumps the integrator must perform
     feedback_control::ControlLaw # the control law in feedback form, that is u(t, x, p, v)
-    ocp::OptimalControlModel{<:TimeDependence, VD} # the optimal control problem
+    ocp::OptimalControlModel{<:TimeDependence,VD} # the optimal control problem
     kwargs_Flow::Any     #
 
     # constructor
@@ -149,18 +149,29 @@ struct OptimalControlFlow{VD} <: AbstractFlow{DCoTangent, CoTangent}
         f::Function,
         rhs!::Function,
         u::ControlLaw,
-        ocp::OptimalControlModel{<:TimeDependence, VD},
+        ocp::OptimalControlModel{<:TimeDependence,VD},
         kwargs_Flow,
-        tstops::Times = Vector{Time}(),
-        jumps::Vector{Tuple{Time, Costate}} = Vector{Tuple{Time, Costate}}(),
-    ) where {VD <: VariableDependence}
+        tstops::Times=Vector{Time}(),
+        jumps::Vector{Tuple{Time,Costate}}=Vector{Tuple{Time,Costate}}(),
+    ) where {VD<:VariableDependence}
         return new{VD}(f, rhs!, tstops, jumps, u, ocp, kwargs_Flow)
     end
 end
 
 # call F.f
-function (F::OptimalControlFlow{Fixed})(t0::Time, x0::State, p0::Costate, tf::Time; kwargs...)
-    F.f(t0, x0, p0, tf; jumps = F.jumps, _t_stops_interne = F.tstops, DiffEqRHS = F.rhs!, kwargs...)
+function (F::OptimalControlFlow{Fixed})(
+    t0::Time, x0::State, p0::Costate, tf::Time; kwargs...
+)
+    return F.f(
+        t0,
+        x0,
+        p0,
+        tf;
+        jumps=F.jumps,
+        _t_stops_interne=F.tstops,
+        DiffEqRHS=F.rhs!,
+        kwargs...,
+    )
 end
 
 function (F::OptimalControlFlow{NonFixed})(
@@ -168,42 +179,40 @@ function (F::OptimalControlFlow{NonFixed})(
     x0::State,
     p0::Costate,
     tf::Time,
-    v::Variable = __variable(t0, x0, p0, tf, F.ocp);
+    v::Variable=__variable(t0, x0, p0, tf, F.ocp);
     kwargs...,
 )
-    F.f(
+    return F.f(
         t0,
         x0,
         p0,
         tf,
         v;
-        jumps = F.jumps,
-        _t_stops_interne = F.tstops,
-        DiffEqRHS = F.rhs!,
+        jumps=F.jumps,
+        _t_stops_interne=F.tstops,
+        DiffEqRHS=F.rhs!,
         kwargs...,
     )
 end
 
 # call F.f and then, construct an optimal control solution
-function (F::OptimalControlFlow{Fixed})(tspan::Tuple{Time, Time}, x0::State, p0::Costate; kwargs...)
+function (F::OptimalControlFlow{Fixed})(
+    tspan::Tuple{Time,Time}, x0::State, p0::Costate; kwargs...
+)
     ode_sol = F.f(
-        tspan,
-        x0,
-        p0;
-        jumps = F.jumps,
-        _t_stops_interne = F.tstops,
-        DiffEqRHS = F.rhs!,
-        kwargs...,
+        tspan, x0, p0; jumps=F.jumps, _t_stops_interne=F.tstops, DiffEqRHS=F.rhs!, kwargs...
+    )
+    flow_sol = OptimalControlFlowSolution(
+        ode_sol, F.feedback_control, F.ocp, __variable(x0, p0)
     )
-    flow_sol = OptimalControlFlowSolution(ode_sol, F.feedback_control, F.ocp, __variable(x0, p0))
     return CTBase.OptimalControlSolution(flow_sol; F.kwargs_Flow..., kwargs...)
 end
 
 function (F::OptimalControlFlow{NonFixed})(
-    tspan::Tuple{Time, Time},
+    tspan::Tuple{Time,Time},
     x0::State,
     p0::Costate,
-    v::Variable = __variable(tspan[1], x0, p0, tspan[2], F.ocp);
+    v::Variable=__variable(tspan[1], x0, p0, tspan[2], F.ocp);
     kwargs...,
 )
     ode_sol = F.f(
@@ -211,9 +220,9 @@ function (F::OptimalControlFlow{NonFixed})(
         x0,
         p0,
         v;
-        jumps = F.jumps,
-        _t_stops_interne = F.tstops,
-        DiffEqRHS = F.rhs!,
+        jumps=F.jumps,
+        _t_stops_interne=F.tstops,
+        DiffEqRHS=F.rhs!,
         kwargs...,
     )
     flow_sol = OptimalControlFlowSolution(ode_sol, F.feedback_control, F.ocp, v)
diff --git a/ext/utils.jl b/ext/utils.jl
index b8d7c4c..3953086 100644
--- a/ext/utils.jl
+++ b/ext/utils.jl
@@ -9,12 +9,15 @@ function __callbacks(callback, jumps, _rg, _t_stops_interne, tstops)
 
         # add the jump η from η_jumps on p with z = (x, p) at time t from t_jumps
         function condition(out, u, t, integrator)
-            out[:] = t_jumps .- t
+            return out[:] = t_jumps .- t
         end
 
         function affect!(integrator, event_index)
-            isnothing(_rg) ? integrator.u += η_jumps[event_index] :
-            integrator.u[_rg] += η_jumps[event_index]
+            return if isnothing(_rg)
+                integrator.u += η_jumps[event_index]
+            else
+                integrator.u[_rg] += η_jumps[event_index]
+            end
         end
 
         cbjumps = VectorContinuousCallback(condition, affect!, size(jumps, 1))
diff --git a/ext/vector_field.jl b/ext/vector_field.jl
index d4b03d3..1a024b3 100644
--- a/ext/vector_field.jl
+++ b/ext/vector_field.jl
@@ -7,14 +7,14 @@ function vector_field_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
 
     # kwargs has priority wrt kwargs_flow
     function f(
-        tspan::Tuple{Time, Time},
+        tspan::Tuple{Time,Time},
         x0::State,
-        v::Variable = __variable(x0);
+        v::Variable=__variable(x0);
         jumps,
         _t_stops_interne,
         DiffEqRHS,
-        tstops = __tstops(),
-        callback = __callback(),
+        tstops=__tstops(),
+        callback=__callback(),
         kwargs...,
     )
 
@@ -26,13 +26,13 @@ function vector_field_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
 
         # solve
         sol = OrdinaryDiffEq.solve(
-            ode,
-            alg = alg,
-            abstol = abstol,
-            reltol = reltol,
-            saveat = saveat,
-            tstops = t_stops_all,
-            callback = cb;
+            ode;
+            alg=alg,
+            abstol=abstol,
+            reltol=reltol,
+            saveat=saveat,
+            tstops=t_stops_all,
+            callback=cb,
             kwargs_Flow...,
             kwargs...,
         )
@@ -40,7 +40,7 @@ function vector_field_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
         return sol
     end
 
-    function f(t0::Time, x0::State, t::Time, v::Variable = __variable(x0); kwargs...)
+    function f(t0::Time, x0::State, t::Time, v::Variable=__variable(x0); kwargs...)
         sol = f((t0, t), x0, v; kwargs...)
         return sol.u[end]
     end
@@ -52,10 +52,10 @@ end
 # Flow of a vector field
 function CTFlows.Flow(
     vf::VectorField;
-    alg = __alg(),
-    abstol = __abstol(),
-    reltol = __reltol(),
-    saveat = __saveat(),
+    alg=__alg(),
+    abstol=__abstol(),
+    reltol=__reltol(),
+    saveat=__saveat(),
     kwargs_Flow...,
 )
     f = vector_field_usage(alg, abstol, reltol, saveat; kwargs_Flow...)
diff --git a/src/optimal_control_problem_utils.jl b/src/optimal_control_problem_utils.jl
index 0e9a575..17e71ed 100644
--- a/src/optimal_control_problem_utils.jl
+++ b/src/optimal_control_problem_utils.jl
@@ -1,22 +1,24 @@
 # ---------------------------------------------------------------------------------------------------
-function __create_hamiltonian(ocp::OptimalControlModel{T, V}, u::ControlLaw{T, V}) where {T, V}
+function __create_hamiltonian(ocp::OptimalControlModel{T,V}, u::ControlLaw{T,V}) where {T,V}
     f, f⁰, p⁰, s = __get_data_for_ocp_flow(ocp) # data
     @assert f ≠ nothing "no dynamics in ocp"
-    h = Hamiltonian(f⁰ ≠ nothing ? makeH(f, u, f⁰, p⁰, s) : makeH(f, u), NonAutonomous, NonFixed)
+    h = Hamiltonian(
+        f⁰ ≠ nothing ? makeH(f, u, f⁰, p⁰, s) : makeH(f, u), NonAutonomous, NonFixed
+    )
     return h, u
 end
 
-function __create_hamiltonian(ocp::OptimalControlModel{T, V}, u::Function) where {T, V}
+function __create_hamiltonian(ocp::OptimalControlModel{T,V}, u::Function) where {T,V}
     return __create_hamiltonian(ocp, ControlLaw(u, T, V))
 end
 
 # ---------------------------------------------------------------------------------------------------
 function __create_hamiltonian(
-    ocp::OptimalControlModel{T, V},
-    u::ControlLaw{T, V},
-    g::MixedConstraint{T, V},
-    μ::Multiplier{T, V},
-) where {T, V}
+    ocp::OptimalControlModel{T,V},
+    u::ControlLaw{T,V},
+    g::MixedConstraint{T,V},
+    μ::Multiplier{T,V},
+) where {T,V}
     f, f⁰, p⁰, s = __get_data_for_ocp_flow(ocp) # data
     @assert f ≠ nothing "no dynamics in ocp"
     h = Hamiltonian(
@@ -27,16 +29,13 @@ function __create_hamiltonian(
     return h, u
 end
 
-function __create_hamiltonian(ocp::OptimalControlModel{T, V}, u::Function, g, μ) where {T, V}
+function __create_hamiltonian(ocp::OptimalControlModel{T,V}, u::Function, g, μ) where {T,V}
     return __create_hamiltonian(ocp, ControlLaw(u, T, V), g, μ)
 end
 
 function __create_hamiltonian(
-    ocp::OptimalControlModel{T, V},
-    u_::FeedbackControl{T, V},
-    g,
-    μ,
-) where {T, V}
+    ocp::OptimalControlModel{T,V}, u_::FeedbackControl{T,V}, g, μ
+) where {T,V}
     u = @match (T, V) begin
         (Autonomous, Fixed) => ControlLaw((x, p) -> u_(x), T, V)
         (Autonomous, NonFixed) => ControlLaw((x, p, v) -> u_(x, v), T, V)
@@ -47,20 +46,14 @@ function __create_hamiltonian(
 end
 
 function __create_hamiltonian(
-    ocp::OptimalControlModel{T, V},
-    u::ControlLaw{T, V},
-    g::Function,
-    μ,
-) where {T, V}
+    ocp::OptimalControlModel{T,V}, u::ControlLaw{T,V}, g::Function, μ
+) where {T,V}
     return __create_hamiltonian(ocp, u, MixedConstraint(g, T, V), μ)
 end
 
 function __create_hamiltonian(
-    ocp::OptimalControlModel{T, V},
-    u::ControlLaw{T, V},
-    g_::StateConstraint{T, V},
-    μ,
-) where {T, V}
+    ocp::OptimalControlModel{T,V}, u::ControlLaw{T,V}, g_::StateConstraint{T,V}, μ
+) where {T,V}
     g = @match (T, V) begin
         (Autonomous, Fixed) => MixedConstraint((x, u) -> g_(x), T, V)
         (Autonomous, NonFixed) => MixedConstraint((x, u, v) -> g_(x, v), T, V)
@@ -71,11 +64,8 @@ function __create_hamiltonian(
 end
 
 function __create_hamiltonian(
-    ocp::OptimalControlModel{T, V},
-    u::ControlLaw{T, V},
-    g::MixedConstraint{T, V},
-    μ::Function,
-) where {T, V}
+    ocp::OptimalControlModel{T,V}, u::ControlLaw{T,V}, g::MixedConstraint{T,V}, μ::Function
+) where {T,V}
     return __create_hamiltonian(ocp, u, g, Multiplier(μ, T, V))
 end
 
@@ -148,7 +138,9 @@ function makeH(
 )
     function H(t, x, p, v)
         u_ = u(t, x, p, v)
-        return p' * f(t, x, u_, v) + s * p⁰ * f⁰(t, x, u_, v) + μ(t, x, p, v)' * g(t, x, u_, v)
+        return p' * f(t, x, u_, v) +
+               s * p⁰ * f⁰(t, x, u_, v) +
+               μ(t, x, p, v)' * g(t, x, u_, v)
     end
     return H
 end
diff --git a/test/manual_plot_goddard.jl b/test/manual_plot_goddard.jl
index 2f7ae21..619bde4 100644
--- a/test/manual_plot_goddard.jl
+++ b/test/manual_plot_goddard.jl
@@ -17,16 +17,16 @@ mf = 0.6
 x0 = [r0, v0, m0]
 
 # OCP model
-ocp = Model(variable = true)
+ocp = Model(; variable=true)
 variable!(ocp, 1)
-time!(ocp, t0 = t0, indf = 1) # if not provided, final time is free
+time!(ocp; t0=t0, indf=1) # if not provided, final time is free
 state!(ocp, 3, "x", ["r", "v", "m"]) # state dim
 control!(ocp, 1) # control dim
-constraint!(ocp, :initial, lb = x0, ub = x0)
-constraint!(ocp, :control, f = (u, v) -> u, lb = 0, ub = 1)
-constraint!(ocp, :state, f = (x, v) -> x[1], lb = r0, ub = Inf, label = :state_con1)
-constraint!(ocp, :state, f = (x, v) -> x[2], lb = 0, ub = vmax, label = :state_con2)
-constraint!(ocp, :state, f = (x, v) -> x[3], lb = m0, ub = mf, label = :state_con3)
+constraint!(ocp, :initial; lb=x0, ub=x0)
+constraint!(ocp, :control; f=(u, v) -> u, lb=0, ub=1)
+constraint!(ocp, :state; f=(x, v) -> x[1], lb=r0, ub=Inf, label=:state_con1)
+constraint!(ocp, :state; f=(x, v) -> x[2], lb=0, ub=vmax, label=:state_con2)
+constraint!(ocp, :state; f=(x, v) -> x[3], lb=m0, ub=mf, label=:state_con3)
 objective!(ocp, :mayer, (x0, xf, v) -> xf[1], :max)
 function F0(x)
     r, v, m = x
@@ -75,7 +75,7 @@ tf = 0.20204744057041196
 f1sb0 = f1 * (t1, fs) * (t2, fb) * (t3, f0) # concatenation of the Hamiltonian flows
 flow_sol = f1sb0((t0, tf), x0, p0)
 
-pp = plot(flow_sol, size = (900, 600))
+pp = plot(flow_sol; size=(900, 600))
 
 # Abstract model
 f(x, u, v) = F0(x) + u * F1(x)
diff --git a/test/manual_plot_integrator.jl b/test/manual_plot_integrator.jl
index 5852a82..eaacbd0 100644
--- a/test/manual_plot_integrator.jl
+++ b/test/manual_plot_integrator.jl
@@ -23,8 +23,8 @@ f = Flow(ocp, u)
 
 xf, pf = f(t0, x0, p0, tf)
 
-sol = f((t0, tf), x0, p0; saveat = range(t0, tf, 101))
+sol = f((t0, tf), x0, p0; saveat=range(t0, tf, 101))
 
-plot(sol, size = (900, 600))
+plot(sol; size=(900, 600))
 
 #nothing
diff --git a/test/runtests.jl b/test/runtests.jl
index eb2939b..2418f38 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -10,7 +10,7 @@ using LinearAlgebra
 const CTFlowsODE = Base.get_extension(CTFlows, :CTFlowsODE) # to test functions from CTFlowsODE not in CTFlows
 
 @testset verbose = true showtiming = true "CTFlows" begin
-    for name ∈ (
+    for name in (
         :aqua,
         :concatenation,
         :default,
diff --git a/test/test_aqua.jl b/test/test_aqua.jl
index 2afee05..6cb84aa 100644
--- a/test/test_aqua.jl
+++ b/test/test_aqua.jl
@@ -2,10 +2,10 @@ function test_aqua()
     @testset "Aqua.jl" begin
         Aqua.test_all(
             CTFlows;
-            ambiguities = false,
+            ambiguities=false,
             #stale_deps=(ignore=[:SomePackage],),
-            deps_compat = (ignore = [:LinearAlgebra, :Unicode],),
-            piracies = true,
+            deps_compat=(ignore=[:LinearAlgebra, :Unicode],),
+            piracies=true,
         )
         # do not warn about ambiguities in dependencies
         Aqua.test_ambiguities(CTFlows)
diff --git a/test/test_concatenation.jl b/test/test_concatenation.jl
index e9aa3a9..f5eb087 100644
--- a/test/test_concatenation.jl
+++ b/test/test_concatenation.jl
@@ -39,7 +39,7 @@ function test_concatenation()
         #
         f1 = Flow(Hamiltonian(H1))
         f2 = Flow(Hamiltonian(H2))
-        f3 = Flow(Hamiltonian(H3, autonomous = false))
+        f3 = Flow(Hamiltonian(H3; autonomous=false))
 
         # one flow is used because t1 > tf
         f = f1 * (2tf, f2)
@@ -69,7 +69,7 @@ function test_concatenation()
         f = f1 * ((t0 + tf) / 4, f1) * ((t0 + tf) / 2, f1)
         N = 100
         saveat = range(t0, tf, N)
-        sol = f((t0, tf), x0, p0, saveat = saveat)
+        sol = f((t0, tf), x0, p0; saveat=saveat)
         xf = sol.u[end][1:n]
         pf = sol.u[end][(n + 1):(2n)]
         Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@@ -84,7 +84,7 @@ function test_concatenation()
         #
         f1 = Flow(HamiltonianVectorField(Hv1))
         f2 = Flow(HamiltonianVectorField(Hv2))
-        f3 = Flow(HamiltonianVectorField(Hv3, autonomous = false))
+        f3 = Flow(HamiltonianVectorField(Hv3; autonomous=false))
 
         # one flow is used because t1 > tf
         f = f1 * (2tf, f2)
@@ -114,7 +114,7 @@ function test_concatenation()
         f = f1 * ((t0 + tf) / 4, f1) * ((t0 + tf) / 2, f1)
         N = 100
         saveat = range(t0, tf, N)
-        sol = f((t0, tf), x0, p0, saveat = saveat)
+        sol = f((t0, tf), x0, p0; saveat=saveat)
         xf = sol.u[end][1:n]
         pf = sol.u[end][(n + 1):(2n)]
         Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@@ -129,7 +129,7 @@ function test_concatenation()
         #
         f1 = Flow(VectorField(V1))
         f2 = Flow(VectorField(V2))
-        f3 = Flow(VectorField(V3, autonomous = false))
+        f3 = Flow(VectorField(V3; autonomous=false))
 
         # one flow is used because t1 > tf
         f = f1 * (2tf, f2)
@@ -155,7 +155,7 @@ function test_concatenation()
         f = f1 * ((t0 + tf) / 4, f1) * ((t0 + tf) / 2, f1)
         N = 100
         saveat = range(t0, tf, N)
-        sol = f((t0, tf), [x0; p0], saveat = saveat)
+        sol = f((t0, tf), [x0; p0]; saveat=saveat)
         xf = sol.u[end][1:n]
         pf = sol.u[end][(n + 1):(2n)]
         Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@@ -170,7 +170,7 @@ function test_concatenation()
         #
         f1 = Flow(V1)
         f2 = Flow(V2)
-        f3 = Flow(V3, autonomous = false)
+        f3 = Flow(V3; autonomous=false)
 
         # one flow is used because t1 > tf
         f = f1 * (2tf, f2)
@@ -196,7 +196,7 @@ function test_concatenation()
         f = f1 * ((t0 + tf) / 4, f1) * ((t0 + tf) / 2, f1)
         N = 100
         saveat = range(t0, tf, N)
-        sol = f((t0, tf), [x0; p0], saveat = saveat)
+        sol = f((t0, tf), [x0; p0]; saveat=saveat)
         xf = sol.u[end][1:n]
         pf = sol.u[end][(n + 1):(2n)]
         Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@@ -211,7 +211,7 @@ function test_concatenation()
         # Hamiltonien
         f1 = Flow(Hamiltonian(H1))
         f2 = Flow(Hamiltonian(H2))
-        f3 = Flow(Hamiltonian(H3, autonomous = false))
+        f3 = Flow(Hamiltonian(H3; autonomous=false))
         f = f1 * ((t0 + tf) / 4, [0, 0], f2) * ((t0 + tf) / 2, f3)
         xf, pf = f(t0, x0, p0, tf + (t0 + tf) / 2)
         Test.@test xf ≈ [x1_sol(tf), x2_sol(tf)] atol = 1e-5
@@ -220,7 +220,7 @@ function test_concatenation()
         # vector field
         f1 = Flow(VectorField(V1))
         f2 = Flow(VectorField(V2))
-        f3 = Flow(VectorField(V3, autonomous = false))
+        f3 = Flow(VectorField(V3; autonomous=false))
         f = f1 * ((t0 + tf) / 4, [0, 0, 0, 0], f2) * ((t0 + tf) / 2, f3)
         zf = f(t0, [x0; p0], tf + (t0 + tf) / 2)
         Test.@test zf ≈ [x1_sol(tf), x2_sol(tf), p1_sol(tf), p2_sol(tf)] atol = 1e-5
@@ -230,7 +230,7 @@ function test_concatenation()
 
         # example from https://docs.sciml.ai/DiffEqDocs/stable/features/callback_functions/#DiscreteCallback-Examples
         function dyn(du, u, p, t)
-            du[1] = -u[1]
+            return du[1] = -u[1]
         end
         u0 = [10.0]
         prob = ODEProblem(dyn, u0, (0.0, 10.0))
@@ -238,7 +238,7 @@ function test_concatenation()
         condition(u, t, integrator) = t ∈ dosetimes
         affect!(integrator) = integrator.u[1] += 10
         cb = DiscreteCallback(condition, affect!)
-        sol = solve(prob, Tsit5(), callback = cb, tstops = dosetimes)
+        sol = solve(prob, Tsit5(); callback=cb, tstops=dosetimes)
 
         #
         x0 = 10
@@ -257,12 +257,20 @@ function test_concatenation()
         # test
         N = 100
         tspan = range(0, 10, N)
-        Test.@test norm([sol(t)[1] - sol2(t) for t ∈ tspan]) / N ≈ 0 atol = 1e-3
-        Test.@test norm([sol(t)[1] - sol3(t) for t ∈ tspan]) / N ≈ 0 atol = 1e-3
+        Test.@test norm([sol(t)[1] - sol2(t) for t in tspan]) / N ≈ 0 atol = 1e-3
+        Test.@test norm([sol(t)[1] - sol3(t) for t in tspan]) / N ≈ 0 atol = 1e-3
 
         # -------
         f = Flow(Hamiltonian((x, p) -> 0.5p^2))
-        fc = f * (1, 1, f) * (1.5, f) * (2, 1, f) * (2.5, f) * (3, 1, f) * (3.5, f) * (4, 1, f)
+        fc =
+            f *
+            (1, 1, f) *
+            (1.5, f) *
+            (2, 1, f) *
+            (2.5, f) *
+            (3, 1, f) *
+            (3.5, f) *
+            (4, 1, f)
         xf, pf = fc(0, 0, 0, 5)
         Test.@test xf ≈ 10 atol = 1e-6
         Test.@test pf ≈ 4 atol = 1e-6
@@ -287,8 +295,8 @@ function test_concatenation()
         ocp = Model()
         state!(ocp, 2)
         control!(ocp, 2)
-        time!(ocp, t0 = 0, tf = 5)
-        constraint!(ocp, :initial, lb = [0, 0], ub = [0, 0])
+        time!(ocp; t0=0, tf=5)
+        constraint!(ocp, :initial; lb=[0, 0], ub=[0, 0])
         dynamics!(ocp, (x, u) -> u)
         objective!(ocp, :mayer, (x0, xf) -> xf)
         f = Flow(ocp, (x, p) -> [p[1] / 2, 0])
@@ -310,10 +318,10 @@ function test_concatenation()
         ocp = Model()
         state!(ocp, 1)
         control!(ocp, 1)
-        time!(ocp, t0 = 0, tf = 1)
-        constraint!(ocp, :initial, lb = -1, ub = -1, label = :initial_constraint)
-        constraint!(ocp, :final, lb = 0, ub = 0, label = :final_constraint)
-        constraint!(ocp, :control, lb = -1, ub = 1, label = :control_constraint)
+        time!(ocp; t0=0, tf=1)
+        constraint!(ocp, :initial; lb=-1, ub=-1, label=:initial_constraint)
+        constraint!(ocp, :final; lb=0, ub=0, label=:final_constraint)
+        constraint!(ocp, :control; lb=-1, ub=1, label=:control_constraint)
         dynamics!(ocp, (x, u) -> -x + u)
         objective!(ocp, :lagrange, (x, u) -> abs(u))
         f0 = Flow(ocp, ControlLaw((x, p) -> 0))
diff --git a/test/test_flow_function.jl b/test/test_flow_function.jl
index 5ed0e12..acd4cc3 100644
--- a/test/test_flow_function.jl
+++ b/test/test_flow_function.jl
@@ -12,7 +12,7 @@ function test_flow_function()
 
     @testset "4D non autonomous, variable" begin
         V(t, z, l) = [z[2], (2 + l) * z[2 + 2], 0.0, -z[2 + 1]]
-        z = Flow(V, autonomous = false, variable = true)
+        z = Flow(V; autonomous=false, variable=true)
         t0 = 0.0
         tf = 1.0
         x0 = [-1.0, 0.0]
@@ -45,12 +45,12 @@ function test_flow_function()
         #
         S = solve(
             ODEProblem((S, _, _) -> ricatti(S), zeros(size(A)), (tf, t0)),
-            Tsit5(),
-            reltol = 1e-12,
-            abstol = 1e-12,
+            Tsit5();
+            reltol=1e-12,
+            abstol=1e-12,
         )
         #
-        f = Flow(ricatti, autonomous = true, variable = false)
+        f = Flow(ricatti; autonomous=true, variable=false)
         SS = f((tf, t0), zeros(size(A)))
         #
         Test.@test S.u[end] ≈ SS.u[end] atol = 1e-5
@@ -60,12 +60,12 @@ function test_flow_function()
         #
         x = solve(
             ODEProblem((x, _, t) -> dyn(t, x), x0, (t0, tf)),
-            Tsit5(),
-            reltol = 1e-8,
-            abstol = 1e-8,
+            Tsit5();
+            reltol=1e-8,
+            abstol=1e-8,
         )
         #
-        f = Flow(dyn, autonomous = false, variable = false)
+        f = Flow(dyn; autonomous=false, variable=false)
         xx = f((t0, tf), x0)
         #
         Test.@test x.u[end] ≈ xx.u[end] atol = 1e-5
@@ -78,12 +78,12 @@ function test_flow_function()
         #
         p = solve(
             ODEProblem((p, _, t) -> ϕ(t, p), zeros(2), (tf, t0)),
-            Tsit5(),
-            reltol = 1e-8,
-            abstol = 1e-8,
+            Tsit5();
+            reltol=1e-8,
+            abstol=1e-8,
         )
         #
-        f = Flow(ϕ, autonomous = false, variable = false)
+        f = Flow(ϕ; autonomous=false, variable=false)
         pp = f((tf, t0), zeros(2))
         #
         Test.@test p.u[end] ≈ pp.u[end] atol = 1e-5
@@ -91,9 +91,11 @@ function test_flow_function()
         # computing objective
         ψ(t) = 0.5 * (x(t)[1]^2 + x(t)[2]^2 + u(t)^2)
         #
-        o = solve(ODEProblem((_, _, t) -> ψ(t), 0, (t0, tf)), Tsit5(), reltol = 1e-8, abstol = 1e-8)
+        o = solve(
+            ODEProblem((_, _, t) -> ψ(t), 0, (t0, tf)), Tsit5(); reltol=1e-8, abstol=1e-8
+        )
         #
-        f = Flow((t, _) -> ψ(t), autonomous = false, variable = false)
+        f = Flow((t, _) -> ψ(t); autonomous=false, variable=false)
         oo = f((t0, tf), 0)
         #
         Test.@test o.u[end] ≈ oo.u[end] atol = 1e-5
diff --git a/test/test_flow_hamiltonian.jl b/test/test_flow_hamiltonian.jl
index 8dd5b98..9339b41 100644
--- a/test/test_flow_hamiltonian.jl
+++ b/test/test_flow_hamiltonian.jl
@@ -13,7 +13,7 @@ function test_flow_hamiltonian()
 
     @testset "2D non autonomous, variable" begin
         H(t, x, p, l) = p[1] * x[2] + p[2] * p[2] + 0.5 * l * p[2]^2
-        z = Flow(Hamiltonian(H, autonomous = false, variable = true))
+        z = Flow(Hamiltonian(H; autonomous=false, variable=true))
         t0 = 0.0
         tf = 1.0
         x0 = [-1.0, 0.0]
diff --git a/test/test_flow_hamiltonian_vector_field.jl b/test/test_flow_hamiltonian_vector_field.jl
index 96d3e49..f0159a8 100644
--- a/test/test_flow_hamiltonian_vector_field.jl
+++ b/test/test_flow_hamiltonian_vector_field.jl
@@ -13,7 +13,7 @@ function test_flow_hamiltonian_vector_field()
 
     @testset "2D non autonomous, variable" begin
         Hv(t, x, p, l) = [x[2], (2 + l) * p[2]], [0.0, -p[1]]
-        z = Flow(HamiltonianVectorField(Hv, autonomous = false, variable = true))
+        z = Flow(HamiltonianVectorField(Hv; autonomous=false, variable=true))
         t0 = 0.0
         tf = 1.0
         x0 = [-1.0, 0.0]
diff --git a/test/test_flow_vector_field.jl b/test/test_flow_vector_field.jl
index 8cada9b..30f95b7 100644
--- a/test/test_flow_vector_field.jl
+++ b/test/test_flow_vector_field.jl
@@ -12,7 +12,7 @@ function test_flow_vector_field()
 
     @testset "4D non autonomous, variable" begin
         V(t, z, l) = [z[2], (2 + l) * z[2 + 2], 0.0, -z[2 + 1]]
-        z = Flow(VectorField(V, autonomous = false, variable = true))
+        z = Flow(VectorField(V; autonomous=false, variable=true))
         t0 = 0.0
         tf = 1.0
         x0 = [-1.0, 0.0]
diff --git a/test/test_optimal_control_problem.jl b/test/test_optimal_control_problem.jl
index e39375b..7dc1e97 100644
--- a/test/test_optimal_control_problem.jl
+++ b/test/test_optimal_control_problem.jl
@@ -23,7 +23,7 @@ function test_optimal_control_problem()
             ∫(0.5u(t)^2) → min
         end
 
-        f = Flow(ocp, (x, p) -> p[2]; alg = BS5())
+        f = Flow(ocp, (x, p) -> p[2]; alg=BS5())
         p0 = [12, 6]
 
         xf_, pf = f(t0, x0, p0, tf)
@@ -58,8 +58,8 @@ function test_optimal_control_problem()
         p0 = 2tf / 3
 
         #
-        f = Flow(ocp, (x, p, tf) -> tf * p; alg = Tsit5())
-        sol = f((t0, tf), x0, p0, tf; alg = BS5())
+        f = Flow(ocp, (x, p, tf) -> tf * p; alg=Tsit5())
+        sol = f((t0, tf), x0, p0, tf; alg=BS5())
         Test.@test plot(sol) isa Plots.Plot
     end
 
@@ -209,9 +209,9 @@ function test_optimal_control_problem()
         ocp = Model()
         state!(ocp, n)   # dimension of the state
         control!(ocp, m) # dimension of the control
-        time!(ocp, t0 = t0, tf = tf)
-        constraint!(ocp, :initial, lb = x0, ub = x0)
-        constraint!(ocp, :mixed, f = (x, u) -> x + u, lb = -Inf, ub = 0)
+        time!(ocp; t0=t0, tf=tf)
+        constraint!(ocp, :initial; lb=x0, ub=x0)
+        constraint!(ocp, :mixed; f=(x, u) -> x + u, lb=-Inf, ub=0)
         dynamics!(ocp, (x, u) -> u)
         objective!(ocp, :lagrange, (x, u) -> -u)
 
@@ -277,10 +277,10 @@ function test_optimal_control_problem()
         ocp = Model()
         state!(ocp, n)   # dimension of the state
         control!(ocp, m) # dimension of the control
-        time!(ocp, t0 = t0, tf = tf)
-        constraint!(ocp, :initial, lb = x0, ub = x0)
-        constraint!(ocp, :final, lb = xf, ub = xf)
-        constraint!(ocp, :state, rg = Index(1), lb = -Inf, ub = l)
+        time!(ocp; t0=t0, tf=tf)
+        constraint!(ocp, :initial; lb=x0, ub=x0)
+        constraint!(ocp, :final; lb=xf, ub=xf)
+        constraint!(ocp, :state; rg=Index(1), lb=-Inf, ub=l)
         A = [
             0 1
             0 0