diff --git a/Project.toml b/Project.toml
index c6f11416..864ef6d1 100644
--- a/Project.toml
+++ b/Project.toml
@@ -10,6 +10,7 @@ DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 HSL = "34c5aeac-e683-54a6-a0e9-6e0fdc586c50"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 NLPModels = "a4795742-8479-5a88-8948-cc11e1c8c1a6"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [weakdeps]
 MadNLP = "2621e9c9-9eb4-46b1-8089-e8c72242dfb6"
@@ -27,4 +28,5 @@ HSL = "0.5"
 MadNLP = "0.8"
 NLPModels = "0.21"
 NLPModelsIpopt = "0.10"
+SparseArrays = "1.10.0"
 julia = "1.10"
diff --git a/docs/src/index.md b/docs/src/index.md
index 81ad3d61..a53ef00b 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -82,15 +82,11 @@ LB \le C(X) \le UB
 \right.
 ```
 
-We use packages from [JuliaSmoothOptimizers](https://github.com/JuliaSmoothOptimizers) to solve the (NLP) problem.
+Solving the (NLP) problem is done using packages from [JuliaSmoothOptimizers](https://github.com/JuliaSmoothOptimizers), with Ipopt as the default solver.
 
-As input of this package we use an [`OptimalControlModel`](@ref) structure from CTBase.
+On the input side of this package, we use an [`OptimalControlModel`](@ref) structure from CTBase to define the (OCP).
 
-!!! note "Current limitations"
-
-    The current implemented is limited to
-    - trapezoidal rule for the ODE discretization
-    - `Ipopt` for the optimization software
+The direct transcription to build the (NLP) can use discretization schemes such as trapeze (default), midpoint, or Gauss-Legendre collocations.
 
 !!! note "Related packages"
 
diff --git a/ext/CTSolveExtIpopt.jl b/ext/CTSolveExtIpopt.jl
index 5e61d4a6..e13aad6a 100644
--- a/ext/CTSolveExtIpopt.jl
+++ b/ext/CTSolveExtIpopt.jl
@@ -59,6 +59,7 @@ function CTDirect.solve_docp(
         tol = tol,
         max_iter = max_iter,
         sb = "yes",
+        #check_derivatives_for_naninf = "yes",
         linear_solver = linear_solver;
         kwargs...,
     )
diff --git a/src/CTDirect.jl b/src/CTDirect.jl
index a30a4e80..9a30661b 100644
--- a/src/CTDirect.jl
+++ b/src/CTDirect.jl
@@ -5,6 +5,7 @@ using DocStringExtensions
 using ADNLPModels               # docp model with AD
 using LinearAlgebra             # norm and misc
 using HSL
+using SparseArrays
 
 import CTBase: OptimalControlSolution, CTBase   # extended
 
diff --git a/src/default.jl b/src/default.jl
index b2d317ae..39058308 100644
--- a/src/default.jl
+++ b/src/default.jl
@@ -24,6 +24,14 @@ The default value is `nothing`.
 """
 __time_grid() = nothing
 
+"""
+$(TYPEDSIGNATURES)
+
+Used to set the default control type for IRK schemes 
+The default value is `:constant`.
+"""
+__control_type() = :constant
+
 """
 $(TYPEDSIGNATURES)
 Used to set the default backend for AD in ADNLPModels.
diff --git a/src/disc/irk.jl b/src/disc/irk.jl
index 0bb9643c..981aa914 100644
--- a/src/disc/irk.jl
+++ b/src/disc/irk.jl
@@ -5,9 +5,10 @@ Internal layout for NLP variables:
  .., 
  X_N-1, U_N-1, K_N-1^1..K_N-1^s,
  X_N, U_N, V]
-with s the stage number and U given by either linear interpolation in [t_i, t_i+1]
-or constant interpolation for 1-stage methods or if specfied (U_N might end up unused)
-Path constraints are all evaluated at time steps
+with s the stage number and U piecewise constant equal to U_i in [t_i, t_i+1]
+or, for methods with s>1, piecewise linear if option control_type set to :linear
+NB. U_N may be removed at some point if we disable piecewise linear control
+Path constraints are all evaluated at time steps, including final time.
 =#
 
 
@@ -57,9 +58,9 @@ struct Gauss_Legendre_2 <: GenericIRK
     _step_variables_block::Int
     _state_stage_eqs_block::Int
     _step_pathcons_block::Int
-    _constant_control::Bool
+    _control_type::Symbol
 
-    function Gauss_Legendre_2(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons, constant_control)
+    function Gauss_Legendre_2(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons, control_type)
         
         stage = 2
 
@@ -70,7 +71,7 @@ struct Gauss_Legendre_2 <: GenericIRK
         [0.5, 0.5],
         [(0.5 - sqrt(3) / 6), (0.5 + sqrt(3) / 6)],
         step_variables_block, state_stage_eqs_block, step_pathcons_block,
-        constant_control
+        control_type
         )
 
         return disc, dim_NLP_variables, dim_NLP_constraints
@@ -93,9 +94,9 @@ struct Gauss_Legendre_3 <: GenericIRK
     _step_variables_block::Int
     _state_stage_eqs_block::Int
     _step_pathcons_block::Int
-    _constant_control::Bool
+    _control_type::Symbol
 
-    function Gauss_Legendre_3(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons, constant_control)
+    function Gauss_Legendre_3(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons, control_type)
         
         stage = 3
 
@@ -107,7 +108,7 @@ struct Gauss_Legendre_3 <: GenericIRK
         (5/36 + sqrt(15) / 30) (2/9 + sqrt(15) / 15) (5.0/36.0)],
         [5.0/18.0, 4.0/9.0, 5.0/18.0],
         [0.5 - 0.1*sqrt(15), 0.5, 0.5 + 0.1*sqrt(15)],
-        step_variables_block, state_stage_eqs_block, step_pathcons_block, constant_control
+        step_variables_block, state_stage_eqs_block, step_pathcons_block, control_type
         )
 
         return disc, dim_NLP_variables, dim_NLP_constraints
@@ -183,7 +184,7 @@ function get_OCP_control_at_time_step(xu, docp::DOCP{ <: GenericIRK, <: ScalVect
     return @view xu[(offset + 1):(offset + docp.dim_NLP_u)]
 end
 function get_OCP_control_at_time_stage(xu, docp::DOCP{ <: GenericIRK, <: ScalVect, ScalVariable, <: ScalVect}, i, cj)
-    if (docp.discretization.stage == 1) || (docp.discretization._constant_control)
+    if (docp.discretization.stage == 1) || (docp.discretization._control_type == :constant)
         # constant interpolation on step
         return get_OCP_control_at_time_step(xu, docp, i)
     else
@@ -194,7 +195,7 @@ function get_OCP_control_at_time_stage(xu, docp::DOCP{ <: GenericIRK, <: ScalVec
     end
 end
 function get_OCP_control_at_time_stage(xu, docp::DOCP{ <: GenericIRK, <: ScalVect, VectVariable, <: ScalVect}, i, cj)
-    if (docp.discretization.stage == 1) || (docp.discretization._constant_control)
+    if (docp.discretization.stage == 1) || (docp.discretization._control_type == :constant)
         # constant interpolation on step
         return get_OCP_control_at_time_step(xu, docp, i)
     else
@@ -251,8 +252,8 @@ $(TYPEDSIGNATURES)
 Set work array for all dynamics and lagrange cost evaluations
 """
 function setWorkArray(docp::DOCP{ <: GenericIRK}, xu, time_grid, v)
-    # work array layout: [x_ij ; sum_bk ; u_ij] ?    
-    work = similar(xu, docp.dim_OCP_x + docp.dim_NLP_x + docp.dim_NLP_u)
+    # work array layout: [x_ij ; sum_bk]   
+    work = similar(xu, docp.dim_OCP_x + docp.dim_NLP_x)
     return work
 end
 
@@ -264,12 +265,9 @@ Convention: 1 <= i <= dim_NLP_steps (+1)
 """
 function setStepConstraints!(docp::DOCP{ <: GenericIRK}, c, xu, v, time_grid, i, work)
 
-    # work array layout: [x_ij ; sum_bk ; u_ij] ?
+    # work array layout: [x_ij ; sum_bk]
     work_xij = @view work[1:docp.dim_OCP_x]
     work_sumbk = @view work[docp.dim_OCP_x+1:docp.dim_OCP_x+docp.dim_NLP_x]
-    #work_sumbk .= zero(eltype(xu)) AD bug when affecting constant values...
-    @views @. work_sumbk[1:docp.dim_NLP_x] = xu[1:docp.dim_NLP_x] * 0.
-    #work_uij ?
 
     # offset for previous steps
     offset = (i-1)*(docp.discretization._state_stage_eqs_block + docp.discretization._step_pathcons_block)
@@ -297,7 +295,11 @@ function setStepConstraints!(docp::DOCP{ <: GenericIRK}, c, xu, v, time_grid, i,
             kij = get_stagevars_at_time_step(xu, docp, i, j)
             
             # update sum b_j k_i^j (w/ lagrange term) for state equation after loop
-            @views @. work_sumbk[1:docp.dim_NLP_x] = work_sumbk[1:docp.dim_NLP_x] + docp.discretization.butcher_b[j] * kij[1:docp.dim_NLP_x]
+            if j == 1
+                @views @. work_sumbk[1:docp.dim_NLP_x] = docp.discretization.butcher_b[j] * kij[1:docp.dim_NLP_x]
+            else
+                @views @. work_sumbk[1:docp.dim_NLP_x] = work_sumbk[1:docp.dim_NLP_x] + docp.discretization.butcher_b[j] * kij[1:docp.dim_NLP_x]
+            end
 
             # state at stage: x_i^j = x_i + h_i sum a_jl k_i^l
             # +++ still some allocations here
@@ -312,8 +314,7 @@ function setStepConstraints!(docp::DOCP{ <: GenericIRK}, c, xu, v, time_grid, i,
                 xij = work_xij
             end
 
-            # control at stage: interpolation between u_i and u_i+1 
-            # +++ use work aray to reduce allocs ?
+            # control at stage
             uij = get_OCP_control_at_time_stage(xu, docp, i, cj)
 
             # stage equations k_i^j = f(t_i^j, x_i^j, u_i, v) as c[] = k - f
@@ -342,3 +343,199 @@ function setStepConstraints!(docp::DOCP{ <: GenericIRK}, c, xu, v, time_grid, i,
     setPathConstraints!(docp, c, ti, xi, ui, v, offset)
 
 end
+
+
+"""
+$(TYPEDSIGNATURES)
+
+Build sparsity pattern for Jacobian of constraints
+"""
+function DOCP_Jacobian_pattern(docp::DOCP{ <: GenericIRK})
+
+    if docp.discretization._control_type != :constant
+        error("Manual Jacobian sparsity pattern not supported for IRK scheme with piecewise linear control")
+    end
+
+    # vector format for sparse matrix
+    Is = Vector{Int}(undef, 0)
+    Js = Vector{Int}(undef, 0)
+
+    s = docp.discretization.stage
+
+    # index alias for v
+    v_start = docp.dim_NLP_variables - docp.dim_NLP_v + 1
+    v_end = docp.dim_NLP_variables
+
+    # 1. main loop over steps
+    for i = 1:docp.dim_NLP_steps
+
+        # constraints block and offset: state equation, stage equations, path constraints
+        c_block = docp.discretization._state_stage_eqs_block + docp.discretization._step_pathcons_block
+        c_offset = (i-1)*c_block
+
+        # contiguous variables blocks will be used when possible
+        # x_i (l_i) u_i k_i x_i+1 (l_i+1)
+        var_offset = (i-1)*docp.discretization._step_variables_block
+        xi_start = var_offset + 1
+        xi_end = var_offset + docp.dim_OCP_x
+        ui_start = var_offset + docp.dim_NLP_x + 1
+        ui_end = var_offset + docp.dim_NLP_x + docp.dim_NLP_u
+        ki_start = var_offset + docp.dim_NLP_x + docp.dim_NLP_u + 1
+        ki_end = var_offset + docp.discretization._step_variables_block
+        xip1_end = var_offset + docp.discretization._step_variables_block + docp.dim_OCP_x
+        li = var_offset + docp.dim_NLP_x
+        lip1 = var_offset + docp.discretization._step_variables_block + docp.dim_NLP_x
+
+        # 1.1 state eq 0 = x_i+1 - (x_i + h_i sum_j b_j k_ij)
+        # depends on x_i, k_ij, x_i+1, and v for h_i in variable times case !
+        # skip l_i, u_i; should skip k_i[n+1] also but annoying...
+        add_nonzero_block!(Is, Js, c_offset+1, c_offset+docp.dim_OCP_x, xi_start, xi_end)
+        add_nonzero_block!(Is, Js, c_offset+1, c_offset+docp.dim_OCP_x, ki_start, xip1_end)
+        add_nonzero_block!(Is, Js, c_offset+1, c_offset+docp.dim_OCP_x, v_start, v_end)
+        # 1.2 lagrange part l_i+1 = l_i + h_i (sum_j b_j k_ij)[n+1]
+        # depends on l_i, k_ij[n+1], l_i+1, and v for h_i in variable times case !
+        if docp.is_lagrange
+            add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x, li)
+            add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x, lip1)
+            for i=1:s
+                kij_l = var_offset + docp.dim_NLP_x + docp.dim_NLP_u + i*docp.dim_NLP_x 
+                add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x, kij_l)
+            end
+            add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x, c_offset+docp.dim_NLP_x, v_start, v_end)
+        end
+
+        # 1.3 stage equations k_ij = f(t_ij, x_ij, u_ij, v) (with lagrange part)
+        # with x_ij = x_i + sum_l a_il k_jl and assuming u_ij = u_i
+        # depends on x_i, u_i, k_i, and v; skip l_i (could skip k_ij[n+1] too...)
+        add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x+1, c_offset+(s+1)*docp.dim_NLP_x, xi_start, xi_end)
+        add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x+1, c_offset+(s+1)*docp.dim_NLP_x, ui_start, ki_end)      
+        add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x+1, c_offset+(s+1)*docp.dim_NLP_x, v_start, v_end)
+
+        # 1.4 path constraint g(t_i, x_i, u_i, v)
+        # depends on x_i, u_i, v; skip l_i
+        add_nonzero_block!(Is, Js, c_offset+(s+1)*docp.dim_NLP_x+1, c_offset+c_block, xi_start, xi_end)
+        add_nonzero_block!(Is, Js, c_offset+(s+1)*docp.dim_NLP_x+1, c_offset+c_block, ui_start, ui_end)
+        add_nonzero_block!(Is, Js, c_offset+(s+1)*docp.dim_NLP_x+1, c_offset+c_block, v_start, v_end)
+    end
+
+    # 2. final path constraints (xf, uf, v)
+    c_offset = docp.dim_NLP_steps * (docp.discretization._state_stage_eqs_block + docp.discretization._step_pathcons_block)
+    c_block = docp.discretization._step_pathcons_block
+    var_offset = docp.dim_NLP_steps*docp.discretization._step_variables_block
+    xf_start = var_offset + 1
+    xf_end = var_offset + docp.dim_OCP_x
+    uf_start = var_offset + docp.dim_NLP_x + 1
+    uf_end = var_offset + docp.dim_NLP_x + docp.dim_NLP_u
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, xf_start, xf_end)
+    add_nonzero_block!(Is, Js, c_offset+1,c_offset+c_block, uf_start, uf_end)
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, v_start, v_end)
+
+    # 3. boundary constraints (x0, xf, v)
+    c_offset = docp.dim_NLP_steps * (docp.discretization._state_stage_eqs_block + docp.discretization._step_pathcons_block) + docp.discretization._step_pathcons_block
+    c_block = docp.dim_boundary_cons + docp.dim_v_cons
+    x0_start = 1
+    x0_end = docp.dim_OCP_x
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, x0_start, x0_end)
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, xf_start, xf_end)
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, v_start, v_end)
+    # 3.4 null initial condition for lagrangian cost state l0
+    if docp.is_lagrange
+        add_nonzero_block!(Is, Js, docp.dim_NLP_constraints, docp.dim_NLP_x)
+    end
+
+    # build and return sparse matrix
+    nnzj = length(Is)
+    Vs = ones(Bool, nnzj)
+    return sparse(Is, Js, Vs, docp.dim_NLP_constraints, docp.dim_NLP_variables)
+end
+
+
+"""
+$(TYPEDSIGNATURES)
+
+Build sparsity pattern for Hessian of Lagrangian
+"""
+function DOCP_Hessian_pattern(docp::DOCP{ <: GenericIRK})
+
+    if docp.discretization._control_type != :constant
+        error("Manual Hessian sparsity pattern not supported for IRK scheme with piecewise linear control")
+    end
+
+    # NB. need to provide full pattern for coloring, not just upper/lower part
+    Is = Vector{Int}(undef, 0)
+    Js = Vector{Int}(undef, 0)
+   
+    s = docp.discretization.stage
+
+    # index alias for v
+    v_start = docp.dim_NLP_variables - docp.dim_NLP_v + 1
+    v_end = docp.dim_NLP_variables
+
+    # 0. objective
+    # 0.1 mayer cost (x0, xf, v) 
+    # -> grouped with term 3. for boundary conditions
+    # 0.2 lagrange case (lf)
+    # -> 2nd order term is zero
+   
+    # 1. main loop over steps
+    # 1.0 v / v term
+    add_nonzero_block!(Is, Js, v_start, v_end, v_start, v_end)
+
+    for i = 1:docp.dim_NLP_steps
+
+        # contiguous variables blocks will be used when possible
+        # x_i (l_i) u_i k_i x_i+1 (l_i+1)
+        var_offset = (i-1)*docp.discretization._step_variables_block
+        xi_start = var_offset + 1
+        xi_end = var_offset + docp.dim_OCP_x
+        ui_start = var_offset + docp.dim_NLP_x + 1
+        ui_end = var_offset + docp.dim_NLP_x + docp.dim_NLP_u
+        ki_start = var_offset + docp.dim_NLP_x + docp.dim_NLP_u + 1
+        ki_end = var_offset + (s+1)*docp.dim_NLP_x + docp.dim_NLP_u
+
+        # 1.1 state eq 0 = x_i+1 - (x_i + h_i sum_j b_j k_ij)
+        # -> 2nd order terms are zero
+        # 1.2 lagrange part 0 = l_i+1 - (l_i + h_i (sum_j b_j k_ij[n+1]))
+        # -> 2nd order terms are zero
+
+        # 1.3 stage equations 0 = k_ij - f(t_ij, x_ij, u_ij, v) (with lagrange part)
+        # with x_ij = x_i + sum_l a_il k_jl and assuming u_ij = u_i
+        # depends on x_i, u_i, k_i, and v; skip l_i (could skip k_ij[n+1] too...)
+        add_nonzero_block!(Is, Js, xi_start, xi_end, xi_start, xi_end)
+        add_nonzero_block!(Is, Js, ui_start, ki_end, ui_start, ki_end)
+        add_nonzero_block!(Is, Js, xi_start, xi_end, ui_start, ki_end; sym=true)
+        add_nonzero_block!(Is, Js, xi_start, xi_end, v_start, v_end; sym=true)
+        add_nonzero_block!(Is, Js, ui_start, ki_end, v_start, v_end; sym=true)
+
+        # 1.4 path constraint g(t_i, x_i, u_i, v)
+        # -> included in 1.3
+    end
+
+    # 2. final path constraints (xf, uf, v) (assume present)
+    var_offset = docp.dim_NLP_steps*docp.discretization._step_variables_block
+    xf_start = var_offset + 1
+    xf_end = var_offset + docp.dim_OCP_x
+    uf_start = var_offset + docp.dim_NLP_x + 1
+    uf_end = var_offset + docp.dim_NLP_x + docp.dim_NLP_u
+    add_nonzero_block!(Is, Js, xf_start, xf_end, xf_start, xf_end)
+    add_nonzero_block!(Is, Js, uf_start, uf_end, uf_start, uf_end)
+    add_nonzero_block!(Is, Js, xf_start, xf_end, uf_start, uf_end; sym=true) 
+    add_nonzero_block!(Is, Js, xf_start, xf_end, v_start, v_end; sym=true)
+    add_nonzero_block!(Is, Js, uf_start, uf_end, v_start, v_end; sym=true)
+
+    # 3. boundary constraints (x0, xf, v) or mayer cost g0(x0, xf, v) (assume present)
+    # -> x0 / x0, x0 / v terms included in first loop iteration
+    # -> xf / xf, xf / v terms included in 2.
+    x0_start = 1
+    x0_end = docp.dim_OCP_x
+    add_nonzero_block!(Is, Js, x0_start, x0_end, xf_start, xf_end; sym=true)
+
+    # 3.1 null initial condition for lagrangian cost state l0
+    # -> 2nd order term is zero
+   
+    # build and return sparse matrix
+    nnzj = length(Is)
+    Vs = ones(Bool, nnzj)
+    return sparse(Is, Js, Vs, docp.dim_NLP_variables, docp.dim_NLP_variables)
+
+end
diff --git a/src/disc/midpoint.jl b/src/disc/midpoint.jl
index 1dc8faee..c2de23e1 100644
--- a/src/disc/midpoint.jl
+++ b/src/disc/midpoint.jl
@@ -19,13 +19,11 @@ struct Midpoint <: Discretization
         # aux variables
         step_variables_block = dim_NLP_x * 2 + dim_NLP_u
         state_stage_eqs_block = dim_NLP_x * 2
+        step_pathcons_block = dim_u_cons + dim_x_cons + dim_xu_cons
 
         # NLP variables size ([state, control]_1..N, final state, variable)
         dim_NLP_variables = dim_NLP_steps * step_variables_block + dim_NLP_x + dim_NLP_v
         
-        # Path constraints (control, state, mixed)
-        step_pathcons_block = dim_u_cons + dim_x_cons + dim_xu_cons
-
         # NLP constraints size ([dynamics, path]_1..N, final path, boundary, variable)
         dim_NLP_constraints = dim_NLP_steps * (state_stage_eqs_block + step_pathcons_block) + step_pathcons_block + dim_boundary_cons + dim_v_cons
 
@@ -188,3 +186,174 @@ function setStepConstraints!(docp::DOCP{Midpoint}, c, xu, v, time_grid, i, work)
 end
 
 
+"""
+$(TYPEDSIGNATURES)
+
+Build sparsity pattern for Jacobian of constraints
+"""
+function DOCP_Jacobian_pattern(docp::DOCP{Midpoint})
+
+    # vector format for sparse matrix
+    Is = Vector{Int}(undef, 0)
+    Js = Vector{Int}(undef, 0)
+
+    # index alias for v
+    v_start = docp.dim_NLP_variables - docp.dim_NLP_v + 1
+    v_end = docp.dim_NLP_variables
+
+    # 1. main loop over steps
+    for i = 1:docp.dim_NLP_steps
+
+        # constraints block and offset: state equation, stage equation, path constraints
+        c_block = docp.discretization._state_stage_eqs_block + docp.discretization._step_pathcons_block
+        c_offset = (i-1)*c_block
+
+        # contiguous variables blocks will be used when possible
+        # x_i (l_i) u_i k_i x_i+1 (l_i+1)
+        var_offset = (i-1)*docp.discretization._step_variables_block
+        xi_start = var_offset + 1
+        xi_end = var_offset + docp.dim_OCP_x
+        ui_start = var_offset + docp.dim_NLP_x + 1
+        ui_end = var_offset + docp.dim_NLP_x + docp.dim_NLP_u
+        ki_start = var_offset + docp.dim_NLP_x + docp.dim_NLP_u + 1
+        xip1_end = var_offset + docp.discretization._step_variables_block + docp.dim_OCP_x
+        li = var_offset + docp.dim_NLP_x
+        lip1 = var_offset + docp.discretization._step_variables_block + docp.dim_NLP_x
+
+        # 1.1 state eq 0 = x_i+1 - (x_i + h_i * k_i)
+        # depends on x_i, k_i, x_i+1, and v for h_i in variable times case !
+        # skip l_i, u_i; should skip k_i[n+1] also but annoying...
+        add_nonzero_block!(Is, Js, c_offset+1, c_offset+docp.dim_OCP_x, xi_start, xi_end)
+        add_nonzero_block!(Is, Js, c_offset+1, c_offset+docp.dim_OCP_x, ki_start, xip1_end)
+        add_nonzero_block!(Is, Js, c_offset+1, c_offset+docp.dim_OCP_x, v_start, v_end)
+        # 1.2 lagrange part 0 = l_i+1 - (l_i + h_i * k_i[n+1])
+        # depends on l_i, k_i[n+1], l_i+1, and v for h_i in variable times case !
+        if docp.is_lagrange
+            add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x, li)
+            add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x, lip1)
+            ki_l = var_offset + docp.dim_NLP_x + docp.dim_NLP_u + docp.dim_NLP_x 
+            add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x, ki_l)
+            add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x, c_offset+docp.dim_NLP_x, v_start, v_end)
+        end
+
+        # 1.3 stage equation 0 = k_i - f(t_s, x_s, u_i, v)
+        # with t_s = (t_i + t_i+1)/2    x_s = (x_i + x_i+1)/2
+        # depends on x_i, u_i, x_i+1, k_i, and v; skip l_i
+        add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x+1, c_offset+2*docp.dim_NLP_x, xi_start, xi_end)
+        add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x+1, c_offset+2*docp.dim_NLP_x, ui_start, xip1_end)      
+        add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x+1, c_offset+2*docp.dim_NLP_x, v_start, v_end)
+
+        # 1.4 path constraint g(t_i, x_i, u_i, v)
+        # depends on x_i, u_i, v; skip l_i
+        add_nonzero_block!(Is, Js, c_offset+2*docp.dim_NLP_x+1, c_offset+c_block, xi_start, xi_end)
+        add_nonzero_block!(Is, Js, c_offset+2*docp.dim_NLP_x+1, c_offset+c_block, ui_start, ui_end)
+        add_nonzero_block!(Is, Js, c_offset+2*docp.dim_NLP_x+1, c_offset+c_block, v_start, v_end)
+    end
+
+    # 2. final path constraints (xf, uf, v)
+    c_offset = docp.dim_NLP_steps * (docp.discretization._state_stage_eqs_block + docp.discretization._step_pathcons_block)
+    c_block = docp.discretization._step_pathcons_block
+    var_offset = docp.dim_NLP_steps*docp.discretization._step_variables_block
+    xf_start = var_offset + 1
+    xf_end = var_offset + docp.dim_OCP_x
+    uf_start = var_offset-docp.discretization._step_variables_block + docp.dim_NLP_x + 1
+    uf_end = var_offset-docp.discretization._step_variables_block + docp.dim_NLP_x + docp.dim_NLP_u
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, xf_start, xf_end)
+    add_nonzero_block!(Is, Js, c_offset+1,c_offset+c_block, uf_start, uf_end)
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, v_start, v_end)
+
+    # 3. boundary constraints (x0, xf, v)
+    c_offset = docp.dim_NLP_steps * (docp.discretization._state_stage_eqs_block + docp.discretization._step_pathcons_block) + docp.discretization._step_pathcons_block
+    c_block = docp.dim_boundary_cons + docp.dim_v_cons
+    x0_start = 1
+    x0_end = docp.dim_OCP_x
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, x0_start, x0_end)
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, xf_start, xf_end)
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, v_start, v_end)
+    # 3.4 null initial condition for lagrangian cost state l0
+    if docp.is_lagrange
+        add_nonzero_block!(Is, Js, docp.dim_NLP_constraints, docp.dim_NLP_x)
+    end
+
+    # build and return sparse matrix
+    nnzj = length(Is)
+    Vs = ones(Bool, nnzj)
+    return sparse(Is, Js, Vs, docp.dim_NLP_constraints, docp.dim_NLP_variables)
+end
+
+
+"""
+$(TYPEDSIGNATURES)
+
+Build sparsity pattern for Hessian of Lagrangian
+"""
+function DOCP_Hessian_pattern(docp::DOCP{Midpoint})
+
+    # NB. need to provide full pattern for coloring, not just upper/lower part
+    Is = Vector{Int}(undef, 0)
+    Js = Vector{Int}(undef, 0)
+
+    # index alias for v
+    v_start = docp.dim_NLP_variables - docp.dim_NLP_v + 1
+    v_end = docp.dim_NLP_variables
+
+    # 0. objective
+    # 0.1 mayer cost (x0, xf, v) 
+    # -> grouped with term 3. for boundary conditions
+    # 0.2 lagrange case (lf)
+    # -> 2nd order term is zero
+   
+    # 1. main loop over steps
+    # 1.0 v / v term
+    add_nonzero_block!(Is, Js, v_start, v_end, v_start, v_end)
+
+    for i = 1:docp.dim_NLP_steps
+
+        # contiguous variables blocks will be used when possible
+        # x_i (l_i) u_i k_i x_i+1 (l_i+1)
+        var_offset = (i-1)*docp.discretization._step_variables_block
+        xi_start = var_offset + 1
+        xi_end = var_offset + docp.dim_OCP_x
+        xip1_end = var_offset + docp.discretization._step_variables_block + docp.dim_NLP_x
+        ui_start = var_offset + docp.dim_NLP_x + 1
+
+        # 1.1 state eq 0 = x_i+1 - (x_i + h_i * k_i)
+        # -> 2nd order terms are zero
+        # 1.2 lagrange part 0 = l_i+1 - (l_i + h_i * k_i[n+1])
+        # -> 2nd order terms are zero
+
+        # 1.3 stage equations 0 = k_i - f(t_s, x_s, u_i, v)
+        # with t_s = (t_i + t_i+1)/2    x_s = (x_i + x_i+1)/2
+        # depends on x_i, u_i, k_i, x_i+1, and v; skip l_i    
+        add_nonzero_block!(Is, Js, xi_start, xi_end, xi_start, xi_end)
+        add_nonzero_block!(Is, Js, ui_start, xip1_end, ui_start, xip1_end)
+        add_nonzero_block!(Is, Js, xi_start, xi_end, ui_start, xip1_end; sym=true)
+        add_nonzero_block!(Is, Js, xi_start, xi_end, v_start, v_end; sym=true)
+        add_nonzero_block!(Is, Js, ui_start, xip1_end, v_start, v_end; sym=true)
+
+        # 1.4 path constraint g(t_i, x_i, u_i, v)
+        # -> included in 1.3
+    end
+
+    # 2. final path constraints (xf, uf, v)
+    # -> included in last loop iteration (with x_i+1 as x_j and u_i as u_f)
+
+    # 3. boundary constraints (x0, xf, v) or mayer cost g0(x0, xf, v) (assume present)
+    # -> x0 / x0, x0 / v terms included in first loop iteration
+    # -> xf / xf, xf / v terms included in last loop iteration (with x_i+1 as x_f)
+    x0_start = 1
+    x0_end = docp.dim_OCP_x
+    var_offset = docp.dim_NLP_steps*docp.discretization._step_variables_block
+    xf_start = var_offset + 1
+    xf_end = var_offset + docp.dim_OCP_x
+    add_nonzero_block!(Is, Js, x0_start, x0_end, xf_start, xf_end; sym=true)
+
+    # 3.1 null initial condition for lagrangian cost state l0
+    # -> 2nd order term is zero
+   
+    # build and return sparse matrix
+    nnzj = length(Is)
+    Vs = ones(Bool, nnzj)
+    return sparse(Is, Js, Vs, docp.dim_NLP_variables, docp.dim_NLP_variables)
+
+end
diff --git a/src/disc/trapeze.jl b/src/disc/trapeze.jl
index 1f8b9e8a..84055748 100644
--- a/src/disc/trapeze.jl
+++ b/src/disc/trapeze.jl
@@ -7,23 +7,25 @@ Internal layout for NLP variables:
 struct Trapeze <: Discretization
 
     info::String
-    _step_pathcons_block::Int
+    _step_variables_block::Int
     _state_stage_eqs_block::Int
+    _step_pathcons_block::Int
 
     # constructor
     function Trapeze(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons)
 
-        # NLP variables size ([state, control]_1..N+1, variable)
-        dim_NLP_variables = (dim_NLP_steps + 1) * (dim_NLP_x + dim_NLP_u) + dim_NLP_v
-
-        # Path constraints (control, state, mixed)
+        # aux variables
+        step_variables_block = dim_NLP_x + dim_NLP_u
         state_stage_eqs_block = dim_NLP_x
         step_pathcons_block = dim_u_cons + dim_x_cons + dim_xu_cons
 
+        # NLP variables size ([state, control]_1..N+1, variable)
+        dim_NLP_variables = (dim_NLP_steps + 1) * step_variables_block + dim_NLP_v
+
         # NLP constraints size ([dynamics, stage, path]_1..N, final path, boundary, variable)
         dim_NLP_constraints = dim_NLP_steps * (state_stage_eqs_block + step_pathcons_block) + step_pathcons_block + dim_boundary_cons + dim_v_cons
 
-        disc = new("Implicit Trapeze aka Crank-Nicolson, 2nd order, A-stable", step_pathcons_block, state_stage_eqs_block)
+        disc = new("Implicit Trapeze aka Crank-Nicolson, 2nd order, A-stable", step_variables_block, state_stage_eqs_block, step_pathcons_block)
 
         return disc, dim_NLP_variables, dim_NLP_constraints
     end
@@ -171,3 +173,197 @@ function setStepConstraints!(docp::DOCP{Trapeze}, c, xu, v, time_grid, i, work)
     setPathConstraints!(docp, c, ti, xi, ui, v, offset)
 
 end
+
+
+"""
+$(TYPEDSIGNATURES)
+
+Build sparsity pattern for Jacobian of constraints
+"""
+function DOCP_Jacobian_pattern(docp::DOCP{Trapeze})
+
+    # vector format for sparse matrix
+    Is = Vector{Int}(undef, 0)
+    Js = Vector{Int}(undef, 0)
+
+    # index alias for v
+    v_start = docp.dim_NLP_variables - docp.dim_NLP_v + 1
+    v_end = docp.dim_NLP_variables
+
+    # 1. main loop over steps
+    for i = 1:docp.dim_NLP_steps
+
+        # constraints block and offset: state equation, path constraints
+        c_block = docp.discretization._state_stage_eqs_block + docp.discretization._step_pathcons_block
+        c_offset = (i-1)*c_block
+
+        # contiguous variables blocks will be used when possible 
+        # x_i (l_i) u_i x_i+1 (l_i+1) u_i+1
+        var_block = docp.discretization._step_variables_block * 2
+        var_offset = (i-1)*docp.discretization._step_variables_block
+        xi_start = var_offset + 1
+        xi_end = var_offset + docp.dim_OCP_x
+        ui_start = var_offset + docp.dim_NLP_x + 1
+        ui_end = var_offset + docp.dim_NLP_x + docp.dim_NLP_u
+        xip1_end = var_offset + docp.dim_NLP_x + docp.dim_NLP_u + docp.dim_OCP_x
+        uip1_start = var_offset + docp.dim_NLP_x*2 + docp.dim_NLP_u + 1
+        uip1_end = var_offset + docp.dim_NLP_x*2 + docp.dim_NLP_u*2
+
+        # 1.1 state eq 0 = x_i+1 - (x_i + h_i/2 (f(t_i,x_i,u_i,v) + f(t_i+1,x_i+1,u_i+1,v)))
+        # depends on x_i, u_i, x_i+1, u_i+1; skip l_i, l_i+1; v cf 1.4
+        add_nonzero_block!(Is, Js, c_offset+1, c_offset+docp.dim_OCP_x, xi_start, xi_end)
+        add_nonzero_block!(Is, Js, c_offset+1, c_offset+docp.dim_OCP_x, ui_start, xip1_end)
+        add_nonzero_block!(Is, Js, c_offset+1, c_offset+docp.dim_OCP_x, uip1_start, uip1_end)
+        # 1.2 lagrange part 0 = l_i+1 - (l_i + h_i/2 (l(t_i,x_i,u_i,v) + l(t_i+1,x_i+1,u_i+1,v)))
+        # depends on x_i, l_i, u_i, x_i+1, l_i+1, u_i+1 ie whole variable block; v cf 1.4
+        if docp.is_lagrange
+            add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x, c_offset+docp.dim_NLP_x, var_offset+1, var_offset+var_block)
+        end
+
+        # 1.3 path constraint g(t_i, x_i, u_i, v) 
+        # depends on x_i, u_i; skip l_i; v cf 1.4
+        add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x+1, c_offset+c_block, xi_start, xi_end)
+        add_nonzero_block!(Is, Js, c_offset+docp.dim_NLP_x+1, c_offset+c_block, ui_start, ui_end)
+
+        # 1.4 whole constraint block depends on v
+        add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, v_start, v_end)
+    end
+
+    # 2. final path constraints (xf, uf, v)
+    c_offset = docp.dim_NLP_steps*(docp.discretization._state_stage_eqs_block + docp.discretization._step_pathcons_block)
+    c_block = docp.discretization._step_pathcons_block
+    var_offset = docp.dim_NLP_steps*docp.discretization._step_variables_block
+    xf_start = var_offset + 1
+    xf_end = var_offset + docp.dim_OCP_x
+    uf_start = var_offset + docp.dim_NLP_x + 1
+    uf_end = var_offset + docp.dim_NLP_x + docp.dim_NLP_u
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, xf_start, xf_end)
+    add_nonzero_block!(Is, Js, c_offset+1,c_offset+c_block, uf_start, uf_end)
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, v_start, v_end)
+
+    # 3. boundary constraints (x0, xf, v)
+    c_offset = docp.dim_NLP_steps * (docp.discretization._state_stage_eqs_block + docp.discretization._step_pathcons_block) + docp.discretization._step_pathcons_block
+    c_block = docp.dim_boundary_cons + docp.dim_v_cons
+    x0_start = 1
+    x0_end = docp.dim_OCP_x
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, x0_start, x0_end)
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, xf_start, xf_end)
+    add_nonzero_block!(Is, Js, c_offset+1, c_offset+c_block, v_start, v_end)
+    # 3.4 null initial condition for lagrangian cost state l0
+    if docp.is_lagrange
+        add_nonzero_block!(Is, Js, docp.dim_NLP_constraints, docp.dim_NLP_x)
+    end
+
+    # build and return sparse matrix
+    nnzj = length(Is)
+    Vs = ones(Bool, nnzj)
+    return sparse(Is, Js, Vs, docp.dim_NLP_constraints, docp.dim_NLP_variables)
+end
+
+
+"""
+$(TYPEDSIGNATURES)
+
+Build sparsity pattern for Hessian of Lagrangian
+"""
+function DOCP_Hessian_pattern(docp::DOCP{Trapeze})
+
+    # NB. need to provide full pattern for coloring, not just upper/lower part
+    Is = Vector{Int}(undef, 0)
+    Js = Vector{Int}(undef, 0)
+
+    # index alias for v
+    v_start = docp.dim_NLP_variables - docp.dim_NLP_v + 1
+    v_end = docp.dim_NLP_variables
+
+    # 0. objective
+    # 0.1 mayer cost (x0, xf, v) 
+    # -> grouped with term 3. for boundary conditions
+    # 0.2 lagrange case (lf)
+    # -> 2nd order term is zero
+   
+    # 1. main loop over steps
+    # 1.0 v / v term
+    add_nonzero_block!(Is, Js, v_start, v_end, v_start, v_end)
+
+    for i = 1:docp.dim_NLP_steps
+
+        # contiguous variables blocks will be used when possible
+        # x_i (l_i) u_i x_i+1 (l_i+1) u_i+1
+        var_block = docp.discretization._step_variables_block * 2
+        var_offset = (i-1)*docp.discretization._step_variables_block
+
+        # 1.1 state eq 0 = x_i+1 - (x_i + h_i/2 (f(t_i,x_i,u_i,v) + f(t_i+1,x_i+1,u_i+1,v)))
+        # depends on x_i, u_i, x_i+1, u_i+1, and v -> included in 1.2
+        # 1.2 lagrange part 0 = l_i+1 - (l_i + h_i/2 (l(t_i,x_i,u_i,v) + l(t_i+1,x_i+1,u_i+1,v)))
+        # depends on x_i, l_i, u_i, x_i+1, l_i+1, u_i+1, and v
+        # -> use single block for all step variables
+        add_nonzero_block!(Is, Js, var_offset+1, var_offset+var_block, var_offset+1, var_offset+var_block)
+        add_nonzero_block!(Is, Js, var_offset+1, var_offset+var_block, v_start, v_end; sym=true)
+
+        # 1.3 path constraint g(t_i, x_i, u_i, v)
+        # -> included in 1.2
+    end
+
+    # 2. final path constraints (xf, uf, v)
+    # -> included in last loop iteration
+
+    # 3. boundary constraints (x0, xf, v)
+    # -> (x0, v) terms included in first loop iteration
+    # -> (xf, v) terms included in last loop iteration
+    if docp.is_mayer || docp.dim_boundary_cons > 0
+        var_offset = docp.dim_NLP_steps*docp.discretization._step_variables_block
+        x0_start = 1
+        x0_end = docp.dim_OCP_x
+        xf_start = var_offset + 1
+        xf_end = var_offset + docp.dim_OCP_x
+        add_nonzero_block!(Is, Js, x0_start, x0_end, xf_start, xf_end; sym=true)
+    end
+    # 3.1 null initial condition for lagrangian cost state l0
+    # -> 2nd order term is zero
+
+    # build and return sparse matrix
+    nnzj = length(Is)
+    Vs = ones(Bool, nnzj)
+    return sparse(Is, Js, Vs, docp.dim_NLP_variables, docp.dim_NLP_variables)
+
+end
+
+
+"""
+$(TYPEDSIGNATURES)
+
+Add block of nonzeros elements to a sparsity pattern 
+Format: boolean matrix (M) or index vectors (Is, Js) 
+Includes a more compact method for single element case
+Option to add the symmetric block also (eg for Hessian)
+Note: independent from discretization scheme
+"""
+function add_nonzero_block!(M, i_start, i_end, j_start, j_end; sym=false)
+    M[i_start:i_end, j_start:j_end] .= true
+    sym && (M[j_start:j_end, i_start:i_end] .= true)
+    return
+end
+function add_nonzero_block!(M, i, j; sym=false)
+    M[i,j] = true
+    sym && (M[j,i] = true)
+    return
+end
+function add_nonzero_block!(Is, Js, i_start, i_end, j_start, j_end; sym=false)
+    for i=i_start:i_end
+        for j=j_start:j_end
+            push!(Is, i)
+            push!(Js, j)
+            sym && push!(Is, j)
+            sym && push!(Js, i)
+        end
+    end
+    return
+end
+function add_nonzero_block!(Is, Js, i, j; sym=false)
+    push!(Is, i)
+    push!(Js, j)
+    sym && push!(Is, j)
+    sym && push!(Js, i)
+    return
+end
diff --git a/src/docp.jl b/src/docp.jl
index 9edf0972..ffcaacec 100644
--- a/src/docp.jl
+++ b/src/docp.jl
@@ -87,7 +87,7 @@ struct DOCP{T <: Discretization, X <: ScalVect, U <: ScalVect, V <: ScalVect, G
     _type_v::V
 
     # constructor
-    function DOCP(ocp::OptimalControlModel; grid_size=__grid_size(), time_grid=__time_grid(), disc_method=__disc_method(), constant_control=false)
+    function DOCP(ocp::OptimalControlModel; grid_size=__grid_size(), time_grid=__time_grid(), disc_method=__disc_method(), control_type=__control_type())
 
         # time grid
         if time_grid == nothing
@@ -189,9 +189,9 @@ struct DOCP{T <: Discretization, X <: ScalVect, U <: ScalVect, V <: ScalVect, G
         elseif disc_method == :gauss_legendre_1
                 discretization, dim_NLP_variables, dim_NLP_constraints = CTDirect.Gauss_Legendre_1(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons)
         elseif disc_method == :gauss_legendre_2
-                discretization, dim_NLP_variables, dim_NLP_constraints = CTDirect.Gauss_Legendre_2(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons, constant_control)
+                discretization, dim_NLP_variables, dim_NLP_constraints = CTDirect.Gauss_Legendre_2(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons, control_type)
         elseif disc_method == :gauss_legendre_3
-                discretization, dim_NLP_variables, dim_NLP_constraints = CTDirect.Gauss_Legendre_3(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons, constant_control)                                 
+                discretization, dim_NLP_variables, dim_NLP_constraints = CTDirect.Gauss_Legendre_3(dim_NLP_steps, dim_NLP_x, dim_NLP_u, dim_NLP_v, dim_u_cons, dim_x_cons, dim_xu_cons, dim_boundary_cons, dim_v_cons, control_type)                                 
         else           
             error("Unknown discretization method: ", disc_method, "\nValid options are disc_method={:trapeze, :midpoint, :gauss_legendre_1, :gauss_legendre_2, :gauss_legendre_3}\n", typeof(disc_method))
         end
@@ -387,7 +387,7 @@ function DOCP_constraints!(c, xu, docp::DOCP)
         #setStepConstraints!(docp, (@view c[offset+1:offset+docp.dim_NLP_x+docp.discretization._step_pathcons_block]), xu, v, time_grid, i, work)
     end
 
-    # point constraints (NB. view on c block could be used with offset here)
+    # point constraints
     setPointConstraints!(docp, c, xu, v)
 
     # NB. the function *needs* to return c for AD...
@@ -515,6 +515,12 @@ function setPointBounds!(docp::DOCP, index::Int, lb, ub)
     return index
 end
 
+
+function DOCP_Jac_pattern(docp::DOCP)
+    error("DOCP_Jac_pattern not implemented for discretization ", typeof(docp.discretization))
+end
+
+
 """
 $(TYPEDSIGNATURES)
 
diff --git a/src/solve.jl b/src/solve.jl
index 79d75e37..bc022d9c 100644
--- a/src/solve.jl
+++ b/src/solve.jl
@@ -13,53 +13,149 @@ function available_methods()
     return algorithms
 end
 
+
+"""
+$(TYPEDSIGNATURES)
+
+Solve an OCP with a direct method
+
+# Arguments
+* ocp: optimal control problem as defined in `CTBase`
+* [description]: can specifiy for instance the NLP model and / or solver
+
+# Keyword arguments (optional)
+* `grid_size`: number of time steps for the discretized problem ([250])
+* `disc_method`: discretization method ([`:trapeze`], `:midpoint`, `gauss_legendre_2`)
+* `time_grid`: explicit time grid (can be non uniform)
+* `init`: info for the starting guess (values or existing solution)
+* `adnlp_backend`: backend for automatic differentiation in ADNLPModels ([`:optimized`], `:manual`, `:default`)
+* `control_type`: ([`:constant`], `:linear`) control piecewise parametrization for IRK methods 
+
+All further keywords are passed to the inner call of `solve_docp`
+"""
+function direct_solve(
+    ocp::OptimalControlModel,
+    description::Symbol...;
+    grid_size::Int = CTDirect.__grid_size(),
+    disc_method = __disc_method(),
+    time_grid = CTDirect.__time_grid(),
+    init = CTBase.__ocp_init(),
+    adnlp_backend = __adnlp_backend(),
+    control_type = __control_type(),
+    kwargs...,
+)
+    method = getFullDescription(description, available_methods())
+
+    # build discretized OCP, including initial guess
+    docp, nlp = direct_transcription(
+        ocp,
+        description;
+        init = init,
+        grid_size = grid_size,
+        time_grid = time_grid,
+        disc_method = disc_method,
+        control_type = control_type,
+        adnlp_backend = adnlp_backend,
+    )
+
+    # solve DOCP
+    if :ipopt ∈ method
+        solver_backend = CTDirect.IpoptBackend()
+    elseif :madnlp ∈ method
+        solver_backend = CTDirect.MadNLPBackend()
+    else
+        error("no known solver in method", method)
+    end
+    docp_solution = CTDirect.solve_docp(solver_backend, docp, nlp; kwargs...)
+
+    # build and return OCP solution
+    return OptimalControlSolution(docp, docp_solution)
+end
+
+
 """
 $(TYPEDSIGNATURES)
 
 Discretize an optimal control problem into a nonlinear optimization problem (ie direct transcription)
+
+# Arguments
+* ocp: optimal control problem as defined in `CTBase`
+* [description]: can specifiy for instance the NLP model and / or solver
+
+# Keyword arguments (optional)
+* `grid_size`: number of time steps for the discretized problem ([250])
+* `disc_method`: discretization method ([`:trapeze`], `:midpoint`, `gauss_legendre_2`)
+* `time_grid`: explicit time grid (can be non uniform)
+* `init`: info for the starting guess (values or existing solution)
+* `adnlp_backend`: backend for automatic differentiation in ADNLPModels ([`:optimized`], `:manual`, `:default`)
+* `control_type`: ([`:constant`], `:linear`) control piecewise parametrization for IRK methods
+* show_time: (:true, [:false]) show timing details from ADNLPModels
+
 """
 function direct_transcription(
     ocp::OptimalControlModel,
     description...;
-    init = CTBase.__ocp_init(),
     grid_size = __grid_size(),
-    time_grid = __time_grid(),
     disc_method = __disc_method(),
-    constant_control = false,
-    adnlp_backend = __adnlp_backend()
+    time_grid = __time_grid(),
+    init = CTBase.__ocp_init(),
+    adnlp_backend = __adnlp_backend(),
+    control_type = __control_type(),
+    show_time = false
 )
 
     # build DOCP
-    docp = DOCP(ocp; grid_size=grid_size, time_grid=time_grid, disc_method=disc_method, constant_control = constant_control)
+    docp = DOCP(ocp; grid_size=grid_size, time_grid=time_grid, disc_method=disc_method, control_type = control_type)
 
     # set bounds in DOCP
     variables_bounds!(docp)
     constraints_bounds!(docp)
 
-    # set initial guess
-    x0 = DOCP_initial_guess(docp,
-        OptimalControlInit(init, state_dim = ocp.state_dimension, control_dim = ocp.control_dimension, variable_dim = ocp.variable_dimension)
-    )
+    # build and set initial guess in DOCP
+    docp_init = OptimalControlInit(init, state_dim = ocp.state_dimension, control_dim = ocp.control_dimension, variable_dim = ocp.variable_dimension)
+    x0 = DOCP_initial_guess(docp, docp_init)
+
+    # redeclare objective and constraints functions
+    f = x -> DOCP_objective(x, docp)
+    c! = (c, x) -> DOCP_constraints!(c, x, docp)
 
     # call NLP problem constructor
-    if adnlp_backend == :no_hessian
+    if adnlp_backend == :manual
+        
+        # build sparsity pattern
+        J_backend = ADNLPModels.SparseADJacobian(docp.dim_NLP_variables, f, docp.dim_NLP_constraints, c!, DOCP_Jacobian_pattern(docp))
+        H_backend = ADNLPModels.SparseReverseADHessian(docp.dim_NLP_variables, f, docp.dim_NLP_constraints, c!, DOCP_Hessian_pattern(docp))
+        
+        # build NLP with given patterns
         nlp = ADNLPModel!(
-            x -> DOCP_objective(x, docp), x0, docp.var_l, docp.var_u,
-            (c, x) -> DOCP_constraints!(c, x, docp), docp.con_l, docp.con_u,
-            backend = __adnlp_backend()
-            )
-        set_adbackend!(nlp, hessian_backend = ADNLPModels.EmptyADbackend, hvprod_backend = ADNLPModels.EmptyADbackend) # directionalsecondderivative)
+        f, x0, docp.var_l, docp.var_u, c!, docp.con_l, docp.con_u,
+        gradient_backend = ADNLPModels.ReverseDiffADGradient,
+        jacobian_backend = J_backend,
+        hessian_backend = H_backend,
+        hprod_backend = ADNLPModels.EmptyADbackend,
+        jtprod_backend = ADNLPModels.EmptyADbackend,
+        jprod_backend = ADNLPModels.EmptyADbackend,
+        ghjvprod_backend = ADNLPModels.EmptyADbackend,
+        show_time = show_time,
+        #excluded_backend = [:jprod_backend, :jtprod_backend, :hprod_backend, :ghjvprod_backend]
+    )
     else
+        # build NLP
         nlp = ADNLPModel!(
-            x -> DOCP_objective(x, docp), x0, docp.var_l, docp.var_u,
-            (c, x) -> DOCP_constraints!(c, x, docp), docp.con_l, docp.con_u,
-            backend = adnlp_backend
+            f, x0, docp.var_l, docp.var_u, c!, docp.con_l, docp.con_u,
+            backend = adnlp_backend, 
+            hprod_backend = ADNLPModels.EmptyADbackend,
+            jtprod_backend = ADNLPModels.EmptyADbackend,
+            jprod_backend = ADNLPModels.EmptyADbackend,
+            ghjvprod_backend = ADNLPModels.EmptyADbackend,       
+            show_time = show_time,
             )
     end
 
     return docp, nlp
 end
 
+
 """
 $(TYPEDSIGNATURES)
 
@@ -78,49 +174,6 @@ function set_initial_guess(docp::DOCP, nlp, init)
     )
 end
 
-"""
-$(TYPEDSIGNATURES)
-
-Solve an OCP with a direct method
-"""
-function direct_solve(
-    ocp::OptimalControlModel,
-    description::Symbol...;
-    init = CTBase.__ocp_init(),
-    grid_size::Int = CTDirect.__grid_size(),
-    time_grid = CTDirect.__time_grid(),
-    disc_method = __disc_method(),
-    constant_control = false,
-    adnlp_backend = __adnlp_backend(),
-    kwargs...,
-)
-    method = getFullDescription(description, available_methods())
-
-    # build discretized OCP, including initial guess
-    docp, nlp = direct_transcription(
-        ocp,
-        description;
-        init = init,
-        grid_size = grid_size,
-        time_grid = time_grid,
-        disc_method = disc_method,
-        constant_control = constant_control,
-        adnlp_backend = adnlp_backend
-    )
-
-    # solve DOCP
-    if :ipopt ∈ method
-        solver_backend = CTDirect.IpoptBackend()
-    elseif :madnlp ∈ method
-        solver_backend = CTDirect.MadNLPBackend()
-    else
-        error("no known solver in method", method)
-    end
-    docp_solution = CTDirect.solve_docp(solver_backend, docp, nlp; kwargs...)
-
-    # build and return OCP solution
-    return OptimalControlSolution(docp, docp_solution)
-end
 
 # placeholders (see CTSolveExt*** extensions)
 abstract type AbstractSolverBackend end
diff --git a/test/ab_jump.jl b/test/ab_jump.jl
index ebec29ea..421aa1d4 100644
--- a/test/ab_jump.jl
+++ b/test/ab_jump.jl
@@ -27,6 +27,7 @@ function algal_bacterial_jump(;grid_size=1000, disc_method=:trapeze, print_level
     set_optimizer_attribute(sys, "max_iter", 1500)
     set_optimizer_attribute(sys, "mu_strategy", "adaptive")
     set_optimizer_attribute(sys, "sb", "yes")
+    set_optimizer_attribute(sys, "print_user_options", "yes")
 
     # Discretization parameters
     N = grid_size
diff --git a/test/benchmark.jl b/test/benchmark.jl
index e59c7b4c..11f6bf55 100644
--- a/test/benchmark.jl
+++ b/test/benchmark.jl
@@ -24,7 +24,7 @@ end
 
 
 
-function bench_list(problem_list; verbose=2, nlp_solver, linear_solver, kwargs...)
+function bench_list(problem_list; verbose=1, nlp_solver, linear_solver, kwargs...)
 
     #######################################################
     # solve examples with timer and objective check
@@ -66,10 +66,11 @@ function bench(;grid_size_list = [250, 500, 1000, 2500, 5000], verbose = 1, nlp_
     verbose > 1 && @printf("Blas config: %s\n", LinearAlgebra.BLAS.lbt_get_config())
 
     # load problems for benchmark
+    # Note that problems may vary significantly in convergence times...  
     if names_list == :default
-        names_list = ["beam", "double_integrator_mintf", "double_integrator_minenergy", "double_integrator_freet0tf", "fuller", "goddard", "goddard_all", "jackson", "robbins", "simple_integrator", "vanderpol"]
+        names_list = ["beam", "double_integrator_mintf", "double_integrator_minenergy", "double_integrator_freet0tf", "fuller", "goddard", "goddard_all", "jackson", "simple_integrator", "vanderpol"]
     elseif names_list == :quick
-        names_list = ["beam", "double_integrator_mintf", "fuller", "jackson", "robbins", "simple_integrator", "vanderpol"]
+        names_list = ["beam", "double_integrator_mintf", "fuller", "jackson", "simple_integrator", "vanderpol"]
     elseif names_list == :all 
         names_list = ["algal_bacterial", "beam", "bioreactor_1day", "bioreactor_Ndays", "bolza_freetf", "double_integrator_mintf", "double_integrator_minenergy", "double_integrator_freet0tf", "fuller", "goddard", "goddard_all", "insurance", "jackson", "robbins", "simple_integrator", "swimmer", "vanderpol"]
     elseif names_list == :hard
@@ -87,9 +88,9 @@ function bench(;grid_size_list = [250, 500, 1000, 2500, 5000], verbose = 1, nlp_
     for grid_size in grid_size_list
         t = bench_list(problem_list; grid_size=grid_size, verbose=verbose, nlp_solver=nlp_solver, linear_solver=linear_solver, kwargs...)
         append!(t_list, t)
-        @printf("Grid size %d: time (s) = %6.1f\n", grid_size, t)
+        @printf("Grid size %6d: time (s) = %6.1f\n", grid_size, t)
     end
-
+    #return t_list
 end
 
 
diff --git a/test/docs/AD_backend.md b/test/docs/AD_backend.md
index 77a27644..78c1f743 100644
--- a/test/docs/AD_backend.md
+++ b/test/docs/AD_backend.md
@@ -1,16 +1,103 @@
-Benchmark CTDirect for different AD backends (keyword adnlp_backend = ...)
-- :default : ForwardDiff
-- :optimized (default for CTDirect) : ForwardDiff for Jacobian / ReverseDiff for Hessian
-- :enzyme : Enzyme
-- :zygote : Zygote
+# Benchmark for different AD backends
+The backend for ADNLPModels can be set in transcription / solve calls with the option `adnlp_backend=`. Possible values include the predefined(*) backends for ADNLPModels:
+- `:optimized`* Default for CTDirect. Forward mode for Jacobian, reverse for Gradient and forward over reverse for Hessian.
+- `:default`* Forward mode only. Significantly slower, but more rugged.
+- `:manual` Explicitely give to ADNLPModels the sparse pattern for Jacobian and Hessian. Uses the same forward / reverse settings as the `:optimized` predefined backend.  
+- `:enzyme`* Enzyme (currently not working).
 
-Using CTDirect benchmark function bench()
-Problem list: ["beam", "double_integrator_mintf", "double_integrator_minenergy", "double_integrator_freet0tf", "fuller", "goddard", "goddard_all", "jackson", "robbins", "simple_integrator", "vanderpol"]
+## Tests:
+```
+julia> include("test/benchmark.jl")
+test_unit (generic function with 1 method)
+
+julia> bench(grid_size_list=[250,500,1000,2500,5000,7500,10000], adnlp_backend=:manual)
+Problem list: ["beam", "double_integrator_mintf", "double_integrator_minenergy", "double_integrator_freet0tf", "fuller", "goddard", "goddard_all", "jackson", "simple_integrator", "vanderpol"]
+```
 
 Takeaways:
-- optimized backend (with ReverseDiff for Hessian) is much better than full ForwardDiff.
+- the `:optimized` backend (with forward over reverse mode for Hessian) is much faster than full forward mode, but does not scale greatly. This is likely due to the increasing cost of computing the Hessian sparsity with SparseConnectivityTracer.jl in terms of allocations and time. Note that the `:default` backend that uses forward mode only may still be useful when having AD errors.
+- manual sparse pattern seems to give better performance for larger problems. See also the comparison with Jump that seems to use a different, less sparse but faster method for the Hessian. The sparsity pattern detection in JuMP relies on the expression tree of the objective and constraints built from its DSL.
+
+![benchmark](AD_backend.png)
+
+Standard benchmark for Trapeze:
+| Trapeze | default | optimized | manual   |
+|---------|---------|-----------|----------|
+| 250     | 51.9    | 0.8       | 1.4      |
+| 500     | 219.0   | 2.0       | 3.1      |
+| 1000    | 858.8   | 5.2       | 5.4      |
+| 2500    | 6932.1  | 20.9      | 17.1     |
+| 5000    |         | 73.0      | 34.1     |
+| 7500    |         | 200.6     | 53.1     |
+| 10000   |         | 415.7     | 70.4     |
+
+* (older version) build sparse matrices from dense boolean matrices
+** build sparse matrices from (i,j,v) vectors
+
+Standard benchmark for Midpoint:
+| Midpoint| optimized | manual |
+|---------|-----------|--------|
+| 250     | 1.4       | 2.1    |
+| 500     | 3.4       | 4.2    |
+| 1000    | 10.0      | 10.2   |
+| 2500    | 45.0      | 30.0   |
+| 5000    | 150.1     | 79.0   |
+| 7500    | 322.6     | 130.7  |
+
+Standard benchmark for Gauss Legendre 2:
+| GL2     | optimized | manual |
+|---------|-----------|--------|
+| 250     | 3.5       | 4.3    |
+| 500     | 9.6       | 11.5   |
+| 1000    | 125.5     | 22.3   |
+| 2500    | 135.0     | 68.3   |
+| 5000    | 527.7     | 156.8  |
+
+Sparsity details: goddard_all Trapeze (1000 and 10000 steps)
+| transcription | optimized | manual     | optimized | manual  |
+|---------------|-----------|------------|-----------|---------|
+| NLP vars      | 4005      | 4005       | 40005     | 40005   |
+| NLP cons      | 6007      | 6007       | 60007     | 60007   |
+| Hess nnz      | 11011     | 30024      | 110011    | 300024  |
+| H sparsity    | 99.86%    | 99.63%     | 99.99%    | 99.96%  |
+| Jac nnz       | 28011     | 42043      | 280011    | 420043  |
+| J sparsity    | 99.88%    | 99.83%     | 99.99%    | 99.98%  |
+| allocs        | 1.2GB     | 92MB       | 71.6GB    | 0.88 GB |
+| time          | 750ms     | 95ms       | 64.7s***  | 2.5s    |
+
+*** hessian accounts for 59 out of total 65s
+```
+julia> direct_transcription(goddard_all().ocp, grid_size=10000, show_time=true);
+gradient backend ADNLPModels.ReverseDiffADGradient: 0.000137972 seconds;
+hprod    backend ADNLPModels.ReverseDiffADHvprod: 0.314931491 seconds;
+jprod    backend ADNLPModels.ForwardDiffADJprod: 2.2412e-5 seconds;
+jtprod   backend ADNLPModels.ReverseDiffADJtprod: 0.612174104 seconds;
+jacobian backend ADNLPModels.SparseADJacobian: 0.425535048 seconds;
+hessian  backend ADNLPModels.SparseReverseADHessian: 58.450146911 seconds;
+ghjvprod backend ADNLPModels.ForwardDiffADGHjvprod: 4.339e-6 seconds.
+```
+
+| solve         | optimized | manual  | optimized | manual  |
+|---------------|-----------|---------|-----------|---------|
+| iterations    | 42        | 28      | 51        | 29      |
+| allocs        | 2.0GB     | 1.2GB   | 87.5GB    | 13.2GB  |
+| time          | 2.5s      | 2.6s    | 151.0s*** | 31.6s   |
+
+*** building the hessian is one third of the total solve time...
+
+
+## Remarks:
+- it is better to build the sparse matrices from the index vectors format rather than a dense boolean matrix (10-20% faster). For larger problems it may not be possible to even allocate the boolean matrix (eg. algal bacterial with GL2 at 5000 steps).
+- disabling the unused backends for the various matrix products (jprod_backend, jtprod_backend, hprod_backend, ghjvprod_backend) gives a slight increase in performance.
+
+## Todo:
+- use automatic differentiation to get the sparsity patterns for first / second derivatives of the OCP functions, and build the Jacobian / Hessian patterns from these instead of assuming full nonzero blocks.
+- some gain may be achieved by preallocating the index vectors
+
+## Errors for Enzyme:
 - enzyme gives correct nonzero counts for Jacobian and Hessian, but fails with
-```ERROR: Constant memory is stored (or returned) to a differentiable variable.
+```
+ERROR: Constant memory is stored (or returned) to a differentiable variable.
 As a result, Enzyme cannot provably ensure correctness and throws this error.
 This might be due to the use of a constant variable as temporary storage for active memory (https://enzyme.mit.edu/julia/stable/faq/#Runtime-Activity).
 If Enzyme should be able to prove this use non-differentable, open an issue!
@@ -18,17 +105,4 @@ To work around this issue, either:
  a) rewrite this variable to not be conditionally active (fastest, but requires a code change), or
  b) set the Enzyme mode to turn on runtime activity (e.g. autodiff(set_runtime_activity(Reverse), ...) ). This will maintain correctness, but may slightly reduce performance.```
  Error apparently occurs when calling the boundary conditions.
-- zygote gives incorrect (huge) nonzero counts then also fails with an error message. 
-
-
-| Trapeze | default | optimized |
-|---------|---------|-----------|
-| 250     | 43.3    | 1.0       |
-| 500     | 176.2   | 2.4       |
-| 1000    | 926.0   | 7.1       |
-| 2500    |         | 31.8      |
-| 5000    |         |           |
-
-Midpoint
-
-GaussLegendre2
+ ```
diff --git a/test/docs/AD_backend.png b/test/docs/AD_backend.png
new file mode 100644
index 00000000..a73e2d70
Binary files /dev/null and b/test/docs/AD_backend.png differ
diff --git a/test/docs/jump_ctdirect.md b/test/docs/jump_ctdirect.md
index 53603e4f..e59058bb 100644
--- a/test/docs/jump_ctdirect.md
+++ b/test/docs/jump_ctdirect.md
@@ -1,88 +1,70 @@
 # Jump / CTDirect comparison - algal bacterial problem
 
-Note that the problem is redefined for each method: jump, ctdirect and ctdirect new model.
-Also, the Gauss Legendre 2 implementations for Jump and CTDirect here use a piecewise constant control (default for CTDirect would have been piecewise linear).
+Note that the problem is redefined for each method, jump and ctdirect.
+Also, the Gauss Legendre 2 implementations use a piecewise constant control.
 
 ## Takeaways
-- CTDirect still allocates at least x10 more memory, worsening for higher problem sizes
-For the biggest allocations, a significant time is passed during the AD phase, before Ipopt. The runs marked with * spend half the time before optimization, likely some swap issue.
-We note that Jump memory appears linear wrt steps for GL2, but a bit superlinear for Trapeze. CTDirect memory always increases superlinearly wrt steps.
-- Hessian seems to be handled differently by Jump, see the higher nonzero values.
-Maybe a less sparse but faster and less memory intensive method is used ? 
-- convergence: objective and trajectory are similar, iterations differ, maybe due to the different hessian handling. Total computation times are similar for Trapeze and x2 to x5 slower for CTDirect for GL2, probably due to the memory effect.
-- for GL2, Jump and CTDirect have slightly different nonzero counts for the Jacobian
-- in terms of control structures, GL2 solutions are clean, Jump Trapeze solutions shows a bit of noise, while CTDirect Trapeze solutions are very noisy.
+- Hessian is handled differently by Jump, with a less sparse but faster method.
+- CTDirect with manual sparsity patterns still allocates 5 to 10 times more memory than Jump, but scales better than the automatic sparsity detection. Manual mode becomes faster than automatic mode for GL2 above 2000 steps, while being slower for Trapeze even at 5000 steps.
+- convergence: objective and trajectory are similar, iterations differ, maybe due to the different hessian handling. Total computation times for Trapeze are similar for Jump and CTDirect (auto). For GL2, Jump is faster than both CTDirect versions.
+- in terms of control structures, GL2 solutions are clean, Jump Trapeze solutions shows a bit of noise, while CTDirect Trapeze solutions are very noisy. How Jump manages to find a cleaner solution with Trapeze is unclear.
 
 ## Todo
-- can we have linear memory wrt steps for Jump / Trapeze ?
-- disable Hessian (in AD model then use ipopt limited memory option ?) and compare memory allocations and convergence. Find more details on the Hessian in Jump.
+- investigate how jump finds a cleaner solution for trapeze discretization (print settings ?)
 
 ## Results: Jump vs CTDirect
 See `test/jump_comparison.jl`
-Ipopt details: `Ipopt version 3.14.17, running with linear solver MUMPS 5.7.3`
+Ipopt details: `This is Ipopt version 3.14.14, running with linear solver MUMPS 5.6.2.`
 Settings: tol=1e-8, mu_strategy=adaptive
-
-```
-Jump trapeze 1000:  17.029 s (7920527 allocations: 351.87 MiB)
-Jump trapeze 2000:  56.928 s (23055273 allocations: 891.64 MiB)
-Jump trapeze 5000:  125.798 s (55226843 allocations: 2.10 GiB)
-Jump gauss_legendre_2 1000:  15.398 s (10988856 allocations: 726.32 MiB)
-Jump gauss_legendre_2 2000:  27.401 s (21345532 allocations: 1.40 GiB)
-Jump gauss_legendre_2 5000:  76.593 s (56269715 allocations: 3.57 GiB)
 ```
+Jump trapeze 1000:  15.633 s (7920529 allocations: 351.87 MiB)
+Jump trapeze 2000:  57.472 s (23055275 allocations: 891.64 MiB)
+Jump trapeze 5000:  124.108 s (55226845 allocations: 2.10 GiB)
+Jump gauss_legendre_2 1000:  15.250 s (10998729 allocations: 726.48 MiB)
+Jump gauss_legendre_2 2000:  26.686 s (21371405 allocations: 1.40 GiB)
+Jump gauss_legendre_2 5000:  76.455 s (56343588 allocations: 3.58 GiB)
 
-```
-CTDirect trapeze 1000:  20.110 s (46501059 allocations: 4.54 GiB)
-CTDirect trapeze 2000:  41.097 s (89302125 allocations: 12.26 GiB)
-CTDirect trapeze 5000:  133.268 s (267989400 allocations: 49.33 GiB)
-```
-GL2 piecewise constant control
-```
-CTDirect gauss_legendre_2 1000:  33.181 s (37843213 allocations: 14.79 GiB)
-CTDirect gauss_legendre_2 2000:  82.605 s (82766476 allocations: 43.19 GiB)
-CTDirect gauss_legendre_2 5000:  356.338 s (221161426 allocations: 312.28 GiB)
-```
-GL2 piecewise linear control
-```
-CTDirect gauss_legendre_2 1000:  37.220 s (39259673 allocations: 15.12 GiB)
-CTDirect gauss_legendre_2 2000:  112.687 s (104950745 allocations: 45.37 GiB)
-CTDirect gauss_legendre_2 5000:  363.224 s (211848763 allocations: 313.02 GiB)
+CTDirect (optimized) trapeze 1000:  17.941 s (45041839 allocations: 4.47 GiB)
+CTDirect (optimized) trapeze 2000:  35.811 s (86384903 allocations: 12.12 GiB)
+CTDirect (optimized) trapeze 5000:  124.451 s (260698176 allocations: 48.99 GiB)
+CTDirect (optimized) gauss_legendre_2 1000:  25.087 s (32172053 allocations: 14.36 GiB)
+CTDirect (optimized) gauss_legendre_2 2000:  76.272 s (77043394 allocations: 42.43 GiB)
+CTDirect (optimized) gauss_legendre_2 5000:  281.000 s (138053972 allocations: 303.56 GiB)
+
+CTDirect (manual) trapeze 1000:  47.442 s (65149511 allocations: 5.33 GiB)
+CTDirect (manual) trapeze 2000:  105.765 s (137360269 allocations: 11.26 GiB)
+CTDirect (manual) trapeze 5000:  324.819 s (410732633 allocations: 33.73 GiB)
+CTDirect (manual) gauss_legendre_2 1000:  38.147 s (51843118 allocations: 4.00 GiB)
+CTDirect (manual) gauss_legendre_2 2000:  87.863 s (118950784 allocations: 9.16 GiB)
+CTDirect (manual) gauss_legendre_2 5000:  187.292 s (241939016 allocations: 18.85 GiB)
 ```
 
 ## Details: Trapeze (1000 and 5000 steps)
 
-|                 | Jump   | CT     | New    | Jump     | CT       | New      |
+|                 | Jump   | CT     | Manual | Jump     | CT       | Manual   |
 |-----------------|--------|--------|--------|----------|----------|----------|
-|nnz jacobian     | 42006  | 42006  |        | 210006   | 210006   |          |
-|nnz hessian      | 74000  | 12012  |        | 370000   | 60012    |          |
-|variables        | 8008   | 8008   |        | 40008    | 40008    |     |
-|lowerbound       | 6006   | 6006   |        | 30006    | 30006    |     |
-|lower/upper      | 2002   | 2002   |        | 10002    | 10002    |     |
-|equality         | 6006   | 6006   |        | 30006    | 30006    |     |
-|iterations       | 334    | 365    |        | 517      | 420      |       |
-|objective        | 5.4522 | 5.4522 |        | 5.4522   | 5.4522   |    |
-|structure        | ok     | noisy  |        | ok       | noisy    |          |
-|allocations      | 352MB  | 4.5GB  |        | 2.1GB    | 49GB     |          |
-|time             | 17     | 20     |        | 126      | 136      |          |
-
+|nnz jacobian     | 42006  | 42006  | 96076  | 210006   | 210006   | 480072   |
+|nnz hessian      | 74000  | 12012  | 100072 | 370000   | 60012    | 500072   |
+|variables        | 8008   | 8008   | 8008   | 40008    | 40008    | 40008    |
+|lowerbound       | 6006   | 6006   | 6006   | 30006    | 30006    | 30006    |
+|lower/upper      | 2002   | 2002   | 2002   | 10002    | 10002    | 10002    |
+|equality         | 6006   | 6006   | 6006   | 30006    | 30006    | 30006    |
+|iterations       | 334    | 365    | 333    | 517      | 420      | 419      |
+|objective        | 5.4522 | 5.4522 | 5.4522 | 5.4522   | 5.4522   | 5.4522   |
+|structure        | ok     | noisy  | noisy  | ok       | noisy    | noisy    |
 
 ## Details: Gauss Legendre 2 (1000 and 5000 steps)
 
-|                 | Jump   | CT     | New    | Jump     | CT       | New      |
+|                 | Jump   | CT     | Manual | Jump     | CT       | Manual   |
 |-----------------|--------|--------|--------|----------|----------|----------|
-|nnz jacobian     | 118006 | 124000 |        | 590006   | 620000   |          |
-|nnz hessian      | 322000 | 63000  |        | 1610000  | 315000   |          |
-|variables        | 20006  | 20008  |   | 100006   | 100008   |    |
-|lowerbound       | 6006   | 6006   |    | 3006     | 30006    |     |
-|lower/upper      | 2000   | 2002   |    | 10000    | 10002    |     |
-|equality         | 18006  | 18006  |   | 90006    | 90006    |     |
-|iterations       | 117    | 96     |        | 146      | 119      |          |
-|objective        | 5.4522 | 5.4522 |  | 5.4522   | 5.4522   |    |
-|structure        | clean  | clean  |        | clean    | clean    |          |
-|allocations      | 726MB  | 14.8GB |        | 3.6GB    | 312GB    |          |
-|time             | 15     | 33     |        | 77       | 356*     |          |
+|nnz jacobian     | 118006 | 118006 | 384072 | 590006   | 590006   | 1920072  |
+|nnz hessian      | 322000 | 63000  | 210072 | 1610000  | 315000   | 1050072  |
+|variables        | 20006  | 20008  | 20008  | 100006   | 100008   | 100008   |
+|lowerbound       | 6006   | 6006   | 6006   | 3006     | 30006    | 30006    |
+|lower/upper      | 2000   | 2002   | 2002   | 10000    | 10002    | 10002    |
+|equality         | 18006  | 18006  | 18006  | 90006    | 90006    | 90006    |
+|iterations       | 117    | 95     | 93     | 146      | 78       | 86       |
+|objective        | 5.4522 | 5.4522 | 5.4522 | 5.4522   | 5.4522   | 5.4522   |
+|structure        | clean  | clean  | clean  | clean    | clean    | clean    |
 
 * half the time is before optimization, swap effect due to huge allocations ?
-
-
-
diff --git a/test/jump_comparison.jl b/test/jump_comparison.jl
index fa1c792f..8d9d4d8f 100644
--- a/test/jump_comparison.jl
+++ b/test/jump_comparison.jl
@@ -6,12 +6,12 @@ using MKL
 using BenchmarkTools
 using Printf
 
-
 jump = true
-ctdirect = false
+ctdirect = true
+adnlp_backend_list = [:manual, :optimized]
 #disc_method_list = [:gauss_legendre_2]
-disc_method_list = [:trapeze]
-#disc_method_list = [:trapeze, :gauss_legendre_2]
+#disc_method_list = [:trapeze]
+disc_method_list = [:trapeze, :gauss_legendre_2]
 grid_size_list = [1000, 2000, 5000] 
 
 # Jump
@@ -28,10 +28,12 @@ end
 # CTDirect
 include("problems/algal_bacterial.jl")
 if ctdirect
-for disc_method in disc_method_list
-    for grid_size in grid_size_list
-        @printf("CTDirect %s %d:", disc_method, grid_size)
-        @btime direct_solve(algal_bacterial().ocp, grid_size=$grid_size, disc_method=$disc_method, print_level=0, constant_control=true)
+    for backend in adnlp_backend_list
+        for disc_method in disc_method_list
+            for grid_size in grid_size_list
+                @printf("CTDirect (%s) %s %d:", backend, disc_method, grid_size)
+                @btime direct_solve(algal_bacterial().ocp, grid_size=$grid_size, disc_method=$disc_method, print_level=0, adnlp_backend=$backend)
+            end
+        end
     end
-end
 end
\ No newline at end of file
diff --git a/test/problems/pattern.jl b/test/problems/pattern.jl
new file mode 100644
index 00000000..bf8ed990
--- /dev/null
+++ b/test/problems/pattern.jl
@@ -0,0 +1,15 @@
+# Duumy problem to visualize sparsity patterns
+
+function pattern()
+    @def ocp begin
+        t ∈ [0, 1], time
+        x ∈ R, state
+        u ∈ R, control
+        v ∈ R, variable
+        x(0) + x(1) + v == 0
+        ẋ(t) == x(t)^2 + u(t)^2 + v^2
+        ∫(u(t)^2 + x(t)^2 + v^2) → min
+    end
+
+    return ((ocp = ocp, obj = nothing, name = "pattern", init = nothing))
+end
\ No newline at end of file
diff --git a/test/suite/test_nlp.jl b/test/suite/test_nlp.jl
index a7a43936..26ad253e 100644
--- a/test/suite/test_nlp.jl
+++ b/test/suite/test_nlp.jl
@@ -17,8 +17,12 @@ end
     @test sol.objective ≈ prob.obj rtol = 1e-2
     sol = direct_solve(prob.ocp, display = false, adnlp_backend = :default)
     @test sol.objective ≈ prob.obj rtol = 1e-2
-    #sol = direct_solve(prob.ocp, display = false, adnlp_backend = :enzyme)
-    #@test sol.objective ≈ prob.obj rtol = 1e-2
+    sol = direct_solve(prob.ocp, display = false, adnlp_backend = :manual)
+    @test sol.objective ≈ prob.obj rtol = 1e-2
+    sol = direct_solve(prob.ocp, display = false, disc_method=:midpoint, adnlp_backend = :manual)
+    @test sol.objective ≈ prob.obj rtol = 1e-2
+    sol = direct_solve(prob.ocp, display = false, disc_method=:gauss_legendre_2, adnlp_backend = :manual)
+    @test sol.objective ≈ prob.obj rtol = 1e-2
 end
 
 # DOCP solving