diff --git a/src/exa_linalg.jl b/src/exa_linalg.jl
index fd3647d..e86c84d 100644
--- a/src/exa_linalg.jl
+++ b/src/exa_linalg.jl
@@ -32,7 +32,8 @@ import Base: sinh, cosh, tanh, csch, sech, coth, asinh, acosh, atanh, acoth
 import Base: ^, zeros, ones
 import LinearAlgebra: dot, Adjoint, det, tr, norm, diag, diagm
 
-export zero, one, zeros, ones, adjoint, transpose, *, +, -, sum, dot, det, tr, norm, diag, diagm
+export zero,
+    one, zeros, ones, adjoint, transpose, *, +, -, sum, dot, det, tr, norm, diag, diagm
 export inv, abs, sqrt, cbrt, abs2, exp, exp2, exp10, log, log2, log10, log1p
 export sin, cos, tan, csc, sec, cot, asin, acos, atan, acot
 export sind, cosd, tand, cscd, secd, cotd, atand, acotd
@@ -60,9 +61,13 @@ const DenseSubArray{T,N} = SubArray{T,N,<:DenseArray}
 # ReshapedArray where parent is a DenseArray (excludes sparse arrays)
 const DenseReshapedArray{T,N} = Base.ReshapedArray{T,N,<:DenseArray}
 
-const VecReal{T<:Real} = Union{Vector{T}, DenseSubArray{T,1}, DenseReshapedArray{T,1}, Base.ReinterpretArray{T,1}}
-const VecNode{T<:ExaModels.AbstractNode} = Union{Vector{T}, DenseSubArray{T,1}, DenseReshapedArray{T,1}}
-const Mat{T} = Union{Matrix{T}, DenseSubArray{T,2}, DenseReshapedArray{T,2}}
+const VecReal{T<:Real} = Union{
+    Vector{T},DenseSubArray{T,1},DenseReshapedArray{T,1},Base.ReinterpretArray{T,1}
+}
+const VecNode{T<:ExaModels.AbstractNode} = Union{
+    Vector{T},DenseSubArray{T,1},DenseReshapedArray{T,1}
+}
+const Mat{T} = Union{Matrix{T},DenseSubArray{T,2},DenseReshapedArray{T,2}}
 
 # ============================================================================
 # Section 1: Canonical Nodes (zero and one)
@@ -111,7 +116,7 @@ one(::ExaModels.AbstractNode) = ExaModels.Null(1)
 Create an array of AbstractNode zeros with the specified dimensions.
 Uses fill with the canonical zero node: Null(0).
 """
-function zeros(::Type{T}, dims::Integer...) where {T <: ExaModels.AbstractNode}
+function zeros(::Type{T}, dims::Integer...) where {T<:ExaModels.AbstractNode}
     return fill(zero(T), dims...)
 end
 
@@ -121,7 +126,7 @@ end
 Create an array of AbstractNode ones with the specified dimensions.
 Uses fill with the canonical one node: Null(1).
 """
-function ones(::Type{T}, dims::Integer...) where {T <: ExaModels.AbstractNode}
+function ones(::Type{T}, dims::Integer...) where {T<:ExaModels.AbstractNode}
     return fill(one(T), dims...)
 end
 
@@ -152,7 +157,9 @@ end
 # ============================================================================
 
 # Addition: Null{T} + Null{S} → Null
-+(x::ExaModels.Null{T}, y::ExaModels.Null{S}) where {T<:Real, S<:Real} = ExaModels.Null(x.value + y.value)
+function +(x::ExaModels.Null{T}, y::ExaModels.Null{S}) where {T<:Real,S<:Real}
+    ExaModels.Null(x.value + y.value)
+end
 # Addition: Null{T} + AbstractNode → unwrap Null, use native +
 +(x::ExaModels.Null{T}, y::ExaModels.AbstractNode) where {T<:Real} = x.value + y
 +(x::ExaModels.AbstractNode, y::ExaModels.Null{T}) where {T<:Real} = x + y.value
@@ -161,19 +168,29 @@ end
 +(x::Real, y::ExaModels.Null{T}) where {T<:Real} = ExaModels.Null(x + y.value)
 
 # Subtraction: Null{T} - Null{S} → Null
--(x::ExaModels.Null{T}, y::ExaModels.Null{S}) where {T<:Real, S<:Real} = ExaModels.Null(x.value - y.value)
+function -(x::ExaModels.Null{T}, y::ExaModels.Null{S}) where {T<:Real,S<:Real}
+    ExaModels.Null(x.value - y.value)
+end
 # Subtraction: Null{T} - AbstractNode → handle 0 - e = -e specially
--(x::ExaModels.Null{T}, y::ExaModels.AbstractNode) where {T<:Real} = iszero(x.value) ? (-y) : (x.value - y)
+function -(x::ExaModels.Null{T}, y::ExaModels.AbstractNode) where {T<:Real}
+    iszero(x.value) ? (-y) : (x.value - y)
+end
 -(x::ExaModels.AbstractNode, y::ExaModels.Null{T}) where {T<:Real} = x - y.value
 # Subtraction: Null{T} - Real → Null (more specific than ExaModels' AbstractNode - Real)
 -(x::ExaModels.Null{T}, y::Real) where {T<:Real} = ExaModels.Null(x.value - y)
 -(x::Real, y::ExaModels.Null{T}) where {T<:Real} = ExaModels.Null(x - y.value)
 
 # Multiplication: Null{T} * Null{S} → Null
-*(x::ExaModels.Null{T}, y::ExaModels.Null{S}) where {T<:Real, S<:Real} = ExaModels.Null(x.value * y.value)
+function *(x::ExaModels.Null{T}, y::ExaModels.Null{S}) where {T<:Real,S<:Real}
+    ExaModels.Null(x.value * y.value)
+end
 # Multiplication: Null{T} * AbstractNode → zero optimization: 0 * e = Null(0)
-*(x::ExaModels.Null{T}, y::ExaModels.AbstractNode) where {T<:Real} = iszero(x.value) ? ExaModels.Null(0) : (x.value * y)
-*(x::ExaModels.AbstractNode, y::ExaModels.Null{T}) where {T<:Real} = iszero(y.value) ? ExaModels.Null(0) : (x * y.value)
+function *(x::ExaModels.Null{T}, y::ExaModels.AbstractNode) where {T<:Real}
+    iszero(x.value) ? ExaModels.Null(0) : (x.value * y)
+end
+function *(x::ExaModels.AbstractNode, y::ExaModels.Null{T}) where {T<:Real}
+    iszero(y.value) ? ExaModels.Null(0) : (x * y.value)
+end
 # Multiplication: Null{T} * Real → zero optimization, more specific than ExaModels' AbstractNode * Real
 *(x::ExaModels.Null{T}, y::Real) where {T<:Real} = ExaModels.Null(x.value * y)
 *(x::Real, y::ExaModels.Null{T}) where {T<:Real} = ExaModels.Null(x * y.value)
@@ -272,7 +289,9 @@ end
 adjoint(x::ExaModels.AbstractNode) = x
 transpose(x::ExaModels.AbstractNode) = x
 
-convert(::Type{ExaModels.AbstractNode}, x::Real) = iszero(x) ? zero(ExaModels.AbstractNode) : ExaModels.Null(x)
+function convert(::Type{ExaModels.AbstractNode}, x::Real)
+    iszero(x) ? zero(ExaModels.AbstractNode) : ExaModels.Null(x)
+end
 
 promote_rule(::Type{<:ExaModels.AbstractNode}, ::Type{<:Real}) = ExaModels.AbstractNode
 
@@ -284,17 +303,19 @@ promote_rule(::Type{<:ExaModels.AbstractNode}, ::Type{<:Real}) = ExaModels.Abstr
 # __dot handles wrapping Real values in Null nodes internally.
 # ============================================================================
 
-function dot(v::VecReal{T}, w::VecNode{S}) where {T<:Real, S<:ExaModels.AbstractNode}
+function dot(v::VecReal{T}, w::VecNode{S}) where {T<:Real,S<:ExaModels.AbstractNode}
     @assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
     return __dot(v, w)
 end
 
-function dot(v::VecNode{T}, w::VecReal{S}) where {T<:ExaModels.AbstractNode, S<:Real}
+function dot(v::VecNode{T}, w::VecReal{S}) where {T<:ExaModels.AbstractNode,S<:Real}
     @assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
     return __dot(v, w)
 end
 
-function dot(v::VecNode{T}, w::VecNode{S}) where {T<:ExaModels.AbstractNode, S<:ExaModels.AbstractNode}
+function dot(
+    v::VecNode{T}, w::VecNode{S}
+) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     @assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
     return __dot(v, w)
 end
@@ -311,19 +332,25 @@ end
 # This eliminates redundant checks in hot loops (e.g., 100 checks for 100×100 matrix×vector).
 # ============================================================================
 
-function __dot(v::AbstractVector{T}, w::AbstractVector{S}) where {T<:Real, S<:Real} # fallback to ensure __dot also accepts vectors of reals
+function __dot(v::AbstractVector{T}, w::AbstractVector{S}) where {T<:Real,S<:Real} # fallback to ensure __dot also accepts vectors of reals
     return sum(v[i] * w[i] for i in eachindex(v))
 end
 
-function __dot(v::AbstractVector{T}, w::AbstractVector{S}) where {T<:Real, S<:ExaModels.AbstractNode}
+function __dot(
+    v::AbstractVector{T}, w::AbstractVector{S}
+) where {T<:Real,S<:ExaModels.AbstractNode}
     return sum(ExaModels.Null(v[i]) * w[i] for i in eachindex(v))
 end
 
-function __dot(v::AbstractVector{T}, w::AbstractVector{S}) where {T<:ExaModels.AbstractNode, S<:Real}
+function __dot(
+    v::AbstractVector{T}, w::AbstractVector{S}
+) where {T<:ExaModels.AbstractNode,S<:Real}
     return sum(v[i] * ExaModels.Null(w[i]) for i in eachindex(v))
 end
 
-function __dot(v::AbstractVector{T}, w::AbstractVector{S}) where {T<:ExaModels.AbstractNode, S<:ExaModels.AbstractNode}
+function __dot(
+    v::AbstractVector{T}, w::AbstractVector{S}
+) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     return sum(v[i] * w[i] for i in eachindex(v))
 end
 
@@ -336,54 +363,54 @@ end
 # ============================================================================
 
 # Scalar × Vector
-function *(a::T, v::VecReal{<:Real}) where {T <: ExaModels.AbstractNode}
+function *(a::T, v::VecReal{<:Real}) where {T<:ExaModels.AbstractNode}
     return [a * ExaModels.Null(vi) for vi in v]
 end
 
-function *(a::Real, v::VecNode{T}) where {T <: ExaModels.AbstractNode}
+function *(a::Real, v::VecNode{T}) where {T<:ExaModels.AbstractNode}
     return [ExaModels.Null(a) * vi for vi in v]
 end
 
-function *(a::T, v::VecNode{S}) where {T <: ExaModels.AbstractNode, S <: ExaModels.AbstractNode}
+function *(a::T, v::VecNode{S}) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     return [a * vi for vi in v]
 end
 
 # Vector × Scalar
-function *(v::VecNode{T}, a::Real) where {T <: ExaModels.AbstractNode}
+function *(v::VecNode{T}, a::Real) where {T<:ExaModels.AbstractNode}
     return [vi * ExaModels.Null(a) for vi in v]
 end
 
-function *(v::VecReal{T}, a::S) where {T <: Real, S <: ExaModels.AbstractNode}
+function *(v::VecReal{T}, a::S) where {T<:Real,S<:ExaModels.AbstractNode}
     return [ExaModels.Null(vi) * a for vi in v]
 end
 
-function *(v::VecNode{T}, a::S) where {T <: ExaModels.AbstractNode, S <: ExaModels.AbstractNode}
+function *(v::VecNode{T}, a::S) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     return [vi * a for vi in v]
 end
 
 # Scalar × Matrix
-function *(a::T, A::Mat{<:Real}) where {T <: ExaModels.AbstractNode}
+function *(a::T, A::Mat{<:Real}) where {T<:ExaModels.AbstractNode}
     return [a * ExaModels.Null(A[i, j]) for i in axes(A, 1), j in axes(A, 2)]
 end
 
-function *(a::Real, A::Mat{T}) where {T <: ExaModels.AbstractNode}
+function *(a::Real, A::Mat{T}) where {T<:ExaModels.AbstractNode}
     return [ExaModels.Null(a) * A[i, j] for i in axes(A, 1), j in axes(A, 2)]
 end
 
-function *(a::T, A::Mat{S}) where {T <: ExaModels.AbstractNode, S <: ExaModels.AbstractNode}
+function *(a::T, A::Mat{S}) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     return [a * A[i, j] for i in axes(A, 1), j in axes(A, 2)]
 end
 
 # Matrix × Scalar
-function *(A::Mat{T}, a::Real) where {T <: ExaModels.AbstractNode}
+function *(A::Mat{T}, a::Real) where {T<:ExaModels.AbstractNode}
     return [A[i, j] * ExaModels.Null(a) for i in axes(A, 1), j in axes(A, 2)]
 end
 
-function *(A::Mat{T}, a::S) where {T <: Real, S <: ExaModels.AbstractNode}
+function *(A::Mat{T}, a::S) where {T<:Real,S<:ExaModels.AbstractNode}
     return [ExaModels.Null(A[i, j]) * a for i in axes(A, 1), j in axes(A, 2)]
 end
 
-function *(A::Mat{T}, a::S) where {T <: ExaModels.AbstractNode, S <: ExaModels.AbstractNode}
+function *(A::Mat{T}, a::S) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     return [A[i, j] * a for i in axes(A, 1), j in axes(A, 2)]
 end
 
@@ -391,19 +418,21 @@ end
 # Section 6: Matrix × Vector Product (uses dot)
 # ============================================================================
 
-function *(A::Mat{<:Real}, x::VecNode{T}) where {T <: ExaModels.AbstractNode}
+function *(A::Mat{<:Real}, x::VecNode{T}) where {T<:ExaModels.AbstractNode}
     m, n = size(A)
     @assert n == length(x) "Dimension mismatch: matrix has $n columns but vector has $(length(x)) elements"
     return [__dot(view(A, i, :), x) for i in 1:m]
 end
 
-function *(A::Mat{T}, x::VecReal{<:Real}) where {T <: ExaModels.AbstractNode}
+function *(A::Mat{T}, x::VecReal{<:Real}) where {T<:ExaModels.AbstractNode}
     m, n = size(A)
     @assert n == length(x) "Dimension mismatch: matrix has $n columns but vector has $(length(x)) elements"
     return [__dot(view(A, i, :), x) for i in 1:m]
 end
 
-function *(A::Mat{T}, x::VecNode{S}) where {T <: ExaModels.AbstractNode, S <: ExaModels.AbstractNode}
+function *(
+    A::Mat{T}, x::VecNode{S}
+) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     m, n = size(A)
     @assert n == length(x) "Dimension mismatch: matrix has $n columns but vector has $(length(x)) elements"
     return [__dot(view(A, i, :), x) for i in 1:m]
@@ -413,21 +442,21 @@ end
 # Section 7: Matrix × Matrix Product (uses dot)
 # ============================================================================
 
-function *(A::Mat{<:Real}, B::Mat{T}) where {T <: ExaModels.AbstractNode}
+function *(A::Mat{<:Real}, B::Mat{T}) where {T<:ExaModels.AbstractNode}
     m, n = size(A)
     p, q = size(B)
     @assert n == p "Dimension mismatch: A has $n columns but B has $p rows"
     return [__dot(view(A, i, :), view(B, :, j)) for i in 1:m, j in 1:q]
 end
 
-function *(A::Mat{T}, B::Mat{<:Real}) where {T <: ExaModels.AbstractNode}
+function *(A::Mat{T}, B::Mat{<:Real}) where {T<:ExaModels.AbstractNode}
     m, n = size(A)
     p, q = size(B)
     @assert n == p "Dimension mismatch: A has $n columns but B has $p rows"
     return [__dot(view(A, i, :), view(B, :, j)) for i in 1:m, j in 1:q]
 end
 
-function *(A::Mat{T}, B::Mat{S}) where {T <: ExaModels.AbstractNode, S <: ExaModels.AbstractNode}
+function *(A::Mat{T}, B::Mat{S}) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     m, n = size(A)
     p, q = size(B)
     @assert n == p "Dimension mismatch: A has $n columns but B has $p rows"
@@ -438,21 +467,23 @@ end
 # Section 8: Adjoint Vector × Matrix Product
 # ============================================================================
 
-function *(p::Adjoint{T, <:VecNode{T}}, A::Mat{<:Real}) where {T <: ExaModels.AbstractNode}
+function *(p::Adjoint{T,<:VecNode{T}}, A::Mat{<:Real}) where {T<:ExaModels.AbstractNode}
     m, n = size(A)
     @assert m == length(p) "Dimension mismatch: vector has $(length(p)) elements but matrix has $m rows"
     v = parent(p)
     return [__dot(v, view(A, :, j)) for j in 1:n]'
 end
 
-function *(p::Adjoint{T, <:VecReal{T}}, A::Mat{S}) where {T <: Real, S <: ExaModels.AbstractNode}
+function *(p::Adjoint{T,<:VecReal{T}}, A::Mat{S}) where {T<:Real,S<:ExaModels.AbstractNode}
     m, n = size(A)
     @assert m == length(p) "Dimension mismatch: vector has $(length(p)) elements but matrix has $m rows"
     v = parent(p)
     return [__dot(v, view(A, :, j)) for j in 1:n]'
 end
 
-function *(p::Adjoint{T, <:VecNode{T}}, A::Mat{S}) where {T <: ExaModels.AbstractNode, S <: ExaModels.AbstractNode}
+function *(
+    p::Adjoint{T,<:VecNode{T}}, A::Mat{S}
+) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     m, n = size(A)
     @assert m == length(p) "Dimension mismatch: vector has $(length(p)) elements but matrix has $m rows"
     v = parent(p)
@@ -463,17 +494,23 @@ end
 # Section 8.5: Adjoint Vector × Vector Product
 # ============================================================================
 
-function *(p::Adjoint{T, <:VecNode{T}}, w::VecReal{S}) where {T<:ExaModels.AbstractNode, S<:Real}
+function *(
+    p::Adjoint{T,<:VecNode{T}}, w::VecReal{S}
+) where {T<:ExaModels.AbstractNode,S<:Real}
     @assert length(p) == length(w) "Vectors must have the same length: got $(length(p)) and $(length(w))"
     return __dot(parent(p), w)
 end
 
-function *(p::Adjoint{T, <:VecReal{T}}, w::VecNode{S}) where {T<:Real, S<:ExaModels.AbstractNode}
+function *(
+    p::Adjoint{T,<:VecReal{T}}, w::VecNode{S}
+) where {T<:Real,S<:ExaModels.AbstractNode}
     @assert length(p) == length(w) "Vectors must have the same length: got $(length(p)) and $(length(w))"
     return __dot(parent(p), w)
 end
 
-function *(p::Adjoint{T, <:VecNode{T}}, w::VecNode{S}) where {T<:ExaModels.AbstractNode, S<:ExaModels.AbstractNode}
+function *(
+    p::Adjoint{T,<:VecNode{T}}, w::VecNode{S}
+) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     @assert length(p) == length(w) "Vectors must have the same length: got $(length(p)) and $(length(w))"
     return __dot(parent(p), w)
 end
@@ -482,11 +519,11 @@ end
 # Section 9: Adjoint and Transpose for Matrices
 # ============================================================================
 
-function adjoint(A::Mat{T}) where {T <: ExaModels.AbstractNode}
+function adjoint(A::Mat{T}) where {T<:ExaModels.AbstractNode}
     return permutedims(A)
 end
 
-function transpose(A::Mat{T}) where {T <: ExaModels.AbstractNode}
+function transpose(A::Mat{T}) where {T<:ExaModels.AbstractNode}
     return permutedims(A)
 end
 
@@ -495,33 +532,35 @@ end
 # ============================================================================
 
 # Vector + Vector
-function +(v::VecNode{T}, w::VecReal{<:Real}) where {T <: ExaModels.AbstractNode}
+function +(v::VecNode{T}, w::VecReal{<:Real}) where {T<:ExaModels.AbstractNode}
     @assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
     return [v[i] + w[i] for i in eachindex(v)]
 end
 
-function +(v::VecReal{<:Real}, w::VecNode{T}) where {T <: ExaModels.AbstractNode}
+function +(v::VecReal{<:Real}, w::VecNode{T}) where {T<:ExaModels.AbstractNode}
     @assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
     return [v[i] + w[i] for i in eachindex(v)]
 end
 
-function +(v::VecNode{T}, w::VecNode{S}) where {T <: ExaModels.AbstractNode, S <: ExaModels.AbstractNode}
+function +(
+    v::VecNode{T}, w::VecNode{S}
+) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     @assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
     return [v[i] + w[i] for i in eachindex(v)]
 end
 
 # Matrix + Matrix
-function +(A::Mat{T}, B::Mat{<:Real}) where {T <: ExaModels.AbstractNode}
+function +(A::Mat{T}, B::Mat{<:Real}) where {T<:ExaModels.AbstractNode}
     @assert size(A) == size(B) "Matrices must have the same size: got $(size(A)) and $(size(B))"
     return [A[i, j] + B[i, j] for i in axes(A, 1), j in axes(A, 2)]
 end
 
-function +(A::Mat{<:Real}, B::Mat{T}) where {T <: ExaModels.AbstractNode}
+function +(A::Mat{<:Real}, B::Mat{T}) where {T<:ExaModels.AbstractNode}
     @assert size(A) == size(B) "Matrices must have the same size: got $(size(A)) and $(size(B))"
     return [A[i, j] + B[i, j] for i in axes(A, 1), j in axes(A, 2)]
 end
 
-function +(A::Mat{T}, B::Mat{S}) where {T <: ExaModels.AbstractNode, S <: ExaModels.AbstractNode}
+function +(A::Mat{T}, B::Mat{S}) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     @assert size(A) == size(B) "Matrices must have the same size: got $(size(A)) and $(size(B))"
     return [A[i, j] + B[i, j] for i in axes(A, 1), j in axes(A, 2)]
 end
@@ -531,33 +570,35 @@ end
 # ============================================================================
 
 # Vector - Vector
-function -(v::VecNode{T}, w::VecReal{<:Real}) where {T <: ExaModels.AbstractNode}
+function -(v::VecNode{T}, w::VecReal{<:Real}) where {T<:ExaModels.AbstractNode}
     @assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
     return [v[i] - w[i] for i in eachindex(v)]
 end
 
-function -(v::VecReal{<:Real}, w::VecNode{T}) where {T <: ExaModels.AbstractNode}
+function -(v::VecReal{<:Real}, w::VecNode{T}) where {T<:ExaModels.AbstractNode}
     @assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
     return [v[i] - w[i] for i in eachindex(v)]
 end
 
-function -(v::VecNode{T}, w::VecNode{S}) where {T <: ExaModels.AbstractNode, S <: ExaModels.AbstractNode}
+function -(
+    v::VecNode{T}, w::VecNode{S}
+) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     @assert length(v) == length(w) "Vectors must have the same length: got $(length(v)) and $(length(w))"
     return [v[i] - w[i] for i in eachindex(v)]
 end
 
 # Matrix - Matrix
-function -(A::Mat{T}, B::Mat{<:Real}) where {T <: ExaModels.AbstractNode}
+function -(A::Mat{T}, B::Mat{<:Real}) where {T<:ExaModels.AbstractNode}
     @assert size(A) == size(B) "Matrices must have the same size: got $(size(A)) and $(size(B))"
     return [A[i, j] - B[i, j] for i in axes(A, 1), j in axes(A, 2)]
 end
 
-function -(A::Mat{<:Real}, B::Mat{T}) where {T <: ExaModels.AbstractNode}
+function -(A::Mat{<:Real}, B::Mat{T}) where {T<:ExaModels.AbstractNode}
     @assert size(A) == size(B) "Matrices must have the same size: got $(size(A)) and $(size(B))"
     return [A[i, j] - B[i, j] for i in axes(A, 1), j in axes(A, 2)]
 end
 
-function -(A::Mat{T}, B::Mat{S}) where {T <: ExaModels.AbstractNode, S <: ExaModels.AbstractNode}
+function -(A::Mat{T}, B::Mat{S}) where {T<:ExaModels.AbstractNode,S<:ExaModels.AbstractNode}
     @assert size(A) == size(B) "Matrices must have the same size: got $(size(A)) and $(size(B))"
     return [A[i, j] - B[i, j] for i in axes(A, 1), j in axes(A, 2)]
 end
@@ -590,7 +631,7 @@ function _det_impl(A)
         # Laplace expansion for n×n matrices (n ≥ 4)
         d = A[1, 1] * det(A[2:end, 2:end])  # Initialize with first term
         for j in 2:n
-            minor = [A[2:end, 1:j-1] A[2:end, j+1:end]]
+            minor = [A[2:end, 1:(j - 1)] A[2:end, (j + 1):end]]
             sign_coeff = iseven(j) ? -1 : 1
             cofactor = sign_coeff * A[1, j] * det(minor)
             d = d + cofactor
@@ -600,14 +641,14 @@ function _det_impl(A)
 end
 
 # Separate methods for each wrapper type
-function det(A::Matrix{T}) where {T <: ExaModels.AbstractNode}
+function det(A::Matrix{T}) where {T<:ExaModels.AbstractNode}
     return _det_impl(A)
 end
 
 # Note: We don't define det for SubArray to avoid potential ambiguities.
 # Users should collect/copy the view first if needed.
 
-function det(A::DenseReshapedArray{T,2}) where {T <: ExaModels.AbstractNode}
+function det(A::DenseReshapedArray{T,2}) where {T<:ExaModels.AbstractNode}
     return _det_impl(A)
 end
 
@@ -616,7 +657,7 @@ end
 # ============================================================================
 
 # Need separate methods to avoid ambiguity with LinearAlgebra.tr(::StridedMatrix)
-function tr(A::Matrix{T}) where {T <: ExaModels.AbstractNode}
+function tr(A::Matrix{T}) where {T<:ExaModels.AbstractNode}
     n, m = size(A)
     @assert n == m "Trace is only defined for square matrices, got $(n)×$(m)"
     return sum(A[i, i] for i in 1:n)
@@ -625,7 +666,7 @@ end
 # Note: We don't define tr for SubArray because it creates unavoidable ambiguity
 # with LinearAlgebra.tr(::StridedMatrix). Users should collect/copy the view first.
 
-function tr(A::DenseReshapedArray{T,2}) where {T <: ExaModels.AbstractNode}
+function tr(A::DenseReshapedArray{T,2}) where {T<:ExaModels.AbstractNode}
     n, m = size(A)
     @assert n == m "Trace is only defined for square matrices, got $(n)×$(m)"
     return sum(A[i, i] for i in 1:n)
@@ -636,12 +677,12 @@ end
 # ============================================================================
 
 # Euclidean norm (2-norm) for vectors
-function norm(v::VecNode{T}) where {T <: ExaModels.AbstractNode}
+function norm(v::VecNode{T}) where {T<:ExaModels.AbstractNode}
     return sqrt(sum(vi * vi for vi in v))
 end
 
 # p-norm for vectors
-function norm(v::VecNode{T}, p::Real) where {T <: ExaModels.AbstractNode}
+function norm(v::VecNode{T}, p::Real) where {T<:ExaModels.AbstractNode}
     if p == Inf
         # Infinity norm: max|vᵢ|
         error("Infinity norm not supported for symbolic AbstractNode vectors")
@@ -658,7 +699,7 @@ function norm(v::VecNode{T}, p::Real) where {T <: ExaModels.AbstractNode}
 end
 
 # Frobenius norm for matrices
-function norm(A::Mat{T}) where {T <: ExaModels.AbstractNode}
+function norm(A::Mat{T}) where {T<:ExaModels.AbstractNode}
     return sqrt(sum(A[i, j] * A[i, j] for i in axes(A, 1), j in axes(A, 2)))
 end
 
@@ -667,14 +708,14 @@ end
 # ============================================================================
 
 # Extract diagonal from matrix
-function diag(A::Mat{T}) where {T <: ExaModels.AbstractNode}
+function diag(A::Mat{T}) where {T<:ExaModels.AbstractNode}
     n, m = size(A)
     k = min(n, m)
     return [A[i, i] for i in 1:k]
 end
 
 # Create diagonal matrix from vector
-function diagm(v::VecNode{T}) where {T <: ExaModels.AbstractNode}
+function diagm(v::VecNode{T}) where {T<:ExaModels.AbstractNode}
     n = length(v)
     # Create a matrix with AbstractNode element type to allow mixed Null types
     D = Matrix{ExaModels.AbstractNode}(undef, n, n)
@@ -685,7 +726,7 @@ function diagm(v::VecNode{T}) where {T <: ExaModels.AbstractNode}
 end
 
 # diagm with pairs (more general form)
-function diagm(kv::Pair{<:Integer, <:VecNode{T}}) where {T <: ExaModels.AbstractNode}
+function diagm(kv::Pair{<:Integer,<:VecNode{T}}) where {T<:ExaModels.AbstractNode}
     k, v = kv
     n = length(v) + abs(k)
     # Create a matrix with AbstractNode element type to allow mixed Null types
diff --git a/src/onepass.jl b/src/onepass.jl
index b6a0ce7..7755d4a 100644
--- a/src/onepass.jl
+++ b/src/onepass.jl
@@ -703,15 +703,15 @@ function p_constraint_exa!(p, p_ocp, e1, e2, e3, c_type, label)
             e2 = replace_call(e2, p.x, p.t0, x0)
             e2 = replace_call(e2, p.x, p.tf, xf)
             e2 = subs2(e2, x0, p.x, 0)
-            e2 = subs(e2, x0, :([$(p.x)[$k, 0] for $k ∈ 1:$(p.dim_x)]))
+            e2 = subs(e2, x0, :([$(p.x)[$k, 0] for $k in 1:($(p.dim_x))]))
             e2 = subs2(e2, xf, p.x, :grid_size)
-            e2 = subs(e2, xf, :([$(p.x)[$k, grid_size] for $k ∈ 1:$(p.dim_x)]))
+            e2 = subs(e2, xf, :([$(p.x)[$k, grid_size] for $k in 1:($(p.dim_x))]))
             quote
                 length($e1) == length($e3) || throw("wrong bound dimension") # (vs. __throw) since raised at runtime
                 if length($e1) == 1
                     $pref.constraint($p_ocp, $e2; lcon=($e1[1]), ucon=($e3[1])) # todo: add _denull
                 else
-                    for $l ∈ 1:length($e1)
+                    for $l in 1:length($e1)
                         $pref.constraint($p_ocp, $e2[$l]; lcon=($e1[$l]), ucon=($e3[$l])) # todo: add _denull
                     end
                 end
@@ -815,17 +815,24 @@ function p_constraint_exa!(p, p_ocp, e1, e2, e3, c_type, label)
             l = __symgen(:l)
             e2 = replace_call(e2, [p.x, p.u], p.t, [xt, ut])
             e2 = subs2(e2, xt, p.x, j)
-            e2 = subs(e2, xt, :([$(p.x)[$k, $j] for $k ∈ 1:$(p.dim_x)]))
+            e2 = subs(e2, xt, :([$(p.x)[$k, $j] for $k in 1:($(p.dim_x))]))
             e2 = subs2(e2, ut, p.u, j)
-            e2 = subs(e2, ut, :([$(p.u)[$k, $j] for $k ∈ 1:$(p.dim_u)]))
+            e2 = subs(e2, ut, :([$(p.u)[$k, $j] for $k in 1:($(p.dim_u))]))
             e2 = subs(e2, p.t, :($(p.t0) + $j * $(p.dt)))
             quote
                 length($e1) == length($e3) || throw("wrong bound dimension") # (vs. __throw) since raised at runtime
                 if length($e1) == 1
-                    $pref.constraint($p_ocp, $e2 for $j in 0:grid_size; lcon=($e1[1]), ucon=($e3[1])) # todo: add _denull
+                    $pref.constraint(
+                        $p_ocp, $e2 for $j in 0:grid_size; lcon=($e1[1]), ucon=($e3[1])
+                    ) # todo: add _denull
                 else
-                    for $l ∈ 1:length($e1)
-                        $pref.constraint($p_ocp, $e2[$l] for $j in 0:grid_size; lcon=($e1[$l]), ucon=($e3[$l])) # todo: add _denull
+                    for $l in 1:length($e1)
+                        $pref.constraint(
+                            $p_ocp,
+                            $e2[$l] for $j in 0:grid_size;
+                            lcon=($e1[$l]),
+                            ucon=($e3[$l]),
+                        ) # todo: add _denull
                     end
                 end
             end
@@ -882,33 +889,47 @@ function p_dynamics_exa!(p, p_ocp, x, t, e)
     j12 = :($j1 + 0.5)
     k = __symgen(:k)
     ej1 = subs2(e, xt, p.x, j1)
-    ej1 = subs(ej1, xt, :([$(p.x)[$k, $j1] for $k ∈ 1:$(p.dim_x)]))
+    ej1 = subs(ej1, xt, :([$(p.x)[$k, $j1] for $k in 1:($(p.dim_x))]))
     ej1 = subs2(ej1, ut, p.u, j1)
-    ej1 = subs(ej1, ut, :([$(p.u)[$k, $j1] for $k ∈ 1:$(p.dim_u)]))
+    ej1 = subs(ej1, ut, :([$(p.u)[$k, $j1] for $k in 1:($(p.dim_u))]))
     ej1 = subs(ej1, p.t, :($(p.t0) + $j1 * $(p.dt)))
     ej2 = subs2(e, xt, p.x, j2)
-    ej2 = subs(ej2, xt, :([$(p.x)[$k, $j2] for $k ∈ 1:$(p.dim_x)]))
+    ej2 = subs(ej2, xt, :([$(p.x)[$k, $j2] for $k in 1:($(p.dim_x))]))
     ej2 = subs2(ej2, ut, p.u, j2)
-    ej2 = subs(ej2, ut, :([$(p.u)[$k, $j2] for $k ∈ 1:$(p.dim_u)]))
+    ej2 = subs(ej2, ut, :([$(p.u)[$k, $j2] for $k in 1:($(p.dim_u))]))
     ej2 = subs(ej2, p.t, :($(p.t0) + $j2 * $(p.dt)))
     ej12 = subs2m(e, xt, p.x, j1)
-    ej12 = subs(ej12, xt, :([(($(p.x)[$k, $j1] + $(p.x)[$k, $j1 + 1]) / 2) for $k ∈ 1:$(p.dim_x)]))
+    ej12 = subs(
+        ej12, xt, :([(($(p.x)[$k, $j1] + $(p.x)[$k, $j1 + 1]) / 2) for $k in 1:($(p.dim_x))
+    ])
+    )
     ej12 = subs2(ej12, ut, p.u, j1)
-    ej12 = subs(ej12, ut, :([$(p.u)[$k, $j1] for $k ∈ 1:$(p.dim_u)]))
+    ej12 = subs(ej12, ut, :([$(p.u)[$k, $j1] for $k in 1:($(p.dim_u))]))
     ej12 = subs(ej12, p.t, :($(p.t0) + $j12 * $(p.dt)))
-    dxj = :([$(p.x)[$k, $j2] - $(p.x)[$k, $j1] for $k ∈ 1:$(p.dim_x)])
+    dxj = :([$(p.x)[$k, $j2] - $(p.x)[$k, $j1] for $k in 1:($(p.dim_x))])
     i = __symgen(:i)
     code = quote
-        for $i ∈ 1:$(p.dim_x)
+        for $i in 1:($(p.dim_x))
             $(p.dyn_con)[$i] = if scheme == :euler # dyn_con already defined outside try catch
-                $pref.constraint($p_ocp, $dxj[$i] - $(p.dt) * $ej1[$i] for $j1 in 0:grid_size-1) # todo: add _denull
+                $pref.constraint(
+                    $p_ocp,
+                    $dxj[$i] - $(p.dt) * $ej1[$i] for $j1 in 0:(grid_size - 1)
+                ) # todo: add _denull
             elseif scheme ∈ (:euler_implicit, :euler_b) # euler_b is deprecated
-                $pref.constraint($p_ocp, $dxj[$i] - $(p.dt) * $ej2[$i] for $j1 in 0:grid_size-1) # todo: add _denull
+                $pref.constraint(
+                    $p_ocp,
+                    $dxj[$i] - $(p.dt) * $ej2[$i] for $j1 in 0:(grid_size - 1)
+                ) # todo: add _denull
             elseif scheme == :midpoint
-                $pref.constraint($p_ocp, $dxj[$i] - $(p.dt) * $ej12[$i] for $j1 in 0:grid_size-1) # todo: add _denull
+                $pref.constraint(
+                    $p_ocp,
+                    $dxj[$i] - $(p.dt) * $ej12[$i] for $j1 in 0:(grid_size - 1)
+                ) # todo: add _denull
             elseif scheme ∈ (:trapeze, :trapezoidal) # trapezoidal is deprecated
                 $pref.constraint(
-                    $p_ocp, $dxj[$i] - $(p.dt) * ($ej1[$i] + $ej2[$i]) / 2 for $j1 in 0:grid_size-1 # todo: add _denull
+                    $p_ocp,
+                    $dxj[$i] - $(p.dt) * ($ej1[$i] + $ej2[$i]) / 2 for
+                    $j1 in 0:(grid_size - 1) # todo: add _denull
                 )
             else
                 throw(
@@ -962,7 +983,8 @@ end
 
 function p_dynamics_coord_exa!(p, p_ocp, x, i::Integer, t, e) # todo: also also add coord = range for :exa
     pref = prefix_exa()
-    i ∈ p.dyn_coords && return __throw("dynamics coordinate $i already defined", p.lnum, p.line)
+    i ∈ p.dyn_coords &&
+        return __throw("dynamics coordinate $i already defined", p.lnum, p.line)
     append!(p.dyn_coords, i)
     xt = __symgen(:xt)
     ut = __symgen(:ut)
@@ -972,31 +994,35 @@ function p_dynamics_coord_exa!(p, p_ocp, x, i::Integer, t, e) # todo: also also
     j12 = :($j1 + 0.5)
     k = __symgen(:k)
     ej1 = subs2(e, xt, p.x, j1)
-    ej1 = subs(ej1, xt, :([$(p.x)[$k, $j1] for $k ∈ 1:$(p.dim_x)]))
+    ej1 = subs(ej1, xt, :([$(p.x)[$k, $j1] for $k in 1:($(p.dim_x))]))
     ej1 = subs2(ej1, ut, p.u, j1)
-    ej1 = subs(ej1, ut, :([$(p.u)[$k, $j1] for $k ∈ 1:$(p.dim_u)]))
+    ej1 = subs(ej1, ut, :([$(p.u)[$k, $j1] for $k in 1:($(p.dim_u))]))
     ej1 = subs(ej1, p.t, :($(p.t0) + $j1 * $(p.dt)))
     ej2 = subs2(e, xt, p.x, j2)
-    ej2 = subs(ej2, xt, :([$(p.x)[$k, $j2] for $k ∈ 1:$(p.dim_x)]))
+    ej2 = subs(ej2, xt, :([$(p.x)[$k, $j2] for $k in 1:($(p.dim_x))]))
     ej2 = subs2(ej2, ut, p.u, j2)
-    ej2 = subs(ej2, ut, :([$(p.u)[$k, $j2] for $k ∈ 1:$(p.dim_u)]))
+    ej2 = subs(ej2, ut, :([$(p.u)[$k, $j2] for $k in 1:($(p.dim_u))]))
     ej2 = subs(ej2, p.t, :($(p.t0) + $j2 * $(p.dt)))
     ej12 = subs2m(e, xt, p.x, j1)
-    ej12 = subs(ej12, xt, :([(($(p.x)[$k, $j1] + $(p.x)[$k, $j1 + 1]) / 2) for $k ∈ 1:$(p.dim_x)]))
+    ej12 = subs(
+        ej12, xt, :([(($(p.x)[$k, $j1] + $(p.x)[$k, $j1 + 1]) / 2) for $k in 1:($(p.dim_x))
+    ])
+    )
     ej12 = subs2(ej12, ut, p.u, j1)
-    ej12 = subs(ej12, ut, :([$(p.u)[$k, $j1] for $k ∈ 1:$(p.dim_u)]))
+    ej12 = subs(ej12, ut, :([$(p.u)[$k, $j1] for $k in 1:($(p.dim_u))]))
     ej12 = subs(ej12, p.t, :($(p.t0) + $j12 * $(p.dt)))
     dxij = :($(p.x)[$i, $j2] - $(p.x)[$i, $j1])
     code = quote
         $(p.dyn_con)[$i] = if scheme == :euler # dyn_con already defined outside try catch
-            $pref.constraint($p_ocp, $dxij - $(p.dt) * $ej1 for $j1 in 0:grid_size-1) # todo: add _denull
+            $pref.constraint($p_ocp, $dxij - $(p.dt) * $ej1 for $j1 in 0:(grid_size - 1)) # todo: add _denull
         elseif scheme ∈ (:euler_implicit, :euler_b) # euler_b is deprecated
-            $pref.constraint($p_ocp, $dxij - $(p.dt) * $ej2 for $j1 in 0:grid_size-1) # todo: add _denull
+            $pref.constraint($p_ocp, $dxij - $(p.dt) * $ej2 for $j1 in 0:(grid_size - 1)) # todo: add _denull
         elseif scheme == :midpoint
-            $pref.constraint($p_ocp, $dxij - $(p.dt) * $ej12 for $j1 in 0:grid_size-1) # todo: add _denull
+            $pref.constraint($p_ocp, $dxij - $(p.dt) * $ej12 for $j1 in 0:(grid_size - 1)) # todo: add _denull
         elseif scheme ∈ (:trapeze, :trapezoidal) # trapezoidal is deprecated
             $pref.constraint(
-                $p_ocp, $dxij - $(p.dt) * ($ej1 + $ej2) / 2 for $j1 in 0:grid_size-1 # todo: add _denull
+                $p_ocp,
+                $dxij - $(p.dt) * ($ej1 + $ej2) / 2 for $j1 in 0:(grid_size - 1) # todo: add _denull
             )
         else
             throw(
@@ -1046,25 +1072,28 @@ function p_lagrange_exa!(p, p_ocp, e, type)
     j12 = :($j1 + 0.5)
     k = __symgen(:k)
     ej1 = subs2(e, xt, p.x, j1)
-    ej1 = subs(ej1, xt, :([$(p.x)[$k, $j1] for $k ∈ 1:$(p.dim_x)]))
+    ej1 = subs(ej1, xt, :([$(p.x)[$k, $j1] for $k in 1:($(p.dim_x))]))
     ej1 = subs2(ej1, ut, p.u, j1)
-    ej1 = subs(ej1, ut, :([$(p.u)[$k, $j1] for $k ∈ 1:$(p.dim_u)]))
+    ej1 = subs(ej1, ut, :([$(p.u)[$k, $j1] for $k in 1:($(p.dim_u))]))
     ej1 = subs(ej1, p.t, :($(p.t0) + $j1 * $(p.dt)))
     ej12 = subs2m(e, xt, p.x, j1)
-    ej12 = subs(ej12, xt, :([(($(p.x)[$k, $j1] + $(p.x)[$k, $j1 + 1]) / 2) for $k ∈ 1:$(p.dim_x)]))
+    ej12 = subs(
+        ej12, xt, :([(($(p.x)[$k, $j1] + $(p.x)[$k, $j1 + 1]) / 2) for $k in 1:($(p.dim_x))
+    ])
+    )
     ej12 = subs2(ej12, ut, p.u, j1)
-    ej12 = subs(ej12, ut, :([$(p.u)[$k, $j1] for $k ∈ 1:$(p.dim_u)]))
+    ej12 = subs(ej12, ut, :([$(p.u)[$k, $j1] for $k in 1:($(p.dim_u))]))
     ej12 = subs(ej12, p.t, :($(p.t0) + $j12 * $(p.dt)))
     code = quote
         if scheme == :euler
-            $pref.objective($p_ocp, $(p.dt) * $ej1 for $j1 in 0:grid_size-1) # todo: add _denull
+            $pref.objective($p_ocp, $(p.dt) * $ej1 for $j1 in 0:(grid_size - 1)) # todo: add _denull
         elseif scheme ∈ (:euler_implicit, :euler_b) # euler_b is deprecated
             $pref.objective($p_ocp, $(p.dt) * $ej1 for $j1 in 1:grid_size) # todo: add _denull
         elseif scheme == :midpoint
-            $pref.objective($p_ocp, $(p.dt) * $ej12 for $j1 in 0:grid_size-1) # todo: add _denull
+            $pref.objective($p_ocp, $(p.dt) * $ej12 for $j1 in 0:(grid_size - 1)) # todo: add _denull
         elseif scheme ∈ (:trapeze, :trapezoidal) # trapezoidal is deprecated
             $pref.objective($p_ocp, $(p.dt) * $ej1 / 2 for $j1 in (0, grid_size)) # todo: add _denull
-            $pref.objective($p_ocp, $(p.dt) * $ej1 for $j1 in 1:grid_size-1) # todo: add _denull
+            $pref.objective($p_ocp, $(p.dt) * $ej1 for $j1 in 1:(grid_size - 1)) # todo: add _denull
         else
             throw(
                 "unknown numerical scheme: $scheme (possible choices are :euler, :euler_implicit, :midpoint, :trapeze)",
@@ -1114,9 +1143,9 @@ function p_mayer_exa!(p, p_ocp, e, type)
     e = replace_call(e, p.x, p.t0, x0)
     e = replace_call(e, p.x, p.tf, xf)
     e = subs2(e, x0, p.x, 0)
-    e = subs(e, x0, :([$(p.x)[$k, 0] for $k ∈ 1:$(p.dim_x)]))
+    e = subs(e, x0, :([$(p.x)[$k, 0] for $k in 1:($(p.dim_x))]))
     e = subs2(e, xf, p.x, :grid_size)
-    e = subs(e, xf, :([$(p.x)[$k, grid_size] for $k ∈ 1:$(p.dim_x)]))
+    e = subs(e, xf, :([$(p.x)[$k, grid_size] for $k in 1:($(p.dim_x))]))
     # now, x[i](t0) has been replaced by x[i, 0] and x[i](tf) by x[i, grid_size]
     code = :($pref.objective($p_ocp, $e)) # todo: add _denull
     return __wrap(code, p.lnum, p.line)
@@ -1370,9 +1399,12 @@ function def_exa(e; log=false)
     p = ParsingInfo()
     code = parse!(p, p_ocp, e; log=log, backend=:exa)
     dyn_check = quote
-        !$p.is_global_dyn && !$p.is_coord_dyn && throw($e_pref.ParsingError("dynamics not defined")) # not $(p.xxxx) as these infos are known statically
+        !$p.is_global_dyn &&
+            !$p.is_coord_dyn &&
+            throw($e_pref.ParsingError("dynamics not defined")) # not $(p.xxxx) as these infos are known statically
         if $p.is_coord_dyn # same: also known statically
-            sort($(p.dyn_coords)) == 1:($(p.dim_x)) || throw($e_pref.ParsingError("some coordinates of dynamics undefined"))
+            sort($(p.dyn_coords)) == 1:($(p.dim_x)) ||
+                throw($e_pref.ParsingError("some coordinates of dynamics undefined"))
         end
     end
     default_scheme = QuoteNode(__default_scheme_exa())
@@ -1434,7 +1466,9 @@ function def_exa(e; log=false)
             $(p.box_u) # lvar and uvar for control
             $(p.box_v) # lvar and uvar for variable (after x and u for compatibility with CTDirect)
             $p_ocp = $pref.ExaCore(
-                base_type; backend=backend, minimize=($p.criterion == :min) # not $(p.xxxx) as this info is known statically
+                base_type;
+                backend=backend,
+                minimize=($p.criterion == :min), # not $(p.xxxx) as this info is known statically
             )
             # _denull(e) = e # todo
             # _denull(e::$pref.Null{T}) where {T<:Real} = e.value # todo: Null must be imported by OC.jl
diff --git a/src/utils.jl b/src/utils.jl
index 5dee3dc..674f416 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -122,16 +122,19 @@ julia> subs2(subs2(e, :x0, :x, 0), :xf, :x, :N)
 :(x0 * (2 * x[3, N]))
 ```
 """
-function subs2(e, x, y, j; k = __symgen(:k))
-    foo(x, y, j) = (h, args...) -> begin
-        f = Expr(h, args...)
-        @match f begin
-            :($xx[$rg]) && if ((xx == x) && is_range(rg)) end => :([$y[$k, $j] for $k ∈ $rg])
-            :($xx[$i]) && if (xx == x) end => :($y[$i, $j])
-            _ => f
+function subs2(e, x, y, j; k=__symgen(:k))
+    foo(x, y, j) =
+        (h, args...) -> begin
+            f = Expr(h, args...)
+            @match f begin
+                :($xx[$rg]) && if ((xx == x) && is_range(rg))
+                end => :([$y[$k, $j] for $k in $rg])
+                :($xx[$i]) && if (xx == x)
+                end => :($y[$i, $j])
+                _ => f
+            end
         end
-    end
-    expr_it(e, foo(x, y, j), x -> x) 
+    expr_it(e, foo(x, y, j), x -> x)
 end
 
 """
@@ -195,13 +198,15 @@ julia> subs2m(e, :x0, :x, 0)
 :([((x[k, 0] + x[k, 0 + 1]) / 2) for k ∈ 1:3])
 ```
 """
-function subs2m(e, x, y, j; k = __symgen(:k))
+function subs2m(e, x, y, j; k=__symgen(:k))
     foo(x, y, j) =
         (h, args...) -> begin
             f = Expr(h, args...)
             @match f begin
-                :($xx[$rg]) && if ((xx == x) && is_range(rg)) end => :([($y[$k, $j] + $y[$k, $j + 1]) / 2 for $k ∈ $rg])
-                :($xx[$i]) && if (xx == x) end => :(($y[$i, $j] + $y[$i, $j + 1]) / 2)
+                :($xx[$rg]) && if ((xx == x) && is_range(rg))
+                end => :([($y[$k, $j] + $y[$k, $j + 1]) / 2 for $k in $rg])
+                :($xx[$i]) && if (xx == x)
+                end => :(($y[$i, $j] + $y[$i, $j + 1]) / 2)
                 _ => f
             end
         end
diff --git a/test/coverage.jl b/test/coverage.jl
index d3d2728..e34fc7a 100644
--- a/test/coverage.jl
+++ b/test/coverage.jl
@@ -14,4 +14,4 @@
 # ==============================================================================
 using Coverage
 using CTBase
-CTBase.postprocess_coverage(; root_dir=dirname(@__DIR__), dest_dir=".coverage")
\ No newline at end of file
+CTBase.postprocess_coverage(; root_dir=dirname(@__DIR__), dest_dir=".coverage")
diff --git a/test/runtests.jl b/test/runtests.jl
index c102dda..46f86d4 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -68,9 +68,7 @@ const SHOWTIMING = true
 CTBase.run_tests(;
     args=String.(ARGS),
     testset_name="CTParser tests",
-    available_tests=(
-        "suite/test_*",
-        ),
+    available_tests=("suite/test_*",),
     filename_builder=name -> "test_$(name).jl",
     funcname_builder=name -> "test_$(name)",
     verbose=VERBOSE,
@@ -89,6 +87,6 @@ if Base.JLOptions().code_coverage != 0
             julia --project -e 'using Pkg; Pkg.test("CTParser"; coverage=true); include("test/coverage.jl")'
 
         ================================================================================
-        """
+        """,
     )
-end
\ No newline at end of file
+end
diff --git a/test/suite/test_dynamics_exa.jl b/test/suite/test_dynamics_exa.jl
index 900900d..22096f1 100644
--- a/test/suite/test_dynamics_exa.jl
+++ b/test/suite/test_dynamics_exa.jl
@@ -7,7 +7,7 @@ activate_backend(:exa) # nota bene: needs to be executed before @def are expande
 
 function test_dynamics_exa()
     l_scheme = [:euler, :euler_implicit, :midpoint, :trapeze]
-    for scheme ∈ l_scheme
+    for scheme in l_scheme
         __test_dynamics_exa(; scheme=scheme)
         CUDA.functional() && __test_dynamics_exa(CUDABackend(); scheme=scheme)
     end
@@ -226,7 +226,8 @@ function __test_dynamics_exa(
             t ∈ [0, 1], time
             x ∈ R³, state
             u ∈ R², control
-            ∂(x)(t) == [x₁(t) + x₂(t) + x₃(t), 2x₁(t) * u₁(t) + x₂(t) * u₂(t), x₁(t) + 2x₂(t) - x₃(t)]
+            ∂(x)(t) ==
+            [x₁(t) + x₂(t) + x₃(t), 2x₁(t) * u₁(t) + x₂(t) * u₂(t), x₁(t) + 2x₂(t) - x₃(t)]
             x(0) == [1, 0, 0]
             x₁(1) => min
         end
@@ -340,14 +341,16 @@ function __test_dynamics_exa(
             u ∈ R³, control
             ∂(x)(t) == [
                 x₂(t),
-                u₁(t) * sin(x₇(t)) * sin(x₉(t)) + u₁(t) * cos(x₇(t)) * sin(x₈(t)) * cos(x₉(t)),
+                u₁(t) * sin(x₇(t)) * sin(x₉(t)) +
+                u₁(t) * cos(x₇(t)) * sin(x₈(t)) * cos(x₉(t)),
                 x₄(t),
-                u₁(t) * sin(x₇(t)) * cos(x₉(t)) - u₁(t) * cos(x₇(t)) * sin(x₈(t)) * sin(x₉(t)),
+                u₁(t) * sin(x₇(t)) * cos(x₉(t)) -
+                u₁(t) * cos(x₇(t)) * sin(x₈(t)) * sin(x₉(t)),
                 x₆(t),
                 u₁(t) * cos(x₇(t)) * cos(x₈(t)) - g,
                 u₂(t) * cos(x₇(t)) / cos(x₈(t)) + u₃(t) * sin(x₇(t)) / cos(x₈(t)),
                 -u₂(t) * sin(x₇(t)) + u₃(t) * cos(x₇(t)),
-                u₂(t) * cos(x₇(t)) * tan(x₈(t)) + u₃(t) * sin(x₇(t)) * tan(x₈(t))
+                u₂(t) * cos(x₇(t)) * tan(x₈(t)) + u₃(t) * sin(x₇(t)) * tan(x₈(t)),
             ]
             x(0) == [0, 0, 0, 0, 0, 0, 0, 0, 0]
             x₁(1) => min
@@ -460,4 +463,4 @@ function __test_dynamics_exa(
         m = discretise_exa(o; backend=backend, scheme=scheme)
         @test m isa ExaModels.ExaModel
     end
-end
\ No newline at end of file
+end
diff --git a/test/suite/test_exa_linalg.jl b/test/suite/test_exa_linalg.jl
index 3e2176a..d5044aa 100644
--- a/test/suite/test_exa_linalg.jl
+++ b/test/suite/test_exa_linalg.jl
@@ -508,7 +508,7 @@ function test_exa_linalg()
             @test v[1] isa ExaModels.AbstractNode
 
             # Create matrix from variable
-            A = [xvar[i, j] for (i, j) ∈ Base.product(1:2, 0:3)]
+            A = [xvar[i, j] for (i, j) in Base.product(1:2, 0:3)]
             @test size(A) == (2, 4)
             @test A isa Matrix
             @test A[1, 1] isa ExaModels.AbstractNode
@@ -541,7 +541,7 @@ function test_exa_linalg()
             @test result6[1] isa ExaModels.AbstractNode
 
             # Extract a 2×2 submatrix for square operations
-            A_square = [xvar[i, j] for (i, j) ∈ Base.product(1:2, 1:2)]
+            A_square = [xvar[i, j] for (i, j) in Base.product(1:2, 1:2)]
             @test size(A_square) == (2, 2)
 
             result7 = det(A_square)
@@ -748,7 +748,7 @@ function test_exa_linalg()
 
                 z3 = zeros(ExaModels.AbstractNode, 2, 2, 2)
                 @test size(z3) == (2, 2, 2)
-                @test z3 isa Array{<:ExaModels.AbstractNode, 3}
+                @test z3 isa Array{<:ExaModels.AbstractNode,3}
                 @test all(is_null_zero.(z3))
 
                 # Test ones with different dimensions
@@ -764,7 +764,7 @@ function test_exa_linalg()
 
                 o3 = ones(ExaModels.AbstractNode, 2, 2, 2)
                 @test size(o3) == (2, 2, 2)
-                @test o3 isa Array{<:ExaModels.AbstractNode, 3}
+                @test o3 isa Array{<:ExaModels.AbstractNode,3}
                 @test all(x -> x isa ExaModels.Null && isone(x.value), o3)
 
                 # Test with specific ExaModels types (Variable <: AbstractNode)
@@ -933,23 +933,23 @@ function test_exa_linalg()
 
                 # AbstractNode matrix + zeros
                 result1 = mat_nodes + [0 0; 0 0]
-                @test result1[1,1].value == x.value
-                @test result1[1,2].value == y.value
-                @test result1[2,1].value == z.value
-                @test result1[2,2].value == w.value
+                @test result1[1, 1].value == x.value
+                @test result1[1, 2].value == y.value
+                @test result1[2, 1].value == z.value
+                @test result1[2, 2].value == w.value
                 @test result1 isa Matrix{<:ExaModels.AbstractNode}
 
                 # Zeros + AbstractNode matrix
                 result2 = [0.0 0.0; 0.0 0.0] + mat_nodes
-                @test result2[1,1].value == x.value
-                @test result2[1,2].value == y.value
+                @test result2[1, 1].value == x.value
+                @test result2[1, 2].value == y.value
 
                 # Mixed: some zeros
                 result3 = mat_nodes + [0 1.0; 0 0]
-                @test result3[1,1].value == x.value  # x + 0 = x
-                @test result3[1,2].value == y.value + 1.0  # y + 1.0
-                @test result3[2,1].value == z.value  # z + 0 = z
-                @test result3[2,2].value == w.value  # w + 0 = w
+                @test result3[1, 1].value == x.value  # x + 0 = x
+                @test result3[1, 2].value == y.value + 1.0  # y + 1.0
+                @test result3[2, 1].value == z.value  # z + 0 = z
+                @test result3[2, 2].value == w.value  # w + 0 = w
             end
         end
 
@@ -985,18 +985,18 @@ function test_exa_linalg()
 
                 # AbstractNode matrix - zeros
                 result1 = mat_nodes - [0 0; 0 0]
-                @test result1[1,1].value == x.value
-                @test result1[1,2].value == y.value
-                @test result1[2,1].value == z.value
-                @test result1[2,2].value == w.value
+                @test result1[1, 1].value == x.value
+                @test result1[1, 2].value == y.value
+                @test result1[2, 1].value == z.value
+                @test result1[2, 2].value == w.value
                 @test result1 isa Matrix{<:ExaModels.AbstractNode}
 
                 # Mixed: some zeros
                 result2 = mat_nodes - [0 1.0; 0 0]
-                @test result2[1,1].value == x.value  # x - 0 = x
-                @test result2[1,2].value == y.value - 1.0  # y - 1.0
-                @test result2[2,1].value == z.value  # z - 0 = z
-                @test result2[2,2].value == w.value  # w - 0 = w
+                @test result2[1, 1].value == x.value  # x - 0 = x
+                @test result2[1, 2].value == y.value - 1.0  # y - 1.0
+                @test result2[2, 1].value == z.value  # z - 0 = z
+                @test result2[2, 2].value == w.value  # w - 0 = w
             end
         end
 
@@ -1091,17 +1091,32 @@ function test_exa_linalg()
 
             @testset "sum with zeros" begin
                 # All AbstractNode zeros: 0 + 0 + 0 = 0
-                result1 = sum([zero(ExaModels.AbstractNode), zero(ExaModels.AbstractNode), zero(ExaModels.AbstractNode)])
+                result1 = sum([
+                    zero(ExaModels.AbstractNode),
+                    zero(ExaModels.AbstractNode),
+                    zero(ExaModels.AbstractNode),
+                ])
                 @test result1 isa ExaModels.Null
                 @test iszero(result1.value)
 
                 # sum([zero(ExaModels.AbstractNode), zero(ExaModels.AbstractNode), x, zero(ExaModels.AbstractNode)]): 0 + 0 + x + 0 = x
-                result2 = sum([zero(ExaModels.AbstractNode), zero(ExaModels.AbstractNode), x, zero(ExaModels.AbstractNode)])
+                result2 = sum([
+                    zero(ExaModels.AbstractNode),
+                    zero(ExaModels.AbstractNode),
+                    x,
+                    zero(ExaModels.AbstractNode),
+                ])
                 @test result2 isa ExaModels.Null
                 @test result2.value == x.value
 
                 # sum with multiple non-zeros interspersed with zeros
-                result3 = sum([zero(ExaModels.AbstractNode), x, zero(ExaModels.AbstractNode), y, zero(ExaModels.AbstractNode)])
+                result3 = sum([
+                    zero(ExaModels.AbstractNode),
+                    x,
+                    zero(ExaModels.AbstractNode),
+                    y,
+                    zero(ExaModels.AbstractNode),
+                ])
                 @test result3 isa ExaModels.Null
                 @test result3.value == x.value + y.value
             end
@@ -1360,7 +1375,7 @@ function test_exa_linalg()
 
                 # Reshape Real vector to matrix
                 A_reshaped = reshape(v_num, 2, 2)
-                @test A_reshaped isa Union{Matrix, Base.ReshapedArray}
+                @test A_reshaped isa Union{Matrix,Base.ReshapedArray}
 
                 # Operations with reshaped array
                 result1 = x * A_reshaped
@@ -1369,7 +1384,7 @@ function test_exa_linalg()
 
                 # Reshape AbstractNode vector to matrix
                 mat_reshaped = reshape(vec_nodes, 2, 2)
-                @test mat_reshaped isa Union{Matrix, Base.ReshapedArray}
+                @test mat_reshaped isa Union{Matrix,Base.ReshapedArray}
 
                 result2 = 2.0 * mat_reshaped
                 @test size(result2) == (2, 2)
@@ -1388,7 +1403,7 @@ function test_exa_linalg()
 
                 # Reshape Real matrix to vector
                 v_reshaped = reshape(A_num, 4)
-                @test v_reshaped isa Union{Vector, Base.ReshapedArray}
+                @test v_reshaped isa Union{Vector,Base.ReshapedArray}
 
                 # Operations with reshaped array
                 result1 = x * v_reshaped
@@ -1397,7 +1412,7 @@ function test_exa_linalg()
 
                 # Reshape AbstractNode matrix to vector
                 vec_reshaped = reshape(mat_nodes, 4)
-                @test vec_reshaped isa Union{Vector, Base.ReshapedArray}
+                @test vec_reshaped isa Union{Vector,Base.ReshapedArray}
 
                 result2 = 2.0 * vec_reshaped
                 @test length(result2) == 4
diff --git a/test/suite/test_onepass_exa.jl b/test/suite/test_onepass_exa.jl
index 7e030c4..2236537 100644
--- a/test/suite/test_onepass_exa.jl
+++ b/test/suite/test_onepass_exa.jl
@@ -6,7 +6,7 @@ activate_backend(:exa) # nota bene: needs to be executed before @def are expande
 function test_onepass_exa()
     l_scheme = [:euler, :euler_implicit, :midpoint, :trapeze]
     #l_scheme = [:midpoint]
-    for scheme ∈ l_scheme
+    for scheme in l_scheme
         __test_onepass_exa(; scheme=scheme)
         CUDA.functional() && __test_onepass_exa(CUDABackend(); scheme=scheme)
     end
@@ -923,7 +923,7 @@ function __test_onepass_exa(
             v ≤ [1, 2, 3]
             v ≥ [1, 2, 3]
             v[1] ≤ 1
-            v[1] ≥ 1 
+            v[1] ≥ 1
             v[1:2] ≤ [1, 2]
             v[1:2] ≥ [1, 2]
             x[2](t) ≤ 1
@@ -1522,20 +1522,20 @@ function __test_onepass_exa(
         f₁(x, u) = 2x[1] * u[1] + x[2] * u[2]
         f₂(x) = x[1] + 2x[2] - x[3]
         f₃(x0, xf) = x0[2]^2 + sum(xf)^2
-        f₄(u) = sum(u.^2)
-        
+        f₄(u) = sum(u .^ 2)
+
         o1 = @def begin
             t ∈ [0, 1], time
             x ∈ R³, state
             u ∈ R², control
 
             x[1:2:3](0) == [1, 3]
- 
+
             ∂(x₁)(t) == sum(x(t))
             ∂(x₂)(t) == f₁(x(t), u(t))
-            ∂(x₃)(t) == f₂(x(t)) 
+            ∂(x₃)(t) == f₂(x(t))
 
-            f₃(x(0), x(1)) + 0.5∫( f₄(u(t)) ) → min
+            f₃(x(0), x(1)) + 0.5∫(f₄(u(t))) → min
         end
 
         N = 250
@@ -1551,12 +1551,12 @@ function __test_onepass_exa(
             u ∈ R², control
 
             x[1:2:3](0) == [1, 3]
- 
+
             ∂(x₁)(t) == x₁(t) + x₂(t) + x₃(t)
             ∂(x₂)(t) == 2x₁(t) * u₁(t) + x₂(t) * u₂(t)
             ∂(x₃)(t) == x₁(t) + 2x₂(t) - x₃(t)
 
-            (x₂(0)^2 + (x₁(1) + x₂(1) + x₃(1))^2) + 0.5∫( u₁(t)^2 + u₂(t)^2 ) → min
+            (x₂(0)^2 + (x₁(1) + x₂(1) + x₃(1))^2) + 0.5∫(u₁(t)^2 + u₂(t)^2) → min
         end
 
         m2, _ = discretise_exa_full(o2; grid_size=N, backend=backend, scheme=scheme)
@@ -1581,14 +1581,14 @@ function __test_onepass_exa(
             x ∈ R⁴, state
             u ∈ R², control
 
-            x(0) == [0, .1, .2, .3]
+            x(0) == [0, 0.1, 0.2, 0.3]
 
             ∂(x₁)(t) == g₁(x[1:2](t))
             ∂(x₂)(t) == g₂(u(t))
             ∂(x₃)(t) == sum(x[2:4](t))
             ∂(x₄)(t) == u₁(t)
 
-            sum(x[1:3](1))^2 + 0.5∫( sum(u(t).^2) ) → min
+            sum(x[1:3](1))^2 + 0.5∫(sum(u(t) .^ 2)) → min
         end
 
         N = 250
@@ -1605,14 +1605,14 @@ function __test_onepass_exa(
             x ∈ R⁴, state
             u ∈ R², control
 
-            x(0) == [0, .1, .2, .3]
+            x(0) == [0, 0.1, 0.2, 0.3]
 
             ∂(x₁)(t) == x₁(t)^2 + x₂(t)^2
             ∂(x₂)(t) == u₁(t) * u₂(t)
             ∂(x₃)(t) == x₂(t) + x₃(t) + x₄(t)
             ∂(x₄)(t) == u₁(t)
 
-            (x₁(1) + x₂(1) + x₃(1))^2 + 0.5∫( u₁(t)^2 + u₂(t)^2 ) → min
+            (x₁(1) + x₂(1) + x₃(1))^2 + 0.5∫(u₁(t)^2 + u₂(t)^2) → min
         end
 
         m2, _ = discretise_exa_full(o1; grid_size=N, backend=backend, scheme=scheme)
@@ -1637,7 +1637,7 @@ function __test_onepass_exa(
             x ∈ R³, state
             u ∈ R², control
 
-            sum(x(0).^2) == 1.5 
+            sum(x(0) .^ 2) == 1.5
             h₁(x(1)) ≤ 200
 
             ∂(x₁)(t) == u₁(t)
@@ -1646,7 +1646,7 @@ function __test_onepass_exa(
 
             h₂(u(t)) ≤ 10
 
-            sum(x(1))^2 + ∫( h₂(u(t)) ) → min
+            sum(x(1))^2 + ∫(h₂(u(t))) → min
         end
 
         N = 250
@@ -1671,7 +1671,7 @@ function __test_onepass_exa(
 
             u₁(t)^2 + u₂(t)^2 ≤ 10
 
-            (x₁(1) + x₂(1) + x₃(1))^2 + ∫( u₁(t)^2 + u₂(t)^2 ) → min
+            (x₁(1) + x₂(1) + x₃(1))^2 + ∫(u₁(t)^2 + u₂(t)^2) → min
         end
 
         m2, _ = discretise_exa_full(o2; grid_size=N, backend=backend, scheme=scheme)
@@ -1681,71 +1681,74 @@ function __test_onepass_exa(
 
         __atol = 1e-9
         @test obj1 - obj2 ≈ 0 atol = __atol
-    end 
+    end
 
     # todo: test below inactived on GPU because run is unstable
-    if isnothing(backend) test_name = "use case no. 7: mixed vectorisation ($backend_name, $scheme)"
-    @testset "$test_name" begin
-        println(test_name)
-
-        # User-defined functions
-        p₁(x, u) = x[1] * u[1] + x[2] * u[2]
-        p₂(x) = x[1]^2 + x[2]^2 + x[3]^2
-
-        # Vectorised version with mixed patterns
-        o1 = @def begin
-            t ∈ [0, 1], time
-            x ∈ R³, state
-            u ∈ R², control
-
-            x[1:2](0) == [0, 0.1]
-            -0.1 ≤ x₃(0) ≤ 0.1
-            sum(x(1)) == 0.2
-
-            ∂(x₁)(t) == p₁(x[1:2](t), u(t))
-            ∂(x₂)(t) == sum(u(t))
-            ∂(x₃)(t) == x₁(t)
-
-            p₂(x(t)) ≤ 50
-
-            (p₂(x(0)) + sum(x[1:2](1))^2) + 0.5∫( sum(u(t).^2) ) → min
+    if isnothing(backend)
+        test_name = "use case no. 7: mixed vectorisation ($backend_name, $scheme)"
+        @testset "$test_name" begin
+            println(test_name)
+
+            # User-defined functions
+            p₁(x, u) = x[1] * u[1] + x[2] * u[2]
+            p₂(x) = x[1]^2 + x[2]^2 + x[3]^2
+
+            # Vectorised version with mixed patterns
+            o1 = @def begin
+                t ∈ [0, 1], time
+                x ∈ R³, state
+                u ∈ R², control
+
+                x[1:2](0) == [0, 0.1]
+                -0.1 ≤ x₃(0) ≤ 0.1
+                sum(x(1)) == 0.2
+
+                ∂(x₁)(t) == p₁(x[1:2](t), u(t))
+                ∂(x₂)(t) == sum(u(t))
+                ∂(x₃)(t) == x₁(t)
+
+                p₂(x(t)) ≤ 50
+
+                (p₂(x(0)) + sum(x[1:2](1))^2) + 0.5∫(sum(u(t) .^ 2)) → min
+            end
+
+            N = 250
+            max_iter = 10
+            m1, _ = discretise_exa_full(o1; grid_size=N, backend=backend, scheme=scheme)
+            @test m1 isa ExaModels.ExaModel
+            sol1 = madnlp(m1; tol=tolerance, max_iter=max_iter, kwargs...)
+            obj1 = sol1.objective
+
+            # Non-vectorised version
+            o2 = @def begin
+                t ∈ [0, 1], time
+                x ∈ R³, state
+                u ∈ R², control
+
+                x₁(0) == 0
+                x₂(0) == 0.1
+                -0.1 ≤ x₃(0) ≤ 0.1
+                x₁(1) + x₂(1) + x₃(1) == 0.2
+
+                ∂(x₁)(t) == x₁(t) * u₁(t) + x₂(t) * u₂(t)
+                ∂(x₂)(t) == u₁(t) + u₂(t)
+                ∂(x₃)(t) == x₁(t)
+
+                x₁(t)^2 + x₂(t)^2 + x₃(t)^2 ≤ 50
+
+                (x₁(0)^2 + x₂(0)^2 + x₃(0)^2 + (x₁(1) + x₂(1))^2) +
+                0.5∫(u₁(t)^2 + u₂(t)^2) → min
+            end
+
+            m2, _ = discretise_exa_full(o2; grid_size=N, backend=backend, scheme=scheme)
+            @test m2 isa ExaModels.ExaModel
+            sol2 = madnlp(m2; tol=tolerance, max_iter=max_iter, kwargs...)
+            obj2 = sol2.objective
+
+            __atol = 1e-9
+            @test obj1 - obj2 ≈ 0 atol = __atol
         end
-
-        N = 250
-        max_iter = 10
-        m1, _ = discretise_exa_full(o1; grid_size=N, backend=backend, scheme=scheme)
-        @test m1 isa ExaModels.ExaModel
-        sol1 = madnlp(m1; tol=tolerance, max_iter=max_iter, kwargs...)
-        obj1 = sol1.objective
-
-        # Non-vectorised version
-        o2 = @def begin
-            t ∈ [0, 1], time
-            x ∈ R³, state
-            u ∈ R², control
-
-            x₁(0) == 0
-            x₂(0) == 0.1
-            -0.1 ≤ x₃(0) ≤ 0.1
-            x₁(1) + x₂(1) + x₃(1) == 0.2 
-
-            ∂(x₁)(t) == x₁(t) * u₁(t) + x₂(t) * u₂(t)
-            ∂(x₂)(t) == u₁(t) + u₂(t)
-            ∂(x₃)(t) == x₁(t)
-
-            x₁(t)^2 + x₂(t)^2 + x₃(t)^2 ≤ 50
-
-            (x₁(0)^2 + x₂(0)^2 + x₃(0)^2 + (x₁(1) + x₂(1))^2) + 0.5∫( u₁(t)^2 + u₂(t)^2 ) → min
-        end
-
-        m2, _ = discretise_exa_full(o2; grid_size=N, backend=backend, scheme=scheme)
-        @test m2 isa ExaModels.ExaModel
-        sol2 = madnlp(m2; tol=tolerance, max_iter=max_iter, kwargs...)
-        obj2 = sol2.objective
-
-        __atol = 1e-9
-        @test obj1 - obj2 ≈ 0 atol = __atol
-    end end
+    end
 
     test_name = "use case no. 8: vectorised dynamics ($backend_name, $scheme)"
     @testset "$test_name" begin
@@ -1770,7 +1773,7 @@ function __test_onepass_exa(
             #∂(x₂)(t) == -x₁(t) + u(t) 
             ∂(x₂)(t) == A[2, :]' * x(t) + u(t) * b[2]
             #0.5∫( x₁(t)^2 + x₂(t)^2 + u(t)^2 ) → min
-            0.5∫( x(t)' * Q * x(t) + u(t)' * R * u(t) ) → min
+            0.5∫(x(t)' * Q * x(t) + u(t)' * R * u(t)) → min
         end
 
         N = 250
@@ -1787,10 +1790,10 @@ function __test_onepass_exa(
             u ∈ R, control
             x(0) == x0
             ∂(x₁)(t) == x₂(t)
-            ∂(x₂)(t) == -x₁(t) + u(t) 
-            0.5∫( x₁(t)^2 + x₂(t)^2 + u(t)^2 ) → min
+            ∂(x₂)(t) == -x₁(t) + u(t)
+            0.5∫(x₁(t)^2 + x₂(t)^2 + u(t)^2) → min
         end
-        
+
         m2, _ = discretise_exa_full(o2; grid_size=N, backend=backend, scheme=scheme)
         @test m2 isa ExaModels.ExaModel
         sol2 = madnlp(m2; tol=tolerance, max_iter=max_iter, kwargs...)
@@ -1798,7 +1801,7 @@ function __test_onepass_exa(
 
         __atol = 1e-9
         @test obj1 - obj2 ≈ 0 atol = __atol
-    end 
+    end
 
     test_name = "use case no. 9: vectorised dynamics ($backend_name, $scheme)"
     @testset "$test_name" begin
@@ -1813,9 +1816,9 @@ function __test_onepass_exa(
             t ∈ [0, 1], time
             x ∈ R², state
             u ∈ R, control
-            c(x(0), x(1)) == [-1, 0, 0, 0] 
-            ẋ(t) == A * x(t) + u(t) * b 
-            0.5∫( u(t)^2 ) → min
+            c(x(0), x(1)) == [-1, 0, 0, 0]
+            ẋ(t) == A * x(t) + u(t) * b
+            0.5∫(u(t)^2) → min
         end
 
         N = 250
@@ -1830,13 +1833,13 @@ function __test_onepass_exa(
             t ∈ [0, 1], time
             x ∈ R², state
             u ∈ R, control
-            x(0) == [-1, 0] 
-            x(1) == [0, 0] 
+            x(0) == [-1, 0]
+            x(1) == [0, 0]
             ∂(x₁)(t) == x₂(t)
-            ∂(x₂)(t) == u(t) 
-            0.5∫( u(t)^2 ) → min
+            ∂(x₂)(t) == u(t)
+            0.5∫(u(t)^2) → min
         end
-        
+
         m2, _ = discretise_exa_full(o2; grid_size=N, backend=backend, scheme=scheme)
         @test m2 isa ExaModels.ExaModel
         sol2 = madnlp(m2; tol=tolerance, max_iter=max_iter, kwargs...)
@@ -1844,53 +1847,51 @@ function __test_onepass_exa(
 
         __atol = 1e-9
         @test obj1 - obj2 ≈ 0 atol = __atol
-    end 
+    end
 
     test_name = "use case no. 10 (Goddard): vectorised dynamics ($backend_name, $scheme)"
     @testset "$test_name" begin
         println(test_name)
 
-        t0 = 0      
-        r0 = 1      
-        v0 = 0      
-        m0 = 1      
-        vmax = 0.1  
-        mf = 0.6    
+        t0 = 0
+        r0 = 1
+        v0 = 0
+        m0 = 1
+        vmax = 0.1
+        mf = 0.6
         Cd = 310
         Tmax = 3.5
         β = 500
         b = 2
-        
+
         o1 = @def begin
-        
             tf ∈ R, variable
             t ∈ [t0, tf], time
             x = (r, v, m) ∈ R³, state
             u ∈ R, control
-        
+
             x(t0) == [r0, v0, m0]
             m(tf) == mf
             0 ≤ u(t) ≤ 1
             r(t) ≥ r0
             0 ≤ v(t) ≤ vmax
-        
+
             ẋ(t) == F0(x(t)) + u(t) * F1(x(t))
-        
+
             -r(tf) → min
-        
         end
-        
+
         F0(x) = begin
             r, v, m = x
             D = Cd * v^2 * exp(-β*(r - 1))
             return [v, -D/m - 1/r^2, 0]
         end
-        
+
         F1(x) = begin
             r, v, m = x
             return [0, Tmax/m, -b*Tmax]
         end
-        
+
         N = 250
         max_iter = 100
         m1, _ = discretise_exa_full(o1; grid_size=N, backend=backend, scheme=scheme)
@@ -1901,5 +1902,4 @@ function __test_onepass_exa(
         obj2 = -1.0125766794048943e+00 # from midpoint with N = 250 
         @test obj1 - obj2 ≈ 0 atol = __atol
     end
-
-end
\ No newline at end of file
+end
diff --git a/test/suite/test_onepass_fun_bis.jl b/test/suite/test_onepass_fun_bis.jl
index d9a7e22..7f9a90b 100644
--- a/test/suite/test_onepass_fun_bis.jl
+++ b/test/suite/test_onepass_fun_bis.jl
@@ -552,9 +552,7 @@ function test_onepass_fun_bis()
         # Test that p_constraint! correctly identifies invalid constraints
         # Mixed initial state and control: x1(0) * u1(t) + u2(t)^2 <= 1
         # This should result in constraint_type returning :other
-        ex = CTParser.p_constraint!(
-            p, p_ocp, nothing, :(x[1](0) * u[1](t) + u[2](t)^2), 1
-        )
+        ex = CTParser.p_constraint!(p, p_ocp, nothing, :(x[1](0) * u[1](t) + u[2](t)^2), 1)
         @test ex isa Expr
         @test_throws ParsingError eval(ex)
 
diff --git a/test/suite/test_utils.jl b/test/suite/test_utils.jl
index 580a11e..9fa382b 100644
--- a/test/suite/test_utils.jl
+++ b/test/suite/test_utils.jl
@@ -54,77 +54,77 @@ function test_utils()
         @testset "range indexing (new)" begin
             # Test 5: Simple range 1:3
             e = :(x0[1:3])
-            result = subs2(e, :x0, :x, 0; k = :k)
-            @test result == :([x[k, 0] for k ∈ 1:3])
+            result = subs2(e, :x0, :x, 0; k=:k)
+            @test result == :([x[k, 0] for k in 1:3])
 
             # Test 6: Range with step 1:2:5
             e = :(x0[1:2:5])
-            result = subs2(e, :x0, :x, 0; k = :k)
-            @test result == :([x[k, 0] for k ∈ 1:2:5])
+            result = subs2(e, :x0, :x, 0; k=:k)
+            @test result == :([x[k, 0] for k in 1:2:5])
 
             # Test 7: Range with symbolic bounds
             e = :(x0[1:n])
-            result = subs2(e, :x0, :x, 0; k = :k)
-            @test result == :([x[k, 0] for k ∈ 1:n])
+            result = subs2(e, :x0, :x, 0; k=:k)
+            @test result == :([x[k, 0] for k in 1:n])
 
             # Test 8: Multiple ranges in same expression
             e = :(x0[1:3] + xf[2:4])
-            result = subs2(subs2(e, :x0, :x, 0; k = :k1), :xf, :x, :N; k = :k2)
-            @test result == :([x[k1, 0] for k1 ∈ 1:3] + [x[k2, N] for k2 ∈ 2:4])
+            result = subs2(subs2(e, :x0, :x, 0; k=:k1), :xf, :x, :N; k=:k2)
+            @test result == :([x[k1, 0] for k1 in 1:3] + [x[k2, N] for k2 in 2:4])
 
             # Test 9: Range inside function call
             e = :(sum(x0[1:n]))
-            result = subs2(e, :x0, :x, 0; k = :k)
-            @test result == :(sum([x[k, 0] for k ∈ 1:n]))
+            result = subs2(e, :x0, :x, 0; k=:k)
+            @test result == :(sum([x[k, 0] for k in 1:n]))
         end
 
         @testset "mixed scalar and range" begin
             # Test 10: Expression with both scalars and ranges
             e = :(x0[1] + x0[2:4] + x0[5])
-            result = subs2(e, :x0, :x, 0; k = :k)
+            result = subs2(e, :x0, :x, 0; k=:k)
             # x0[1] → x[1, 0]
             # x0[2:4] → [x[k, 0] for k ∈ 2:4]
             # x0[5] → x[5, 0]
-            @test result == :(x[1, 0] + [x[k, 0] for k ∈ 2:4] + x[5, 0])
+            @test result == :(x[1, 0] + [x[k, 0] for k in 2:4] + x[5, 0])
         end
 
         @testset "nested and complex expressions" begin
             # Test 11: Nested function calls with ranges
             e = :(norm(x0[1:3]) + cos(x0[4]))
-            result = subs2(e, :x0, :x, 0; k = :k)
-            @test result == :(norm([x[k, 0] for k ∈ 1:3]) + cos(x[4, 0]))
+            result = subs2(e, :x0, :x, 0; k=:k)
+            @test result == :(norm([x[k, 0] for k in 1:3]) + cos(x[4, 0]))
 
             # Test 12: Range in matrix operations
             e = :(A * x0[1:n])
-            result = subs2(e, :x0, :x, 0; k = :k)
-            @test result == :(A * [x[k, 0] for k ∈ 1:n])
+            result = subs2(e, :x0, :x, 0; k=:k)
+            @test result == :(A * [x[k, 0] for k in 1:n])
 
             # Test 13: Multiple substitutions with symbolic j
             e = :(x0[1:3] + xf[2:4])
-            result = subs2(subs2(e, :x0, :x, :j; k = :k1), :xf, :x, :(j+1); k = :k2)
-            @test result == :([x[k1, j] for k1 ∈ 1:3] + [x[k2, j+1] for k2 ∈ 2:4])
+            result = subs2(subs2(e, :x0, :x, :j; k=:k1), :xf, :x, :(j+1); k=:k2)
+            @test result == :([x[k1, j] for k1 in 1:3] + [x[k2, j + 1] for k2 in 2:4])
         end
 
         @testset "edge cases" begin
             # Test 14: Single-element range (should still create comprehension)
             e = :(x0[1:1])
-            result = subs2(e, :x0, :x, 0; k = :k)
-            @test result == :([x[k, 0] for k ∈ 1:1])
+            result = subs2(e, :x0, :x, 0; k=:k)
+            @test result == :([x[k, 0] for k in 1:1])
 
             # Test 15: Wrong variable name (should not substitute)
             e = :(y0[1:3])
-            result = subs2(e, :x0, :x, 0; k = :k)
+            result = subs2(e, :x0, :x, 0; k=:k)
             @test result == e  # Unchanged
 
             # Test 16: Complex symbolic j expression
             e = :(x0[1:3])
-            result = subs2(e, :x0, :x, :grid_size; k = :k)
-            @test result == :([x[k, grid_size] for k ∈ 1:3])
+            result = subs2(e, :x0, :x, :grid_size; k=:k)
+            @test result == :([x[k, grid_size] for k in 1:3])
 
             # Test 17: Scalar index that is a range expression (should not match)
             # This tests that we properly distinguish i (scalar) from 1:3 (range)
             e = :(x0[i])
-            result = subs2(e, :x0, :x, 0; k = :k)
+            result = subs2(e, :x0, :x, 0; k=:k)
             @test result == :(x[i, 0])  # Scalar behavior
         end
 
@@ -155,36 +155,35 @@ function test_utils()
         @testset "range indexing" begin
             # Test 1: Basic range substitution
             e = :(x0[1:3])
-            result = subs2m(e, :x0, :x, 0; k = :k)
-            @test result == :([((x[k, 0] + x[k, 0 + 1]) / 2) for k ∈ 1:3])
+            result = subs2m(e, :x0, :x, 0; k=:k)
+            @test result == :([((x[k, 0] + x[k, 0 + 1]) / 2) for k in 1:3])
 
             # Test 2: Range with step
             e = :(x0[1:2:5])
-            result = subs2m(e, :x0, :x, 0; k = :k)
-            @test result == :([((x[k, 0] + x[k, 0 + 1]) / 2) for k ∈ 1:2:5])
+            result = subs2m(e, :x0, :x, 0; k=:k)
+            @test result == :([((x[k, 0] + x[k, 0 + 1]) / 2) for k in 1:2:5])
 
             # Test 3: Range in arithmetic expression
             e = :(2 * x0[1:3])
-            result = subs2m(e, :x0, :x, 0; k = :k)
-            @test result == :(2 * [((x[k, 0] + x[k, 0 + 1]) / 2) for k ∈ 1:3])
+            result = subs2m(e, :x0, :x, 0; k=:k)
+            @test result == :(2 * [((x[k, 0] + x[k, 0 + 1]) / 2) for k in 1:3])
 
             # Test 4: Multiple ranges in same expression
             e = :(x0[1:2] + xf[2:4])
-            result = subs2m(subs2m(e, :x0, :x, 0; k = :k), :xf, :x, :N; k = :k)
+            result = subs2m(subs2m(e, :x0, :x, 0; k=:k), :xf, :x, :N; k=:k)
             @test result == :(
-                [((x[k, 0] + x[k, 0 + 1]) / 2) for k ∈ 1:2] +
-                [((x[k, N] + x[k, N + 1]) / 2) for k ∈ 2:4]
+                [((x[k, 0] + x[k, 0 + 1]) / 2) for k in 1:2] + [((x[k, N] + x[k, N + 1]) / 2) for k in 2:4]
             )
 
             # Test 5: Range with symbolic j
             e = :(x0[1:3])
-            result = subs2m(e, :x0, :x, :j; k = :k)
-            @test result == :([((x[k, j] + x[k, j + 1]) / 2) for k ∈ 1:3])
+            result = subs2m(e, :x0, :x, :j; k=:k)
+            @test result == :([((x[k, j] + x[k, j + 1]) / 2) for k in 1:3])
 
             # Test 6: Single-element range
             e = :(x0[2:2])
-            result = subs2m(e, :x0, :x, 0; k = :k)
-            @test result == :([((x[k, 0] + x[k, 0 + 1]) / 2) for k ∈ 2:2])
+            result = subs2m(e, :x0, :x, 0; k=:k)
+            @test result == :([((x[k, 0] + x[k, 0 + 1]) / 2) for k in 2:2])
         end
 
         @testset "backward compatibility" begin
diff --git a/test/utils.jl b/test/utils.jl
index 619652d..0e693f9 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -33,4 +33,4 @@ function discretise_exa_full(
     return build_exa(;
         scheme=scheme, grid_size=grid_size, backend=backend, init=init, base_type=base_type
     )
-end
\ No newline at end of file
+end