diff --git a/Project.toml b/Project.toml
index 9af68ca..6eced7f 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,9 +14,6 @@ ExaModels = "1037b233-b668-4ce9-9b63-f9f681f55dd2"
 COPSBenchmarkExaModels = "ExaModels"
 COPSBenchmarkJuMP = "JuMP"
 
-[sources]
-ExaModels = {url = "https://github.com/exanauts/ExaModels.jl", rev = "main"}
-
 [compat]
 JuMP = "^1.19"
 julia = "1.11"
diff --git a/ext/COPSBenchmarkExaModels/catmix.jl b/ext/COPSBenchmarkExaModels/catmix.jl
index cd4ba87..34d419a 100644
--- a/ext/COPSBenchmarkExaModels/catmix.jl
+++ b/ext/COPSBenchmarkExaModels/catmix.jl
@@ -7,51 +7,53 @@
     ne = 2
     nc = 3
 
-    tf = 1
-    h = tf / nh   # Final time
+    h = T(1) / T(nh)   # Final time / nh (was Int/Int → Float64)
 
-    rho = [
+    rho = T[
         0.11270166537926,
         0.50000000000000,
         0.88729833462074,
     ]
-    bc = [1.0, 0.0]  # Boundary conditions for x
-    alpha = 0.0      # Smoothing parameter
+    bc = T[1, 0]   # Boundary conditions for x
+    alpha = T(0)   # Smoothing parameter
+    neg_one = -T(1)
+    ten_T = T(10)
+    one_T = T(1)
     rho_index = [(i, rho[i]) for i in 1:nc]
 
     c = ExaModels.ExaCore(T; backend = backend, concrete = Val(true))
-    ExaModels.@add_var(c, u, nh, nc; lvar = zeros(nh, nc), uvar = ones(nh, nc), start = zeros(nh, nc))
-    ExaModels.@add_var(c, v, nh, ne; start = [mod(j, ne) for i in 1:nh, j in 1:ne])
-    ExaModels.@add_var(c, w, nh, nc, ne; start = zeros(nh, nc, ne))
-    ExaModels.@add_var(c, pp, nh, nc, ne; start = [mod(k, ne) for i in 1:nh, j in 1:nc, k in 1:ne])
-    ExaModels.@add_var(c, Dpp, nh, nc, ne; start = zeros(nh, nc, ne))
-    ExaModels.@add_var(c, ppf, ne; start = [mod(i,ne) for i in 1:ne])
-
-    ExaModels.@add_obj(c, -1.0 + ppf[1] + ppf[2])
+    ExaModels.@add_var(c, u, nh, nc; lvar = zeros(T, nh, nc), uvar = ones(T, nh, nc), start = zeros(T, nh, nc))
+    ExaModels.@add_var(c, v, nh, ne; start = T[mod(j, ne) for i in 1:nh, j in 1:ne])
+    ExaModels.@add_var(c, w, nh, nc, ne; start = zeros(T, nh, nc, ne))
+    ExaModels.@add_var(c, pp, nh, nc, ne; start = T[mod(k, ne) for i in 1:nh, j in 1:nc, k in 1:ne])
+    ExaModels.@add_var(c, Dpp, nh, nc, ne; start = zeros(T, nh, nc, ne))
+    ExaModels.@add_var(c, ppf, ne; start = T[mod(i,ne) for i in 1:ne])
+
+    ExaModels.@add_obj(c, neg_one + ppf[1] + ppf[2])
     ExaModels.@add_obj(c, alpha/h*(u[i+1, j] - u[i, j])^2 for i in 1:nh-1, j in 1:nc)
 
     ExaModels.@add_con(
         c,
         c1,
-        pp[i, k, s] - v[i, s] - h*sum(w[i, j, s]*(rho^j/factorial(j)) for j in 1:nc) for i=1:nh,  (k, rho) in rho_index, s=1:ne
+        pp[i, k, s] - v[i, s] - h*sum(w[i, j, s]*(rho^j/T(factorial(j))) for j in 1:nc) for i=1:nh,  (k, rho) in rho_index, s=1:ne
     )
 
     ExaModels.@add_con(
         c,
         c2,
-        Dpp[i, k, s] - sum(w[i, j, s]*(rho^(j-1)/factorial(j-1)) for j in 1:nc) for i=1:nh, (k, rho) in rho_index, s=1:ne
+        Dpp[i, k, s] - sum(w[i, j, s]*(rho^(j-1)/T(factorial(j-1))) for j in 1:nc) for i=1:nh, (k, rho) in rho_index, s=1:ne
     )
 
     ExaModels.@add_con(
         c,
         c3,
-        ppf[s] - v[nh, s] - h * sum(w[nh, j, s] / factorial(j) for j in 1:nc) for s in 1:ne
+        ppf[s] - v[nh, s] - h * sum(w[nh, j, s] / T(factorial(j)) for j in 1:nc) for s in 1:ne
     )
 
     ExaModels.@add_con(
         c,
         c4,
-        v[i, s] + sum(w[i, j, s] * h / factorial(j) for j in 1:nc) - v[i+1, s] for i in 1:nh-1, s in 1:ne
+        v[i, s] + sum(w[i, j, s] * h / T(factorial(j)) for j in 1:nc) - v[i+1, s] for i in 1:nh-1, s in 1:ne
     )
 
 
@@ -59,13 +61,13 @@
     ExaModels.@add_con(
         c,
         c5,
-        Dpp[i,j,1] - u[i,j] * (10.0*pp[i,j,2] - pp[i,j,1]) for i=1:nh, j=1:nc
+        Dpp[i,j,1] - u[i,j] * (ten_T*pp[i,j,2] - pp[i,j,1]) for i=1:nh, j=1:nc
     )
 
     ExaModels.@add_con(
         c,
         c6,
-        Dpp[i,j,2] - u[i,j] * (pp[i,j,1] - 10.0*pp[i,j,2]) + (1 - u[i,j])*pp[i,j,2] for i=1:nh, j=1:nc
+        Dpp[i,j,2] - u[i,j] * (pp[i,j,1] - ten_T*pp[i,j,2]) + (one_T - u[i,j])*pp[i,j,2] for i=1:nh, j=1:nc
     )
 
 
@@ -77,5 +79,3 @@
 
     return ExaModels.ExaModel(c; kwargs...)
 end
-
-
diff --git a/ext/COPSBenchmarkExaModels/gasoil.jl b/ext/COPSBenchmarkExaModels/gasoil.jl
index 03c2284..8b49476 100644
--- a/ext/COPSBenchmarkExaModels/gasoil.jl
+++ b/ext/COPSBenchmarkExaModels/gasoil.jl
@@ -12,18 +12,19 @@
     nm = 21       # number of measurements
 
     # roots of k-th degree Legendre polynomial
-    rho = [0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703]
+    rho = T[0.06943184420297, 0.33000947820757, 0.66999052179243, 0.93056815579703]
     # ODE initial conditions
     bc = [1, 1, 2, 0]
     # times at which observations made
-    tau = [0.0, 0.025, 0.05, 0.075, 0.10, 0.125, 0.150, 0.175, 0.20, 0.225, 0.250, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.65, 0.75, 0.85, 0.95]
+    tau = T[0.0, 0.025, 0.05, 0.075, 0.10, 0.125, 0.150, 0.175, 0.20, 0.225, 0.250, 0.30, 0.35, 0.40, 0.45, 0.50, 0.55, 0.65, 0.75, 0.85, 0.95]
     # ODEs defined in [0,tf]
     tf = tau[nm]
     # uniform interval length
-    h = tf / nh
-    t = [(i-1)*h for i in 1:nh+1]
+    h = tf / T(nh)
+    t = T[T(i-1)*h for i in 1:nh+1]
+    zero_T = T(0)
 
-    itau = Int[min(nh, floor(tau[i]/h)+1) for i in 1:nm]
+    itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm]
 
     # Concentrations
     z = reshape(T[
@@ -50,10 +51,10 @@
         0.0690, 0.0100,
     ], ne, nm)'
 
-    v0 = zeros(nh, ne)
+    v0 = zeros(T, nh, ne)
     # Starting-value
     for i in 1:itau[1], s in 1:ne
-        v0[i, s] = bc[s]
+        v0[i, s] = T(bc[s])
     end
     for j in 2:nm, i =itau[j-1]+1:itau[j], s in 1:ne
         v0[i, s] = z[j, s]
@@ -65,31 +66,31 @@
     core = ExaModels.ExaCore(T; backend= backend, concrete = Val(true))
 
     # ODE parameters
-    ExaModels.@add_var(core, theta, 1:np; lvar = 0.0, start=0.0)
+    ExaModels.@add_var(core, theta, 1:np; lvar = zero_T, start=zero_T)
     # The collocation approximation u is defined by the parameters v and w.
     # uc and Duc are, respectively, u and u' evaluated at the collocation points.
     ExaModels.@add_var(core, v, 1:nh, 1:ne; start=[v0[i, s] for i =1:nh, s = 1:ne])
-    ExaModels.@add_var(core, w, 1:nh, 1:nc, 1:ne; start=0.0)
+    ExaModels.@add_var(core, w, 1:nh, 1:nc, 1:ne; start=zero_T)
     ExaModels.@add_var(core, uc, 1:nh, 1:nc, 1:ne; start=[v0[i, s] for i =1:nh, j=1:nc, s = 1:ne])
-    ExaModels.@add_var(core, Duc, 1:nh, 1:nc, 1:ne; start=0.0)
+    ExaModels.@add_var(core, Duc, 1:nh, 1:nc, 1:ne; start=zero_T)
 
     itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne]
     itr2 = [(j,rho[j]) for j=1:nc]
 
     # L2 error
-    ExaModels.@add_obj(core, (v[itauj,s] + sum(w[itauj,k,s]*(tauj-tj)^k/(factorial(k)*h^(k-1)) for k in 1:nc) - zjs)^2 for (j,s,itauj, tauj, tj, zjs) in itr)
+    ExaModels.@add_obj(core, (v[itauj,s] + sum(w[itauj,k,s]*(tauj-tj)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - zjs)^2 for (j,s,itauj, tauj, tj, zjs) in itr)
 
     # Collocation model
     ExaModels.@add_con(
         core,
         c1,
-        - uc[i, j, s] + v[i,s] + h*sum(w[i,k,s]*(rhoj^k/factorial(k)) for k in 1:nc)
+        - uc[i, j, s] + v[i,s] + h*sum(w[i,k,s]*(rhoj^k/T(factorial(k))) for k in 1:nc)
         for i=1:nh, (j,rhoj) in itr2, s=1:ne
     )
     ExaModels.@add_con(
         core,
         c2,
-        - Duc[i, j, s] + sum(w[i,k,s]*(rhoj^(k-1)/factorial(k-1)) for k in 1:nc)
+        - Duc[i, j, s] + sum(w[i,k,s]*(rhoj^(k-1)/T(factorial(k-1))) for k in 1:nc)
         for i=1:nh, (j,rhoj) in itr2, s=1:ne
     )
 
@@ -100,7 +101,7 @@
     ExaModels.@add_con(
         core,
         c4,
-        v[i, s] + sum(w[i, j, s]*h/factorial(j) for j in 1:nc) - v[i+1, s]
+        v[i, s] + sum(w[i, j, s]*h/T(factorial(j)) for j in 1:nc) - v[i+1, s]
         for i=1:nh-1, s=1:ne
     )
     ExaModels.@add_con(
@@ -117,4 +118,3 @@
     )
     return ExaModels.ExaModel(core; kwargs...)
 end
-
diff --git a/ext/COPSBenchmarkExaModels/glider.jl b/ext/COPSBenchmarkExaModels/glider.jl
index 322660a..9000e0c 100644
--- a/ext/COPSBenchmarkExaModels/glider.jl
+++ b/ext/COPSBenchmarkExaModels/glider.jl
@@ -6,52 +6,56 @@
 # COPS 3.1 - March 2004
 
 @inline function COPSBenchmark.glider_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...)
-    # Design parameters
-    x_0 = 0.0
-    y_0 = 1000.0
-    y_f = 900.0
-    vx_0 = 13.23
-    vx_f = 13.23
-    vy_0 = -1.288
-    vy_f = -1.288
-    u_c = 2.5
-    r_0 = 100.0
-    m = 100.0
-    g = 9.81
-    cd0 = 0.034
-    cd1 = 0.069662
-    S = 14.0
-    rho = 1.13
-    cL_min = 0.0
-    cL_max = 1.4
-    cL0 = cL_max / 2
+    # Design parameters (T-typed)
+    x_0 = T(0)
+    y_0 = T(1000)
+    y_f = T(900)
+    vx_0 = T(13.23)
+    vx_f = T(13.23)
+    vy_0 = T(-1.288)
+    vy_f = T(-1.288)
+    u_c = T(2.5)
+    r_0 = T(100)
+    m = T(100)
+    g = T(9.81)
+    cd0 = T(0.034)
+    cd1 = T(0.069662)
+    S = T(14)
+    rho = T(1.13)
+    cL_min = T(0)
+    cL_max = T(1.4)
+    cL0 = cL_max / T(2)
+    half = T(0.5)
+    inv_nh = T(1) / T(nh)
+    r_offset = T(2.5)
+    one_T = T(1)
 
     c = ExaModels.ExaCore(T; backend = backend, concrete = Val(true))
 
-    ExaModels.@add_var(c, t_f, 1; lvar = 0.0, start = 1.0)
-    ExaModels.@add_var(c, x, nh+1; lvar = zeros(nh+1), start = [x_0 + vx_0*(k/nh) for k in 0:nh])
-    ExaModels.@add_var(c, y, nh+1; start = [y_0 + (k/nh)*(y_f - y_0) for k in 0:nh])
-    ExaModels.@add_var(c, vx, nh+1; lvar = zeros(nh+1), start = fill(vx_0, nh+1))
+    ExaModels.@add_var(c, t_f, 1; lvar = T(0), start = T(1))
+    ExaModels.@add_var(c, x, nh+1; lvar = zeros(T, nh+1), start = [x_0 + vx_0*(T(k)*inv_nh) for k in 0:nh])
+    ExaModels.@add_var(c, y, nh+1; start = [y_0 + (T(k)*inv_nh)*(y_f - y_0) for k in 0:nh])
+    ExaModels.@add_var(c, vx, nh+1; lvar = zeros(T, nh+1), start = fill(vx_0, nh+1))
     ExaModels.@add_var(c, vy, nh+1; start = fill(vy_0, nh+1))
     ExaModels.@add_var(c, cL, nh+1; lvar = fill(cL_min,nh+1), uvar = fill(cL_max, nh+1), start = fill(cL0, nh+1))
 
     # Expressions (matching JuMP @expressions)
-    ExaModels.@add_expr(c, r, (x[i]/r_0 - 2.5)^2 for i in 1:nh+1)
-    ExaModels.@add_expr(c, u, u_c*(1 - r[i])*exp(-r[i]) for i in 1:nh+1)
+    ExaModels.@add_expr(c, r, (x[i]/r_0 - r_offset)^2 for i in 1:nh+1)
+    ExaModels.@add_expr(c, u, u_c*(one_T - r[i])*exp(-r[i]) for i in 1:nh+1)
     ExaModels.@add_expr(c, w, vy[i] - u[i] for i in 1:nh+1)
     ExaModels.@add_expr(c, v, sqrt(vx[i]^2 + w[i]^2) for i in 1:nh+1)
-    ExaModels.@add_expr(c, D, 0.5*(cd0+cd1*cL[i]^2)*rho*S*v[i]^2 for i in 1:nh+1)
-    ExaModels.@add_expr(c, L, 0.5*cL[i]*rho*S*v[i]^2 for i in 1:nh+1)
+    ExaModels.@add_expr(c, D, half*(cd0+cd1*cL[i]^2)*rho*S*v[i]^2 for i in 1:nh+1)
+    ExaModels.@add_expr(c, L, half*cL[i]*rho*S*v[i]^2 for i in 1:nh+1)
     ExaModels.@add_expr(c, vx_dot, (-L[i]*(w[i]/v[i]) - D[i]*(vx[i]/v[i]))/m for i in 1:nh+1)
     ExaModels.@add_expr(c, vy_dot, (L[i]*(vx[i]/v[i]) - D[i]*(w[i]/v[i]))/m - g for i in 1:nh+1)
 
     ExaModels.@add_obj(c, -x[nh+1])
 
     # Dynamics
-    ExaModels.@add_con(c, c1, x[j] - (x[j-1] + 0.5 * t_f[1]/nh * (vx[j] + vx[j-1])) for j in 2:nh+1)
-    ExaModels.@add_con(c, c2, y[j] - (y[j-1] + 0.5 * t_f[1]/nh * (vy[j] + vy[j-1])) for j in 2:nh+1)
-    ExaModels.@add_con(c, c3, vx[j] - (vx[j-1] + 0.5 * t_f[1]/nh * (vx_dot[j] + vx_dot[j-1])) for j in 2:nh+1)
-    ExaModels.@add_con(c, c4, vy[j] - (vy[j-1] + 0.5 * t_f[1]/nh * (vy_dot[j] + vy_dot[j-1])) for j in 2:nh+1)
+    ExaModels.@add_con(c, c1, x[j] - (x[j-1] + half * t_f[1]*inv_nh * (vx[j] + vx[j-1])) for j in 2:nh+1)
+    ExaModels.@add_con(c, c2, y[j] - (y[j-1] + half * t_f[1]*inv_nh * (vy[j] + vy[j-1])) for j in 2:nh+1)
+    ExaModels.@add_con(c, c3, vx[j] - (vx[j-1] + half * t_f[1]*inv_nh * (vx_dot[j] + vx_dot[j-1])) for j in 2:nh+1)
+    ExaModels.@add_con(c, c4, vy[j] - (vy[j-1] + half * t_f[1]*inv_nh * (vy_dot[j] + vy_dot[j-1])) for j in 2:nh+1)
 
     # Boundary constraints
     ExaModels.@add_con(c, c5, x[1] - x_0)
diff --git a/ext/COPSBenchmarkExaModels/marine.jl b/ext/COPSBenchmarkExaModels/marine.jl
index d297adb..51da740 100644
--- a/ext/COPSBenchmarkExaModels/marine.jl
+++ b/ext/COPSBenchmarkExaModels/marine.jl
@@ -10,14 +10,15 @@
     ne = 8     # number of differential equations
     nm = 21    # number of measurements
 
-    rho = [0.5]                            # roots of k-th degree Legendre polynomial
-    tau = collect(range(0.0, 10.0, 21))    # times at which observations made
-    tf = tau[nm]                           # ODEs defined in [0,tf]
-    h = tf / nh                            # uniform interval length
-    t = [(i-1)*h for i in 1:nh+1]          # partition
+    rho = T[0.5]                                  # roots of k-th degree Legendre polynomial
+    tau = collect(T, range(T(0), T(10), 21))      # times at which observations made
+    tf = tau[nm]                                  # ODEs defined in [0,tf]
+    h = tf / T(nh)                                # uniform interval length
+    t = T[T(i-1)*h for i in 1:nh+1]               # partition
+    zero_T = T(0)
 
     # itau[i] is the largest integer k with t[k] <= tau[i]
-    itau = Int[min(nh, floor(tau[i]/h)+1) for i in 1:nm]
+    itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm]
 
     # Observation
     z = reshape(T[
@@ -44,7 +45,7 @@
         1.0,    12.0,   198.0,   707.0,  2562.0,  3163.0,  3232.0,  5566.0,
     ], ne, nm)'
 
-    v0 = zeros(nh, ne)
+    v0 = zeros(T, nh, ne)
     # Starting-value
     for i in 1:itau[1], s in 1:ne
         v0[i, s] = z[1, s]
@@ -59,15 +60,15 @@
     core = ExaModels.ExaCore(T; backend= backend, concrete = Val(true))
 
     # Growth rates
-    ExaModels.@add_var(core, g, 1:ne-1; lvar =  0.0)
+    ExaModels.@add_var(core, g, 1:ne-1; lvar = zero_T)
     # Mortality rates
-    ExaModels.@add_var(core, m, 1:ne; lvar = 0.0)
+    ExaModels.@add_var(core, m, 1:ne; lvar = zero_T)
     # The collocation approximation u is defined by the parameters v and w.
     # uc and Duc are, respectively, u and u' evaluated at the collocation points.
     ExaModels.@add_var(core, v, 1:nh, 1:ne; start=[v0[i, s] for i=1:nh, s=1:ne])
-    ExaModels.@add_var(core, w, 1:nh, 1:nc, 1:ne; start=0.0)
+    ExaModels.@add_var(core, w, 1:nh, 1:nc, 1:ne; start=zero_T)
     ExaModels.@add_var(core, uc, 1:nh, 1:nc, 1:ne; start=[v0[i, s] for i=1:nh, s=1:ne])
-    ExaModels.@add_var(core, Duc, 1:nh, 1:nc, 1:ne; start=0.0)
+    ExaModels.@add_var(core, Duc, 1:nh, 1:nc, 1:ne; start=zero_T)
 
     # error
 
@@ -76,26 +77,26 @@
     itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne]
     itr2 = [(j,rho[j]) for j=1:nc]
 
-    ExaModels.@add_obj(core, (v[itauj,s] + sum(w[itauj,k,s]*(tauj-tj)^k/(factorial(k)*h^(k-1)) for k in 1:nc) - zjs)^2 for (j,s,itauj,tauj, tj, zjs) in itr)
+    ExaModels.@add_obj(core, (v[itauj,s] + sum(w[itauj,k,s]*(tauj-tj)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - zjs)^2 for (j,s,itauj,tauj, tj, zjs) in itr)
 
     # Collocation model
     ExaModels.@add_con(
         core,
         c1,
-        - uc[i, j, s] + v[i,s] + h*sum(w[i,k,s]*(rhoj^k/factorial(k)) for k in 1:nc)
+        - uc[i, j, s] + v[i,s] + h*sum(w[i,k,s]*(rhoj^k/T(factorial(k))) for k in 1:nc)
         for i=1:nh, (j,rhoj) in itr2, s=1:ne
     )
     ExaModels.@add_con(
         core,
         c2,
-        - Duc[i, j, s] + sum(w[i,k,s]*(rhoj^(k-1)/factorial(k-1)) for k in 1:nc)
+        - Duc[i, j, s] + sum(w[i,k,s]*(rhoj^(k-1)/T(factorial(k-1))) for k in 1:nc)
         for i=1:nh, (j,rhoj) in itr2, s=1:ne
     )
     # Continuity
     ExaModels.@add_con(
         core,
         c3,
-        v[i, s] + sum(w[i, j, s]*h/factorial(j) for j in 1:nc) - v[i+1, s]
+        v[i, s] + sum(w[i, j, s]*h/T(factorial(j)) for j in 1:nc) - v[i+1, s]
         for i=1:nh-1, s=1:ne
     )
     # Boundary conditions
@@ -121,5 +122,3 @@
 
     return ExaModels.ExaModel(core; kwargs...)
 end
-
-
diff --git a/ext/COPSBenchmarkExaModels/methanol.jl b/ext/COPSBenchmarkExaModels/methanol.jl
index 1c49551..4a9835d 100644
--- a/ext/COPSBenchmarkExaModels/methanol.jl
+++ b/ext/COPSBenchmarkExaModels/methanol.jl
@@ -11,9 +11,9 @@
     nc = 3
     nm = 17
 
-    rho = [0.11270166537926, 0.5, 0.88729833462074]
+    rho = T[0.11270166537926, 0.5, 0.88729833462074]
     # times at which observations made
-    tau = [
+    tau = T[
         0.,
         0.050,
         0.065,
@@ -32,12 +32,14 @@
         0.937,
         1.122,
     ]
-    tf = tau[nm]                   # ODEs defined in [0,tf]
-    h = tf / nh                    # uniform interval length
-    t = [(i-1)*h for i in 1:nh+1]  # partition
+    tf = tau[nm]                          # ODEs defined in [0,tf]
+    h = tf / T(nh)                        # uniform interval length
+    t = T[T(i-1)*h for i in 1:nh+1]       # partition
+    zero_T = T(0)
+    two_T = T(2)
 
     # itau[i] is the largest integer k with t[k] <= tau[i]
-    itau = Int[min(nh, floor(tau[i]/h)+1) for i in 1:nm]
+    itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm]
 
     # Concentrations
     z = reshape(T[
@@ -60,10 +62,10 @@
         0.0006, 0.3698, 0.2899,
     ], ne, nm)'
     con1_matrix = [(j, s, itau[j], tau[j], z[j,s], t[itau[j]]) for j in 1:nm, s in 1:ne]
-    bc = [1.0, 0.0, 0.0]
+    bc = T[1, 0, 0]
 
     # Starting-value
-    v0 = zeros(nh, ne)
+    v0 = zeros(T, nh, ne)
     for i in 1:itau[1], s in 1:ne
         v0[i, s] = bc[s]
     end
@@ -73,27 +75,27 @@
     for i in itau[nm]+1:nh, s in 1:ne
         v0[i, s] = z[nm, s]
     end
-    v0 .= 0.001
+    v0 .= T(0.001)
 
     c = ExaModels.ExaCore(T; backend = backend, concrete = Val(true))
-    ExaModels.@add_var(c, theta, np; lvar = 0, start = fill(1, np))
+    ExaModels.@add_var(c, theta, np; lvar = zero_T, start = fill(T(1), np))
     ExaModels.@add_var(c, v, nh, ne; start = v0)
-    ExaModels.@add_var(c, w, nh, nc, ne; start = 0)
+    ExaModels.@add_var(c, w, nh, nc, ne; start = zero_T)
     ExaModels.@add_var(c, uc, nh, nc, ne; start = [v0[i,s] for i=1:nh, j=1:nc, s=1:ne])
-    ExaModels.@add_var(c, Duc, nh, nc, ne; start = 0)
+    ExaModels.@add_var(c, Duc, nh, nc, ne; start = zero_T)
 
-    ExaModels.@add_obj(c, (v[itau,s] + sum(w[itau,k,s]*(tau-t)^k/(factorial(k)*h^(k-1)) for k in 1:nc) - z)^2 for (j,s,itau,tau,z,t) in con1_matrix)
+    ExaModels.@add_obj(c, (v[itau,s] + sum(w[itau,k,s]*(tau-t)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - z)^2 for (j,s,itau,tau,z,t) in con1_matrix)
 
     ExaModels.@add_con(
         c,
         c1,
-        uc[i, j, s] - v[i,s] - h*sum(w[i,k,s]*(rho^k/factorial(k)) for k in 1:nc) for i=1:nh, (j,rho) in [(j, rho[j]) for j in 1:nc], s=1:ne
+        uc[i, j, s] - v[i,s] - h*sum(w[i,k,s]*(rho^k/T(factorial(k))) for k in 1:nc) for i=1:nh, (j,rho) in [(j, rho[j]) for j in 1:nc], s=1:ne
     )
 
     ExaModels.@add_con(
         c,
         c2,
-        Duc[i, j, s] - sum(w[i,k,s]*(rho^(k-1)/factorial(k-1)) for k in 1:nc) for i=1:nh, (j,rho) in [(j, rho[j]) for j in 1:nc], s=1:ne
+        Duc[i, j, s] - sum(w[i,k,s]*(rho^(k-1)/T(factorial(k-1))) for k in 1:nc) for i=1:nh, (j,rho) in [(j, rho[j]) for j in 1:nc], s=1:ne
     )
 
     ExaModels.@add_con(
@@ -106,13 +108,13 @@
     ExaModels.@add_con(
         c,
         c4,
-        v[i, s] + sum(w[i, j, s]*h/factorial(j) for j in 1:nc) - v[i+1, s] for i=1:nh-1, s=1:ne
+        v[i, s] + sum(w[i, j, s]*h/T(factorial(j)) for j in 1:nc) - v[i+1, s] for i=1:nh-1, s=1:ne
     )
 
     ExaModels.@add_con(
         c,
         c5,
-        Duc[i,j,1] + ((2*theta[2] - (theta[1]*uc[i,j,2])/((theta[2]+theta[5])*uc[i,j,1]+uc[i,j,2]) +
+        Duc[i,j,1] + ((two_T*theta[2] - (theta[1]*uc[i,j,2])/((theta[2]+theta[5])*uc[i,j,1]+uc[i,j,2]) +
                          theta[3] + theta[4])*uc[i,j,1]) for i=1:nh, j=1:nc
     )
 
diff --git a/ext/COPSBenchmarkExaModels/pde_models.jl b/ext/COPSBenchmarkExaModels/pde_models.jl
index 645d542..7a08854 100644
--- a/ext/COPSBenchmarkExaModels/pde_models.jl
+++ b/ext/COPSBenchmarkExaModels/pde_models.jl
@@ -1,6 +1,12 @@
 
 @inline function COPSBenchmark.transition_state_model(::ExaModelsBackend, problem, dom::COPSBenchmark.PDEDiscretizationDomain; T = Float64, backend = nothing, kwargs...)
     a, b, c, d, p = problem.a, problem.b, problem.c, problem.d, problem.p
+    # Precompute T-typed constants so the kernel expressions stay fp32-clean:
+    # the scalar `a` and the integer divisors `/2`, `1/(DIMEN+1)` would
+    # otherwise inject Float64 into the GPU AD kernels.
+    aT = T(a)
+    half = T(1) / T(2)
+    inv_dim = T(1) / T(dom.DIMEN + 1)
     x0 = COPSBenchmark._initial_position!(problem, dom, 10)
     array1 = [
         (
@@ -75,7 +81,7 @@
         core,
         c3,
         AREA*(
-            a / (8*AREA^2)*(
+            aT / (8*AREA^2)*(
                 u[b1,TRIANG1]^2*(EDGE_21^2 + EDGE_22^2) +
                 u[b1,TRIANG2]^2*(EDGE_31^2 + EDGE_32^2) +
                 u[b1,TRIANG3]^2*(EDGE_11^2 + EDGE_12^2) +
@@ -91,8 +97,8 @@
     ExaModels.@add_con!(
         core,
         c3,
-        (b1, e1) => AREA* 1 / (dom.DIMEN+1) *
-                    (b*u[b1,TRIANG]^2/2- c*u[b1,TRIANG]^(p+1)/(p+1)+ d*u[b1, TRIANG])
+        (b1, e1) => AREA * inv_dim *
+                    (b*u[b1,TRIANG]^2*half - c*u[b1,TRIANG]^(p+1)/(p+1)+ d*u[b1, TRIANG])
                     for b1 in 1:dom.BREAK+2, (e1, AREA, b, c, d, p, TRIANG) in array2
     )
 
diff --git a/ext/COPSBenchmarkExaModels/pinene.jl b/ext/COPSBenchmarkExaModels/pinene.jl
index 0bde8f2..485666d 100644
--- a/ext/COPSBenchmarkExaModels/pinene.jl
+++ b/ext/COPSBenchmarkExaModels/pinene.jl
@@ -12,17 +12,18 @@
     nm = 8        # number of measurements
 
     # roots of k-th degree Legendre polynomial
-    rho = [0.11270166537926, 0.5, 0.88729833462074]
+    rho = T[0.11270166537926, 0.5, 0.88729833462074]
     # boundary conditions
-    bc = [100.0, 0.0, 0.0, 0.0, 0.0]
+    bc = T[100, 0, 0, 0, 0]
     # times at which observations made
-    tau = [1230.0, 3060.0, 4920.0, 7800.0, 10680.0, 15030.0, 22620.0, 36420.0]
-    tf = tau[nm]                       # ODEs defined in [0,tf]
-    h = tf / nh                        # uniform interval length
-    t = [(i-1)*h for i in 1:nh+1]      # partition
+    tau = T[1230, 3060, 4920, 7800, 10680, 15030, 22620, 36420]
+    tf = tau[nm]                            # ODEs defined in [0,tf]
+    h = tf / T(nh)                          # uniform interval length
+    t = T[T(i-1)*h for i in 1:nh+1]         # partition
+    zero_T = T(0)
 
     # itau[i] is the largest integer k with t[k] <= tau[i]
-    itau = Int[min(nh, floor(tau[i]/h)+1) for i in 1:nm]
+    itau = Int[min(nh, Int(floor(tau[i]/h))+1) for i in 1:nm]
 
     # Observations
     z = reshape(T[
@@ -36,7 +37,7 @@
         4.5,  63.1,  3.8,  2.9, 25.7,
     ], ne, nm)'
 
-    v0 = zeros(nh, ne)
+    v0 = zeros(T, nh, ne)
     # Starting-value
     for i in 1:itau[1], s in 1:ne
         v0[i, s] = bc[s]
@@ -50,30 +51,30 @@
 
     core = ExaModels.ExaCore(T; backend= backend, concrete = Val(true))
 
-    ExaModels.@add_var(core, theta, 1:np; lvar = 0.0, start=0.0)
+    ExaModels.@add_var(core, theta, 1:np; lvar = zero_T, start=zero_T)
     # The collocation approximation u is defined by the parameters v and w.
     # uc and Duc are, respectively, u and u' evaluated at the collocation points.
     ExaModels.@add_var(core, v, 1:nh, 1:ne; start=[v0[i, s] for i=1:nh, s=1:ne])
-    ExaModels.@add_var(core, w, 1:nh, 1:nc, 1:ne; start=0.0)
+    ExaModels.@add_var(core, w, 1:nh, 1:nc, 1:ne; start=zero_T)
     ExaModels.@add_var(core, uc, 1:nh, 1:nc, 1:ne; start=[v0[i,s] for i=1:nh, j=1:nc, s=1:ne])
-    ExaModels.@add_var(core, Duc, 1:nh, 1:nc, 1:ne; start=0.0)
+    ExaModels.@add_var(core, Duc, 1:nh, 1:nc, 1:ne; start=zero_T)
 
     itr = [(j, s, itau[j], tau[j], t[itau[j]], z[j,s]) for j=1:nm, s in 1:ne]
     # l2 error
-    ExaModels.@add_obj(core, (v[it,s] + sum(w[it,k,s]*(tj-ti)^k/(factorial(k)*h^(k-1)) for k in 1:nc) - zjs)^2 for (j, s, it, tj, ti, zjs) in itr)
+    ExaModels.@add_obj(core, (v[it,s] + sum(w[it,k,s]*(tj-ti)^k/(T(factorial(k))*h^(k-1)) for k in 1:nc) - zjs)^2 for (j, s, it, tj, ti, zjs) in itr)
 
     # Collocation model
     itr2 = [(j,rho[j]) for j=1:nc]
     ExaModels.@add_con(
         core,
         c1,
-        - uc[i, j, s] + v[i,s] + h*sum(w[i,k,s]*(rhoj^k/factorial(k)) for k in 1:nc)
+        - uc[i, j, s] + v[i,s] + h*sum(w[i,k,s]*(rhoj^k/T(factorial(k))) for k in 1:nc)
         for i=1:nh, (j,rhoj) in itr2, s=1:ne
     )
     ExaModels.@add_con(
         core,
         c2,
-        - Duc[i, j, s] + sum(w[i,k,s]*(rhoj^(k-1)/factorial(k-1)) for k in 1:nc)
+        - Duc[i, j, s] + sum(w[i,k,s]*(rhoj^(k-1)/T(factorial(k-1))) for k in 1:nc)
         for i=1:nh, (j,rhoj) in itr2, s=1:ne
     )
     # Boundary
@@ -82,7 +83,7 @@
     ExaModels.@add_con(
         core,
         c4,
-        v[i, s] + sum(w[i, j, s]*h/factorial(j) for j in 1:nc) - v[i+1, s]
+        v[i, s] + sum(w[i, j, s]*h/T(factorial(j)) for j in 1:nc) - v[i+1, s]
         for i=1:nh-1, s=1:ne
     )
     ExaModels.@add_con(core, c5, -Duc[i,j,1] - (theta[1]+theta[2])*uc[i,j,1] for i=1:nh, j=1:nc)
@@ -93,5 +94,3 @@
 
     return ExaModels.ExaModel(core; kwargs...)
 end
-
-
diff --git a/ext/COPSBenchmarkExaModels/robot.jl b/ext/COPSBenchmarkExaModels/robot.jl
index 4794f2c..dcbf135 100644
--- a/ext/COPSBenchmarkExaModels/robot.jl
+++ b/ext/COPSBenchmarkExaModels/robot.jl
@@ -8,59 +8,66 @@
 @inline function COPSBenchmark.robot_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...)
 
     # total length of arm
-    L = 5.0
+    L = T(5)
 
     # Upper bounds on the controls
-    max_u_rho = 1.0
-    max_u_the = 1.0
-    max_u_phi = 1.0
+    max_u_rho = T(1)
+    max_u_the = T(1)
+    max_u_phi = T(1)
     # Initial positions of the length and the angles for the robot arm
-    rho0 = 4.5
-    phi0 = pi /4
+    rho0 = T(4.5)
+    pi_T = T(pi)
+    phi0 = pi_T / T(4)
+    half = T(0.5)
+    third = T(1) / T(3)
+    inv_nh = T(1) / T(nh)
+    two_pi_3 = T(2) * pi_T / T(3)
+    four_pi_3 = T(4) * pi_T / T(3)
+    zero_T = T(0)
 
     core = ExaModels.ExaCore(T; backend= backend, concrete = Val(true))
 
-    ExaModels.@add_var(core, rho, nh+1; start=rho0, lvar =  0 , uvar = L)
-    ExaModels.@add_var(core, the, nh+1; start=[2*pi/3*(k/nh)^2 for k=1:nh+1], lvar =  -pi , uvar = pi)
-    ExaModels.@add_var(core, phi, nh+1; start=phi0, lvar = 0, uvar = pi)
+    ExaModels.@add_var(core, rho, nh+1; start=rho0, lvar = zero_T, uvar = L)
+    ExaModels.@add_var(core, the, nh+1; start=[two_pi_3*(T(k)*inv_nh)^2 for k=1:nh+1], lvar = -pi_T, uvar = pi_T)
+    ExaModels.@add_var(core, phi, nh+1; start=phi0, lvar = zero_T, uvar = pi_T)
     # Derivatives
-    ExaModels.@add_var(core, rho_dot, nh+1; start=0.0)
-    ExaModels.@add_var(core, the_dot, nh+1; start=[4*pi/3*(k/nh) for k=1:nh+1])
-    ExaModels.@add_var(core, phi_dot, nh+1; start=0.0)
+    ExaModels.@add_var(core, rho_dot, nh+1; start=zero_T)
+    ExaModels.@add_var(core, the_dot, nh+1; start=[four_pi_3*(T(k)*inv_nh) for k=1:nh+1])
+    ExaModels.@add_var(core, phi_dot, nh+1; start=zero_T)
     # Control
-    ExaModels.@add_var(core, u_rho, nh+1; start=0.0, lvar = -max_u_rho, uvar = max_u_rho)
-    ExaModels.@add_var(core, u_the, nh+1; start=0.0, lvar = -max_u_the, uvar = max_u_the)
-    ExaModels.@add_var(core, u_phi, nh+1; start=0.0, lvar = -max_u_phi, uvar = max_u_phi)
+    ExaModels.@add_var(core, u_rho, nh+1; start=zero_T, lvar = -max_u_rho, uvar = max_u_rho)
+    ExaModels.@add_var(core, u_the, nh+1; start=zero_T, lvar = -max_u_the, uvar = max_u_the)
+    ExaModels.@add_var(core, u_phi, nh+1; start=zero_T, lvar = -max_u_phi, uvar = max_u_phi)
     # Final time
-    ExaModels.@add_var(core, tf, 1; start = 1.0, lvar = 0.0)
+    ExaModels.@add_var(core, tf, 1; start = T(1), lvar = zero_T)
 
     # Expressions (matching JuMP @expressions)
-    ExaModels.@add_expr(core, I_the, ((L-rho[i])^3+rho[i]^3)*(sin(phi[i]))^2/3.0 for i=1:nh+1)
-    ExaModels.@add_expr(core, I_phi, ((L-rho[i])^3+rho[i]^3)/3.0 for i=1:nh+1)
+    ExaModels.@add_expr(core, I_the, ((L-rho[i])^3+rho[i]^3)*(sin(phi[i]))^2*third for i=1:nh+1)
+    ExaModels.@add_expr(core, I_phi, ((L-rho[i])^3+rho[i]^3)*third for i=1:nh+1)
 
     ExaModels.@add_obj(core, tf[1] for i=1:1)
 
     # Dynamics
-    ExaModels.@add_con(core, c1, - rho[j] + rho[j-1] + 0.5 * tf[1]/nh * (rho_dot[j] + rho_dot[j-1]) for j=2:nh+1)
-    ExaModels.@add_con(core, c2, - phi[j] + phi[j-1] + 0.5 * tf[1]/nh * (phi_dot[j] + phi_dot[j-1]) for j=2:nh+1)
-    ExaModels.@add_con(core, c3, - the[j] + the[j-1] + 0.5 * tf[1]/nh * (the_dot[j] + the_dot[j-1]) for j=2:nh+1)
-    ExaModels.@add_con(core, c4, - rho_dot[j] + rho_dot[j-1] + 0.5 * tf[1]/nh * (u_rho[j] + u_rho[j-1]) / L for j=2:nh+1)
-    ExaModels.@add_con(core, c5, - the_dot[j] + the_dot[j-1] + 0.5 * tf[1]/nh * (u_the[j] / I_the[j] + u_the[j-1] / I_the[j-1]) for j=2:nh+1)
-    ExaModels.@add_con(core, c6, - phi_dot[j] + phi_dot[j-1] + 0.5 * tf[1]/nh * (u_phi[j] / I_phi[j] + u_phi[j-1] / I_phi[j-1]) for j=2:nh+1)
+    ExaModels.@add_con(core, c1, - rho[j] + rho[j-1] + half * tf[1]*inv_nh * (rho_dot[j] + rho_dot[j-1]) for j=2:nh+1)
+    ExaModels.@add_con(core, c2, - phi[j] + phi[j-1] + half * tf[1]*inv_nh * (phi_dot[j] + phi_dot[j-1]) for j=2:nh+1)
+    ExaModels.@add_con(core, c3, - the[j] + the[j-1] + half * tf[1]*inv_nh * (the_dot[j] + the_dot[j-1]) for j=2:nh+1)
+    ExaModels.@add_con(core, c4, - rho_dot[j] + rho_dot[j-1] + half * tf[1]*inv_nh * (u_rho[j] + u_rho[j-1]) / L for j=2:nh+1)
+    ExaModels.@add_con(core, c5, - the_dot[j] + the_dot[j-1] + half * tf[1]*inv_nh * (u_the[j] / I_the[j] + u_the[j-1] / I_the[j-1]) for j=2:nh+1)
+    ExaModels.@add_con(core, c6, - phi_dot[j] + phi_dot[j-1] + half * tf[1]*inv_nh * (u_phi[j] / I_phi[j] + u_phi[j-1] / I_phi[j-1]) for j=2:nh+1)
 
     # Boundary conditions
-    ExaModels.@add_con(core, c7, - rho[1] + 4.5)
-    ExaModels.@add_con(core, c8, - the[1] + 0.0)
-    ExaModels.@add_con(core, c9, - phi[1] + pi / 4.0)
-    ExaModels.@add_con(core, c10, - rho[nh+1] + 4.5)
-    ExaModels.@add_con(core, c11, - the[nh+1] + 2.0 * pi / 3)
-    ExaModels.@add_con(core, c12, - phi[nh+1] + pi / 4.0)
-    ExaModels.@add_con(core, c13, - rho_dot[1] + 0.0)
-    ExaModels.@add_con(core, c14, - the_dot[1] + 0.0)
-    ExaModels.@add_con(core, c15, - phi_dot[1] + 0.0)
-    ExaModels.@add_con(core, c16, - rho_dot[nh+1] + 0.0)
-    ExaModels.@add_con(core, c17, - the_dot[nh+1] + 0.0)
-    ExaModels.@add_con(core, c18, - phi_dot[nh+1] + 0.0)
+    ExaModels.@add_con(core, c7, - rho[1] + rho0)
+    ExaModels.@add_con(core, c8, - the[1] + zero_T)
+    ExaModels.@add_con(core, c9, - phi[1] + phi0)
+    ExaModels.@add_con(core, c10, - rho[nh+1] + rho0)
+    ExaModels.@add_con(core, c11, - the[nh+1] + two_pi_3)
+    ExaModels.@add_con(core, c12, - phi[nh+1] + phi0)
+    ExaModels.@add_con(core, c13, - rho_dot[1] + zero_T)
+    ExaModels.@add_con(core, c14, - the_dot[1] + zero_T)
+    ExaModels.@add_con(core, c15, - phi_dot[1] + zero_T)
+    ExaModels.@add_con(core, c16, - rho_dot[nh+1] + zero_T)
+    ExaModels.@add_con(core, c17, - the_dot[nh+1] + zero_T)
+    ExaModels.@add_con(core, c18, - phi_dot[nh+1] + zero_T)
 
     return ExaModels.ExaModel(core; kwargs...)
 end
diff --git a/ext/COPSBenchmarkExaModels/rocket.jl b/ext/COPSBenchmarkExaModels/rocket.jl
index 5e7ed22..9b3fcaf 100644
--- a/ext/COPSBenchmarkExaModels/rocket.jl
+++ b/ext/COPSBenchmarkExaModels/rocket.jl
@@ -5,34 +5,43 @@
 # COPS 3.1 - March 2004
 
 @inline function COPSBenchmark.rocket_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...)
-    h_0 = 1.0
-    v_0 = 0.0
-    m_0 = 1.0
-    g_0 = 1.0
-    T_c = 3.5
-    h_c = 500.0
-    v_c = 620.0
-    m_c = 0.6
-
-    c = 0.5*sqrt(g_0 * h_0)
+    h_0 = T(1)
+    v_0 = T(0)
+    m_0 = T(1)
+    g_0 = T(1)
+    T_c = T(3.5)
+    h_c = T(500)
+    v_c = T(620)
+    m_c = T(0.6)
+    half = T(0.5)
+    inv_nh = T(1) / T(nh)
+
+    c = half*sqrt(g_0 * h_0)
     m_f = m_c * m_0
-    D_c = 0.5 * v_c * (m_0 / g_0)
+    D_c = half * v_c * (m_0 / g_0)
     T_max = T_c * m_0 * g_0
 
     core = ExaModels.ExaCore(T; backend= backend, minimize=false, concrete = Val(true))
 
-    ExaModels.@add_var(core, h, 0:nh; start=1.0, lvar = 1.0)
-    ExaModels.@add_var(core, v, 0:nh; start=(i/nh*(1.0 - i/nh) for i=0:nh), lvar = 0.0)
-    ExaModels.@add_var(core, m, 0:nh; start=((m_f - m_0)*(i/nh) + m_0 for i=0:nh), lvar = m_f, uvar = m_0)
-    ExaModels.@add_var(core, T, 0:nh; start=T_max/2.0, lvar = 0.0, uvar = T_max)
-    ExaModels.@add_var(core, step, 1; start=1/nh, lvar = 0.0)
+    # Precompute generator scalars so the broadcast closure stays isbits
+    # (Metal rejects Type{T} captures). And bind the thrust-variable name as
+    # Th rather than T to avoid shadowing the Type parameter T inside the
+    # function — that shadowing caused Julia to Box the captured T value.
+    s_v_start  = T[T(i)*inv_nh*(T(1) - T(i)*inv_nh) for i=0:nh]
+    s_m_start  = T[(m_f - m_0)*(T(i)*inv_nh) + m_0 for i=0:nh]
+    s_Th_start = T_max*half
+    ExaModels.@add_var(core, h, 0:nh; start=h_0, lvar = h_0)
+    ExaModels.@add_var(core, v, 0:nh; start=s_v_start, lvar = v_0)
+    ExaModels.@add_var(core, m, 0:nh; start=s_m_start, lvar = m_f, uvar = m_0)
+    ExaModels.@add_var(core, Th, 0:nh; start=s_Th_start, lvar = v_0, uvar = T_max)
+    ExaModels.@add_var(core, step, 1; start=inv_nh, lvar = v_0)
 
     ExaModels.@add_obj(core, h[nh])
 
     # Dynamics
-    ExaModels.@add_con(core, c1, - h[i] + h[i-1] + 0.5 * step[1] * (v[i] + v[i-1]) for i=1:nh)
-    ExaModels.@add_con(core, c2, - v[i] + v[i-1] + 0.5 * step[1] * ((T[i] - D_c*v[i]^2*exp(-h_c*(h[i] - h_0))/h_0 - m[i] * g_0 * (h_0 / h[i])^2) / m[i] + (T[i-1] - D_c*v[i-1]^2*exp(-h_c*(h[i-1] - h_0))/h_0 - m[i-1] * g_0 * (h_0 / h[i-1])^2) / m[i-1]) for i=1:nh)
-    ExaModels.@add_con(core, c3, - m[i] + m[i-1] + 0.5 * step[1] * (-T[i]/c + -T[i-1]/c) for i=1:nh)
+    ExaModels.@add_con(core, c1, - h[i] + h[i-1] + half * step[1] * (v[i] + v[i-1]) for i=1:nh)
+    ExaModels.@add_con(core, c2, - v[i] + v[i-1] + half * step[1] * ((Th[i] - D_c*v[i]^2*exp(-h_c*(h[i] - h_0))/h_0 - m[i] * g_0 * (h_0 / h[i])^2) / m[i] + (Th[i-1] - D_c*v[i-1]^2*exp(-h_c*(h[i-1] - h_0))/h_0 - m[i-1] * g_0 * (h_0 / h[i-1])^2) / m[i-1]) for i=1:nh)
+    ExaModels.@add_con(core, c3, - m[i] + m[i-1] + half * step[1] * (-Th[i]/c + -Th[i-1]/c) for i=1:nh)
 
     # Boundary ExaModels.constraints
     ExaModels.@add_con(core, c4, h[0] - h_0)
@@ -42,4 +51,3 @@
 
     return ExaModels.ExaModel(core; kwargs...)
 end
-
diff --git a/ext/COPSBenchmarkExaModels/steering.jl b/ext/COPSBenchmarkExaModels/steering.jl
index 8930c1d..e28f8a2 100644
--- a/ext/COPSBenchmarkExaModels/steering.jl
+++ b/ext/COPSBenchmarkExaModels/steering.jl
@@ -5,43 +5,43 @@
 # COPS 3.1 - March 2004
 
 @inline function COPSBenchmark.steering_model(::ExaModelsBackend, nh; T = Float64, backend = nothing, kwargs...)
-    a = 100.0  # Magnitude of force.
+    a = T(100)  # Magnitude of force.
     # Bounds on the control
-    u_min, u_max = -pi/2.0, pi/2.0
-    xs = zeros(4)
-    xf = [NaN, 5.0, 45.0, 0.0]
+    u_min, u_max = -T(pi)/T(2), T(pi)/T(2)
+    xs = zeros(T, 4)
+    xf = [T(NaN), T(5), T(45), T(0)]
+    half = T(0.5)
+    inv_nh = T(1) / T(nh)
 
     function gen_x0(k, i)
         if i == 1 || i == 4
-            return 0.0
+            return T(0)
         elseif i == 2
-            return 5*k/nh
+            return T(5)*T(k)*inv_nh
         elseif i == 3
-            return 45.0*k/nh
+            return T(45)*T(k)*inv_nh
         else
-            return 0.0
+            return T(0)
         end
     end
 
     core = ExaModels.ExaCore(T; backend= backend, concrete = Val(true))
 
-    ExaModels.@add_var(core, u, 1:nh+1; lvar = u_min, uvar =  u_max, start=0.0)   # control
+    ExaModels.@add_var(core, u, 1:nh+1; lvar = u_min, uvar =  u_max, start=T(0))   # control
     ExaModels.@add_var(core, x, 1:nh+1, 1:4; start=[gen_x0(i, j) for i=1:nh+1, j=1:4])     # state
-    ExaModels.@add_var(core, tf, 1; start=1.0)                 # final time
+    ExaModels.@add_var(core, tf, 1; start=T(1))                 # final time
 
     ExaModels.@add_obj(core, tf[1])
 
-    ExaModels.@add_con(core, c0, tf[1]; lcon = 0., ucon= Inf)
+    ExaModels.@add_con(core, c0, tf[1]; lcon = T(0), ucon= T(Inf))
     # Dynamics
-    ExaModels.@add_con(core, c1, -x[i+1,1] + x[i,1] + 0.5*(tf[1] / nh)*(x[i,3] + x[i+1,3]) for i=1:nh)
-    ExaModels.@add_con(core, c2, -x[i+1,2] + x[i,2] + 0.5*(tf[1] / nh)*(x[i,4] + x[i+1,4]) for i=1:nh)
-    ExaModels.@add_con(core, c3, -x[i+1,3] + x[i,3] + 0.5*(tf[1] / nh)*(a*cos(u[i]) + a*cos(u[i+1])) for i=1:nh)
-    ExaModels.@add_con(core, c4, -x[i+1,4] + x[i,4] + 0.5*(tf[1] / nh)*(a*sin(u[i]) + a*sin(u[i+1])) for i=1:nh)
+    ExaModels.@add_con(core, c1, -x[i+1,1] + x[i,1] + half*(tf[1]*inv_nh)*(x[i,3] + x[i+1,3]) for i=1:nh)
+    ExaModels.@add_con(core, c2, -x[i+1,2] + x[i,2] + half*(tf[1]*inv_nh)*(x[i,4] + x[i+1,4]) for i=1:nh)
+    ExaModels.@add_con(core, c3, -x[i+1,3] + x[i,3] + half*(tf[1]*inv_nh)*(a*cos(u[i]) + a*cos(u[i+1])) for i=1:nh)
+    ExaModels.@add_con(core, c4, -x[i+1,4] + x[i,4] + half*(tf[1]*inv_nh)*(a*sin(u[i]) + a*sin(u[i+1])) for i=1:nh)
     # Boundary conditions
     ExaModels.@add_con(core, c5, -x[1, j] + s for (j,s) in enumerate(xs))
     ExaModels.@add_con(core, c6, -x[nh+1, j] + f for (j,f) in zip(2:4, xf[2:4]))
 
     return ExaModels.ExaModel(core; kwargs...)
 end
-
-
diff --git a/test/runtests.jl b/test/runtests.jl
index e47d30f..bbc6891 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -135,8 +135,8 @@ end
     @test m == get_ncon(exa_model)
 
     # Find variable permutation by matching (lvar, uvar, x0) tuples
-    lj, uj, x0j = get_lvar(jump_nlp), get_uvar(jump_nlp), jump_nlp.meta.x0
-    le, ue, x0e = get_lvar(exa_model), get_uvar(exa_model), exa_model.meta.x0
+    lj, uj, x0j = NLPModels.get_lvar(jump_nlp), NLPModels.get_uvar(jump_nlp), jump_nlp.meta.x0
+    le, ue, x0e = NLPModels.get_lvar(exa_model), NLPModels.get_uvar(exa_model), exa_model.meta.x0
     var_keys_j = collect(zip(lj, uj, x0j))
     var_keys_e = collect(zip(le, ue, x0e))
     Pv = find_permutation(var_keys_j, var_keys_e)
@@ -157,8 +157,8 @@ end
     @test obj(jump_nlp, x0j) ≈ obj(exa_model, x0e) rtol = 1e-6
 
     # Find constraint permutation
-    lcj, ucj = get_lcon(jump_nlp), get_ucon(jump_nlp)
-    lce, uce = get_lcon(exa_model), get_ucon(exa_model)
+    lcj, ucj = NLPModels.get_lcon(jump_nlp), NLPModels.get_ucon(jump_nlp)
+    lce, uce = NLPModels.get_lcon(exa_model), NLPModels.get_ucon(exa_model)
     cj0 = cons(jump_nlp, x0j)
     ce0 = cons(exa_model, x0e)
     con_keys_j = collect(zip(lcj, ucj, cj0))