diff --git a/src/Expressions.jl b/src/Expressions.jl
index 17a769d..da8dddf 100644
--- a/src/Expressions.jl
+++ b/src/Expressions.jl
@@ -370,66 +370,251 @@ $(TYPEDEF)
 """
 abstract type AbstractExpressionFunction{T, N, D} <: Function end
 
+########################################################
+# Flat expression-tree representation
+#
+# `ScalarExpressionFunction` stores its expression tree as a fixed-size
+# `NTuple` of `FlatNode` records — a struct of plain integer fields and a
+# value of type `T`.  The whole thing is `isbits`, which means:
+#
+#   * It passes through KernelAbstractions kernels as a value argument and
+#     can be called on the GPU directly — no host↔device round-trip per
+#     time step for BC evaluation.
+#   * It survives `juliac --trim`: no closures, no `@eval`, no
+#     `RuntimeGeneratedFunctions`, just data.
+#
+# `FEC_EXPR_MAX_NODES` is the maximum tree size.  Carina's largest
+# inlined expression today is ~25 nodes, but second-derivative trees from
+# the symbolic differentiator can grow to ~100 nodes for Gaussian-pulse
+# BCs; 256 gives 2-3× headroom at ~6 KB per function — negligible memory
+# for typical BC/IC counts.  Trees that exceed the cap raise at
+# construction time.
+########################################################
+
+const FEC_EXPR_MAX_NODES = 256
+
 """
+$(TYPEDEF)
 $(TYPEDFIELDS)
+
+One node of a flattened expression tree.  Children are referenced by
+1-based index into the parent array; index 0 marks "no child".  Leaves
+carry either a constant `val` or a `feature` (1-based variable index).
+"""
+struct FlatNode{T <: Number}
+    degree::UInt8     # 0 = leaf, 1 = unary, 2 = binary
+    op::UInt8         # FUNC_* or BINARY_* op code; 0 for leaves
+    constant::Bool    # leaf form: true ⇒ use `val`, false ⇒ use `feature`
+    val::T            # leaf value when `constant`
+    feature::UInt16   # leaf variable index (1-based) when !`constant`
+    l_idx::UInt16     # index of left  child (0 if leaf)
+    r_idx::UInt16     # index of right child (0 if leaf or unary)
+end
+
+@inline FlatNode{T}() where T <: Number =
+    FlatNode{T}(UInt8(0), UInt8(0), false, zero(T), UInt16(0), UInt16(0), UInt16(0))
+
+# Flatten a recursive `Node{T, D}` (produced by FEC's Pratt parser) into a
+# preorder NTuple of `FlatNode{T}` records.  The first entry is the root;
+# children follow in DFS order.
+function _flatten_visit!(buf::Vector{FlatNode{T}}, node::Node{T, D})::UInt16 where {T, D}
+    if node.degree == 0
+        if node.constant
+            push!(buf, FlatNode{T}(UInt8(0), UInt8(0), true,
+                                   T(node.val), UInt16(0),
+                                   UInt16(0), UInt16(0)))
+        else
+            push!(buf, FlatNode{T}(UInt8(0), UInt8(0), false,
+                                   zero(T), UInt16(node.feature),
+                                   UInt16(0), UInt16(0)))
+        end
+        return UInt16(length(buf))
+    elseif node.degree == 1
+        push!(buf, FlatNode{T}())                   # reserve own slot
+        my_idx = length(buf)
+        l = _flatten_visit!(buf, node.l)
+        buf[my_idx] = FlatNode{T}(UInt8(1), UInt8(node.op), false,
+                                  zero(T), UInt16(0), l, UInt16(0))
+        return UInt16(my_idx)
+    else                                            # degree == 2
+        push!(buf, FlatNode{T}())                   # reserve own slot
+        my_idx = length(buf)
+        l = _flatten_visit!(buf, node.l)
+        r = _flatten_visit!(buf, node.r)
+        buf[my_idx] = FlatNode{T}(UInt8(2), UInt8(node.op), false,
+                                  zero(T), UInt16(0), l, r)
+        return UInt16(my_idx)
+    end
+end
+
+function _flatten(root::Node{T, D}) where {T, D}
+    buf = FlatNode{T}[]
+    sizehint!(buf, FEC_EXPR_MAX_NODES)
+    _flatten_visit!(buf, root)
+    n_active = length(buf)
+    n_active <= FEC_EXPR_MAX_NODES || error(
+        "expression too large: $n_active nodes (max $FEC_EXPR_MAX_NODES)"
+    )
+    # Pad to fixed length with default nodes so the resulting NTuple type
+    # has a constant size at the type level.
+    while length(buf) < FEC_EXPR_MAX_NODES
+        push!(buf, FlatNode{T}())
+    end
+    nodes = NTuple{FEC_EXPR_MAX_NODES, FlatNode{T}}(buf)
+    return nodes, UInt16(n_active)
+end
+
+# Inverse used by `differentiate`: rebuild a recursive `Node{T, D}` from a
+# flat NTuple so the existing recursive symbolic differentiator can run on
+# it.  Called once per `differentiate` invocation, off the GPU.
+function _unflatten(nodes::NTuple{N, FlatNode{T}}, idx::Integer = 1) where {N, T}
+    n = nodes[idx]
+    if n.degree == 0
+        if n.constant
+            return Node{T, DEFAULT_MAX_DEGREE}(; val = n.val)
+        else
+            return Node{T, DEFAULT_MAX_DEGREE}(; feature = Int(n.feature))
+        end
+    elseif n.degree == 1
+        l = _unflatten(nodes, n.l_idx)
+        return Node{T, DEFAULT_MAX_DEGREE}(; op = Int(n.op), l = l)
+    else
+        l = _unflatten(nodes, n.l_idx)
+        r = _unflatten(nodes, n.r_idx)
+        return Node{T, DEFAULT_MAX_DEGREE}(; op = Int(n.op), l = l, r = r)
+    end
+end
+
+# Op-code dispatch — open-coded if/elseif so the compiler can inline
+# branches inside KA kernels without resorting to a function table.
+@inline function _apply_unary_op(::Type{T}, op::UInt8, u::T) where T <: Number
+    if     op == UInt8(FUNC_MINUS); return -u
+    elseif op == UInt8(FUNC_COS);   return cos(u)
+    elseif op == UInt8(FUNC_COSH);  return cosh(u)
+    elseif op == UInt8(FUNC_EXP);   return exp(u)
+    elseif op == UInt8(FUNC_LOG);   return log(u)
+    elseif op == UInt8(FUNC_SIN);   return sin(u)
+    elseif op == UInt8(FUNC_SINH);  return sinh(u)
+    elseif op == UInt8(FUNC_SQRT);  return sqrt(u)
+    elseif op == UInt8(FUNC_TAN);   return tan(u)
+    elseif op == UInt8(FUNC_TANH);  return tanh(u)
+    end
+    return T(NaN)
+end
+
+@inline function _apply_binary_op(::Type{T}, op::UInt8, u::T, v::T) where T <: Number
+    if     op == UInt8(BINARY_PLUS);     return u + v
+    elseif op == UInt8(BINARY_MINUS);    return u - v
+    elseif op == UInt8(BINARY_MULTIPLY); return u * v
+    elseif op == UInt8(BINARY_DIVIDE);   return u / v
+    elseif op == UInt8(BINARY_POWER);    return u ^ v
+    end
+    return T(NaN)
+end
+
+# Recursive evaluator over the flat NTuple.  Depth is bounded by the
+# expression tree height (≤ ~10 for the expressions Carina uses today),
+# so GPUCompiler handles the recursion without stack pressure.
+function _eval_node(nodes::NTuple{N, FlatNode{T}}, idx::UInt16,
+                    vars) where {N, T}
+    n = nodes[idx]
+    if n.degree == 0
+        if n.constant
+            return n.val
+        else
+            return T(vars[n.feature])
+        end
+    elseif n.degree == 1
+        u = _eval_node(nodes, n.l_idx, vars)
+        return _apply_unary_op(T, n.op, u)
+    else
+        u = _eval_node(nodes, n.l_idx, vars)
+        v = _eval_node(nodes, n.r_idx, vars)
+        return _apply_binary_op(T, n.op, u, v)
+    end
+end
+
+"""
 $(TYPEDEF)
+$(TYPEDFIELDS)
+
+Scalar expression function as a flat, `isbits` value — usable as a
+KernelAbstractions kernel argument and trim-mode safe under `juliac`.
+The trailing variable is conventionally time; FEC's juliac-safe
+`DirichletBCs` constructor uses `num_vars` as the time-derivative index.
 """
-struct ScalarExpressionFunction{T <: Number} <: AbstractExpressionFunction{T, Node{T, DEFAULT_MAX_DEGREE}, ntuple_type}
-    expr::Expression{T, Node{T, DEFAULT_MAX_DEGREE}, ntuple_type}
-    num_vars::Int
+struct ScalarExpressionFunction{T <: Number} <: AbstractExpressionFunction{T, FlatNode{T}, ntuple_type}
+    nodes::NTuple{FEC_EXPR_MAX_NODES, FlatNode{T}}
+    n_active::UInt16
+    num_vars::UInt8
 
     """
     $(TYPEDSIGNATURES)
+
+    Parse `string` as an expression in the variable namespace `var_names`
+    and store the resulting tree in flat form.  `var_names` is consumed by
+    the parser to bind identifiers to feature indices; it is not retained
+    on the resulting function.
     """
     function ScalarExpressionFunction{T}(string::String, var_names::Vector{String}) where T <: Number
         p = Parser{T}(string, var_names)
-        # params = _find_parameters(p)
         _reset!(p)
         ast = _parse_statement(p, 0)
-        expr = Expression(ast; operators, var_names)
-        new{T}(expr, length(var_names))
+        nodes, n_active = _flatten(ast)
+        new{T}(nodes, n_active, UInt8(length(var_names)))
     end
 
     """
     $(TYPEDSIGNATURES)
 
-    Build a `ScalarExpressionFunction` directly from a prebuilt
-    `DynamicExpressions.Expression` — used by [`differentiate`](@ref) to wrap
-    the result of a tree rewrite without round-tripping through the parser.
+    Build a `ScalarExpressionFunction` from a prebuilt flat NTuple — used
+    internally by [`differentiate`](@ref) to wrap the result of a tree
+    rewrite without round-tripping through the parser.
     """
     function ScalarExpressionFunction{T}(
-        expr::Expression{T, Node{T, DEFAULT_MAX_DEGREE}, ntuple_type},
-        num_vars::Int
+        nodes::NTuple{FEC_EXPR_MAX_NODES, FlatNode{T}},
+        n_active::UInt16,
+        num_vars::Integer
     ) where T <: Number
-        new{T}(expr, num_vars)
+        new{T}(nodes, n_active, UInt8(num_vars))
     end
 end
 
 Base.eltype(::ScalarExpressionFunction{T}) where T <: Number = T
 
-function (f::ScalarExpressionFunction)(var::T) where T <: Number
+# Single scalar
+function (f::ScalarExpressionFunction{T})(var::T) where T <: Number
     @assert f.num_vars == 1
-    return f.expr(SMatrix{1, 1, T, 1}(var))[1]
+    return _eval_node(f.nodes, UInt16(1), SVector{1, T}(var))
 end
 
-# fall back function, not that efficient though
-function (f::ScalarExpressionFunction)(vars::AbstractVector{T}) where T <: Number
-    vars = reshape(vars, length(vars), 1)
-    return f.expr(vars)[1]
+# Vector of variable values (no `t` overload)
+function (f::ScalarExpressionFunction{T})(vars::AbstractVector{T}) where T <: Number
+    @assert length(vars) == Int(f.num_vars) "expected $(Int(f.num_vars)) variables, got $(length(vars))"
+    return _eval_node(f.nodes, UInt16(1), vars)
 end
 
-# for ic type funcs
-function (f::ScalarExpressionFunction)(X::SVector{ND, T}) where {ND, T <: Number}
-    @assert f.num_vars == ND "You need $ND variables for this function"
-    X = SMatrix{ND, 1, T, ND}(X.data)
-    return f.expr(X)[1]
+# IC-style call: ND spatial coords (no time variable in the expression).
+function (f::ScalarExpressionFunction{T})(X::SVector{ND, T}) where {ND, T <: Number}
+    @assert Int(f.num_vars) == ND "You need $ND variables for this function"
+    return _eval_node(f.nodes, UInt16(1), X)
 end
 
-# for bc type funcs
-function (f::ScalarExpressionFunction)(X::SVector{ND, T}, t::T) where {ND, T <: Number}
-    @assert f.num_vars == ND + 1 "You need $(ND + 1) variables for this function"
-    vars = SMatrix{ND + 1, 1, T, ND + 1}(X..., t)
-    return f.expr(vars)[1]
+# BC-style call: ND spatial coords + scalar time, packed into a stack-
+# allocated SVector so the call survives KA kernels.
+function (f::ScalarExpressionFunction{T})(X::SVector{ND, T}, t::T) where {ND, T <: Number}
+    @assert Int(f.num_vars) == ND + 1 "You need $(ND + 1) variables for this function"
+    if ND == 1
+        vars = SVector{2, T}(X[1], t)
+    elseif ND == 2
+        vars = SVector{3, T}(X[1], X[2], t)
+    elseif ND == 3
+        vars = SVector{4, T}(X[1], X[2], X[3], t)
+    else
+        # Generic path (very unlikely; covered for completeness).  Allocates.
+        vars = T[X...; t]
+    end
+    return _eval_node(f.nodes, UInt16(1), vars)
 end
 
 """
@@ -471,12 +656,13 @@ function (f::VectorExpressionFunction)(X::SVector{ND, T}, t::T) where {ND, T <:
 end
 
 ########################################################
-# Symbolic differentiation on Node{T, D} trees.
+# Symbolic differentiation on the recursive Node form.
 #
-# The grammar is finite (10 unary + 5 binary operators), so the chain rule
-# can be applied by tree rewriting in ~80 lines of pure Julia.  This avoids
-# pulling in ForwardDiff/Zygote/Symbolics, so the result survives
-# `juliac --trim`.  All helpers operate on values; no closures are formed.
+# The differentiator is a pure tree-rewrite over FEC's closed grammar
+# (10 unary + 5 binary ops).  It operates on `Node{T, D}` because the
+# constant-folding helpers compose recursively; the boundary with
+# `ScalarExpressionFunction` is `_unflatten` / `_flatten`, called once
+# per `differentiate` invocation.
 ########################################################
 
 @inline _is_const(n::Node) = n.degree == 0 && n.constant
@@ -484,7 +670,7 @@ end
 @inline _is_one(n::Node)   = _is_const(n) && isone(n.val)
 
 # Smart constructors that fold trivial constants so derivative trees stay
-# proportional in size to the input.  Each returns a fresh `Node{T, D}`.
+# proportional in size to the input.
 function _add(a::Node{T, D}, b::Node{T, D}) where {T, D}
     _is_zero(a) && return b
     _is_zero(b) && return a
@@ -516,18 +702,12 @@ function _neg(a::Node{T, D}) where {T, D}
 end
 
 function _pow_int(a::Node{T, D}, k::Int) where {T, D}
-    # Build a^k for a small positive integer constant k (>= 1).
     k == 1 && return a
     return Node{T, D}(; op = BINARY_POWER, l = a, r = Node{T, D}(; val = T(k)))
 end
 
-"""
-$(TYPEDSIGNATURES)
-
-Pure tree-rewrite differentiator.  Recurses over `node` returning a new
-`Node{T, D}` representing `∂ node / ∂ x_{var_idx}`, where `var_idx` is the
-1-based feature index of the variable to differentiate with respect to.
-"""
+# Recursive tree-rewrite differentiator.  Returns a new tree representing
+# ∂node/∂x_{var_idx}.
 function _differentiate(node::Node{T, D}, var_idx::Int) where {T, D}
     if node.degree == 0
         if node.constant
@@ -585,8 +765,6 @@ function _differentiate(node::Node{T, D}, var_idx::Int) where {T, D}
             num = _sub(_mul(du, v), _mul(u, dv))
             return _div(num, _pow_int(v, 2))
         elseif op == BINARY_POWER
-            # Common cases: constant exponent or constant base get clean
-            # derivatives.  General case uses u^v * (v' log u + v u' / u).
             if _is_const(v)
                 # d/dx(u^c) = c · u^(c-1) · u'
                 du     = _differentiate(u, var_idx)
@@ -617,31 +795,36 @@ end
 """
 $(TYPEDSIGNATURES)
 
-Return the symbolic derivative of `f` with respect to variable `var_name`.
+Return the symbolic derivative of `f` with respect to the variable whose
+1-based feature index is `var_idx`.  Differentiation is implemented as a
+recursive tree rewrite over FEC's closed grammar (10 unary + 5 binary
+operators) — no dependency on ForwardDiff, Zygote, or Symbolics.
 
-The returned function is a `ScalarExpressionFunction` over the same variable
-list as `f` — only the underlying expression tree changes.  Differentiation
-is implemented as a recursive tree rewrite over FEC's closed grammar (10
-unary + 5 binary operators), so there is no dependency on ForwardDiff,
-Zygote, or Symbolics; the result survives `juliac --trim`.
+The result is a fresh `ScalarExpressionFunction` over the same variable
+slots; the trailing variable is conventionally time.
+"""
+function differentiate(f::ScalarExpressionFunction{T}, var_idx::Integer) where T
+    @assert 1 <= Int(var_idx) <= Int(f.num_vars) "var_idx $(var_idx) out of range 1..$(Int(f.num_vars))"
+    tree        = _unflatten(f.nodes)
+    deriv_tree  = _differentiate(tree, Int(var_idx))
+    nodes, n_active = _flatten(deriv_tree)
+    return ScalarExpressionFunction{T}(nodes, n_active, f.num_vars)
+end
 
-Supported operators: cos, cosh, exp, log, sin, sinh, sqrt, tan, tanh, unary
-minus, +, -, *, /, ^.
+"""
+$(TYPEDSIGNATURES)
 
-```julia
-f = ScalarExpressionFunction{Float64}("a * exp(-(t - tc)^2 / (2 * τ^2))",
-                                      ["a", "tc", "τ", "t"])
-f_dot     = differentiate(f, "t")
-f_dot_dot = differentiate(f_dot, "t")
-```
+Convenience overload that resolves `var_name` against an explicit
+`var_names` list, then delegates to the integer form.  Useful when the
+caller still has the var-name list in scope (typical at TOML parse time);
+runtime hot paths should call the integer form directly.
 """
-function differentiate(f::ScalarExpressionFunction{T}, var_name::String) where T
-    var_names = f.expr.metadata.var_names
-    idx       = findfirst(==(var_name), var_names)
+function differentiate(f::ScalarExpressionFunction{T},
+                       var_names::AbstractVector{<:AbstractString},
+                       var_name::AbstractString) where T
+    idx = findfirst(==(var_name), var_names)
     @assert idx !== nothing "variable \"$var_name\" not in $(var_names)"
-    deriv_tree = _differentiate(f.expr.tree, idx)
-    deriv_expr = Expression(deriv_tree; operators, var_names)
-    return ScalarExpressionFunction{T}(deriv_expr, f.num_vars)
+    return differentiate(f, idx)
 end
 
 end # module
diff --git a/src/bcs/DirichletBCs.jl b/src/bcs/DirichletBCs.jl
index 662682f..8e53379 100644
--- a/src/bcs/DirichletBCs.jl
+++ b/src/bcs/DirichletBCs.jl
@@ -162,6 +162,22 @@ struct DirichletBCFunction{F1, F2, F3} <: AbstractBCFunction{F1}
     new{F, typeof(func_dot), typeof(func_dot_dot)}(func, func_dot, func_dot_dot)
   end
 
+  # Specialization for the flat, juliac-safe scalar expression: take time
+  # derivatives symbolically via [`Expressions.differentiate`](@ref) instead
+  # of ForwardDiff, so the resulting `DirichletBCFunction` is wholly isbits
+  # (no captured closures) and stays `juliac --trim` compatible.  Convention
+  # is that the last variable is time, so the derivative index is
+  # `func.num_vars`.  This is what `DirichletBCs(mesh, dof, bcs)` (the
+  # untyped, non-juliac container constructor) ends up calling when bcs
+  # carry `ScalarExpressionFunction` funcs, giving callers symbolic
+  # derivatives "for free" without needing to thread an `F` type parameter.
+  function DirichletBCFunction(func::F) where F <: Expressions.ScalarExpressionFunction
+    t_idx        = Int(func.num_vars)
+    func_dot     = Expressions.differentiate(func, t_idx)
+    func_dot_dot = Expressions.differentiate(func_dot, t_idx)
+    new{F, F, F}(func, func_dot, func_dot_dot)
+  end
+
   function DirichletBCFunction{F1, F2, F3}(
     f::F1, f_dot::F2, f_dot_dot::F3
   ) where {
@@ -264,10 +280,14 @@ struct DirichletBCs{
   end
 
   # juliac-safe path: derives `func_dot` and `func_dot_dot` symbolically via
-  # [`differentiate`](@ref) on the user's expression tree, so dynamic
-  # Dirichlet BCs work under `juliac --trim` without requiring ForwardDiff,
-  # Zygote, Symbolics, or user-supplied derivatives.  F is expected to be
-  # `Expressions.ScalarExpressionFunction{T}`.
+  # [`Expressions.differentiate`](@ref) on the user's expression tree, so
+  # dynamic Dirichlet BCs work under `juliac --trim` without requiring
+  # ForwardDiff, Zygote, Symbolics, or user-supplied derivatives.
+  #
+  # F is expected to be `Expressions.ScalarExpressionFunction{T}`.  The
+  # convention is that the last variable in the user's expression is time,
+  # so the time-derivative index is `bc.func.num_vars` (e.g. 4 for the
+  # standard `["x", "y", "z", "t"]` namespace).
   function DirichletBCs{F}(mesh::AbstractMesh, dof, bcs_input) where {F <: Function}
     bc_funcs = DirichletBCFunction{F, F, F}[]
     if length(bcs_input) == 0
@@ -278,8 +298,9 @@ struct DirichletBCs{
     end
 
     for bc in bcs_input
-      func_dot     = Expressions.differentiate(bc.func, "t")
-      func_dot_dot = Expressions.differentiate(func_dot, "t")
+      t_idx        = Int(bc.func.num_vars)        # convention: t is last
+      func_dot     = Expressions.differentiate(bc.func, t_idx)
+      func_dot_dot = Expressions.differentiate(func_dot, t_idx)
       push!(bc_funcs,
             DirichletBCFunction{F, F, F}(bc.func, func_dot, func_dot_dot))
     end
diff --git a/test/TestBCs.jl b/test/TestBCs.jl
index 6d255f9..ddad2cc 100644
--- a/test/TestBCs.jl
+++ b/test/TestBCs.jl
@@ -105,10 +105,11 @@ end
 end
 
 @testitem "BCs - juliac-safe dynamic DirichletBCs (symbolic derivatives)" setup=[BCHelper] begin
-  # Mirrors `test_dirichlet_update_bc_values!` but goes through the
-  # `DirichletBCs{F}` juliac-safe constructor, which computes the BC's first
-  # and second time derivatives symbolically via Expressions.differentiate.
-  # No ForwardDiff, no user-supplied derivatives, no Symbolics.
+  # Walks the `DirichletBCs{F}` juliac-safe constructor, which now computes
+  # the BC's first and second time derivatives symbolically via
+  # Expressions.differentiate on the flat tree.  No ForwardDiff, no user-
+  # supplied derivatives, no Symbolics — and the resulting BC function is
+  # `isbits`, so it passes through KA kernels on any backend.
   import FiniteElementContainers: update_bc_values!
   import FiniteElementContainers.Expressions: ScalarExpressionFunction
 
@@ -136,9 +137,9 @@ end
 end
 
 @testitem "BCs - juliac-safe DirichletBCs Gaussian pulse" setup=[BCHelper] begin
-  # Exercises a realistic Gaussian-pulse BC of the form used in the
-  # Norma-ported clamped-bar test: g(t) = a · exp(-(t-tc)^2 / (2 τ^2)).
-  # All three of g, g', g'' are produced by symbolic differentiation alone.
+  # Realistic Gaussian-pulse BC of the form used in the Norma-ported
+  # clamped-bar test.  All three of g, g', g'' come from symbolic
+  # differentiation on the flat tree alone.
   import FiniteElementContainers: update_bc_values!
   import FiniteElementContainers.Expressions: ScalarExpressionFunction
 
diff --git a/test/TestExpressions.jl b/test/TestExpressions.jl
index b5b280b..195fd3f 100644
--- a/test/TestExpressions.jl
+++ b/test/TestExpressions.jl
@@ -185,6 +185,15 @@ end
     @test val[3] ≈ 15.0 * exp(5.0)
 end
 
+@testitem "ScalarExpressionFunction - is isbits (GPU + juliac safe)" begin
+    import FiniteElementContainers.Expressions: ScalarExpressionFunction, FlatNode
+    @test isbitstype(FlatNode{Float64})
+    @test isbitstype(ScalarExpressionFunction{Float64})
+    f = ScalarExpressionFunction{Float64}("a*exp(-(t-tc)^2 / (2*tau^2))",
+                                          ["a", "tc", "tau", "t"])
+    @test isbits(f)
+end
+
 @testitem "ScalarExpressionFunction - parser precedence: unary minus vs ^" begin
     import FiniteElementContainers.Expressions: ScalarExpressionFunction
     # Standard math: `-t^2` parses as `-(t^2)`, not `(-t)^2`.
@@ -192,38 +201,27 @@ end
     @test f([3.0]) ≈ -9.0
     @test f([-3.0]) ≈ -9.0
 
-    # Numeric base: -2^2 = -(2^2) = -4
     g = ScalarExpressionFunction{Float64}("-2^2", String[])
     @test g(Float64[]) ≈ -4.0
 
-    # Right-side unary minus inside power: 2^-2 = 0.25
     h = ScalarExpressionFunction{Float64}("2^-2", String[])
     @test h(Float64[]) ≈ 0.25
-
-    # Multiplicative remains unaffected (commutativity hides any change)
-    p = ScalarExpressionFunction{Float64}("-x*y", ["x", "y"])
-    @test p([2.0, 3.0]) ≈ -6.0
 end
 
 @testitem "differentiate - constants and variables" begin
     import FiniteElementContainers.Expressions: ScalarExpressionFunction, differentiate
-    # d/dx of a constant is zero
     f = ScalarExpressionFunction{Float64}("3.14", ["x"])
-    fp = differentiate(f, "x")
+    fp = differentiate(f, 1)
     @test fp([1.7]) ≈ 0.0
     @test fp([-9.3]) ≈ 0.0
 
-    # d/dx of x is 1, of y is 0
     g = ScalarExpressionFunction{Float64}("x", ["x", "y"])
-    gp = differentiate(g, "x")
-    @test gp([1.0, 2.0]) ≈ 1.0
-    gq = differentiate(g, "y")
-    @test gq([1.0, 2.0]) ≈ 0.0
+    @test differentiate(g, 1)([1.0, 2.0]) ≈ 1.0
+    @test differentiate(g, 2)([1.0, 2.0]) ≈ 0.0
 end
 
 @testitem "differentiate - each unary operator" begin
     import FiniteElementContainers.Expressions: ScalarExpressionFunction, differentiate
-    # (func name, analytical derivative at x)
     cases = [
         ("cos(x)",  x -> -sin(x)),
         ("cosh(x)", x ->  sinh(x)),
@@ -237,63 +235,56 @@ end
     ]
     for (expr, dexpr) in cases
         f  = ScalarExpressionFunction{Float64}(expr, ["x"])
-        fp = differentiate(f, "x")
+        fp = differentiate(f, 1)
         for x in (0.3, 0.7, 1.4, 2.6)
             @test fp([x]) ≈ dexpr(x) rtol=1e-12
         end
     end
 
-    # Unary minus: d/dx(-x) = -1
+    # Unary minus
     f  = ScalarExpressionFunction{Float64}("-x", ["x"])
-    fp = differentiate(f, "x")
-    @test fp([5.0]) ≈ -1.0
+    @test differentiate(f, 1)([5.0]) ≈ -1.0
 end
 
 @testitem "differentiate - binary operators" begin
     import FiniteElementContainers.Expressions: ScalarExpressionFunction, differentiate
 
-    # +, -, *, /
-    f  = ScalarExpressionFunction{Float64}("x + y", ["x", "y"])
-    @test differentiate(f, "x")([2.0, 3.0]) ≈ 1.0
-    @test differentiate(f, "y")([2.0, 3.0]) ≈ 1.0
-
-    f  = ScalarExpressionFunction{Float64}("x - y", ["x", "y"])
-    @test differentiate(f, "x")([2.0, 3.0]) ≈ 1.0
-    @test differentiate(f, "y")([2.0, 3.0]) ≈ -1.0
-
-    f  = ScalarExpressionFunction{Float64}("x * y", ["x", "y"])
-    @test differentiate(f, "x")([2.0, 3.0]) ≈ 3.0
-    @test differentiate(f, "y")([2.0, 3.0]) ≈ 2.0
-
-    f  = ScalarExpressionFunction{Float64}("x / y", ["x", "y"])
-    @test differentiate(f, "x")([2.0, 3.0]) ≈ 1/3
-    @test differentiate(f, "y")([2.0, 3.0]) ≈ -2/9
-
-    # Power: constant exponent — d/dx(x^c) = c x^{c-1}
-    f  = ScalarExpressionFunction{Float64}("x^3", ["x"])
-    fp = differentiate(f, "x")
-    @test fp([2.0]) ≈ 12.0
-
-    # Power: constant base — d/dx(c^x) = c^x log(c)
-    f  = ScalarExpressionFunction{Float64}("2^x", ["x"])
-    fp = differentiate(f, "x")
-    @test fp([3.0]) ≈ 8.0 * log(2.0)
-
-    # Power: general — d/dx(x^x) = x^x (log(x) + 1)
-    f  = ScalarExpressionFunction{Float64}("x^x", ["x"])
-    fp = differentiate(f, "x")
-    @test fp([2.0]) ≈ 2.0^2.0 * (log(2.0) + 1.0)
+    f = ScalarExpressionFunction{Float64}("x + y", ["x", "y"])
+    @test differentiate(f, 1)([2.0, 3.0]) ≈ 1.0
+    @test differentiate(f, 2)([2.0, 3.0]) ≈ 1.0
+
+    f = ScalarExpressionFunction{Float64}("x - y", ["x", "y"])
+    @test differentiate(f, 1)([2.0, 3.0]) ≈ 1.0
+    @test differentiate(f, 2)([2.0, 3.0]) ≈ -1.0
+
+    f = ScalarExpressionFunction{Float64}("x * y", ["x", "y"])
+    @test differentiate(f, 1)([2.0, 3.0]) ≈ 3.0
+    @test differentiate(f, 2)([2.0, 3.0]) ≈ 2.0
+
+    f = ScalarExpressionFunction{Float64}("x / y", ["x", "y"])
+    @test differentiate(f, 1)([2.0, 3.0]) ≈ 1/3
+    @test differentiate(f, 2)([2.0, 3.0]) ≈ -2/9
+
+    # Power: constant exponent
+    f = ScalarExpressionFunction{Float64}("x^3", ["x"])
+    @test differentiate(f, 1)([2.0]) ≈ 12.0
+
+    # Power: constant base
+    f = ScalarExpressionFunction{Float64}("2^x", ["x"])
+    @test differentiate(f, 1)([3.0]) ≈ 8.0 * log(2.0)
+
+    # Power: general
+    f = ScalarExpressionFunction{Float64}("x^x", ["x"])
+    @test differentiate(f, 1)([2.0]) ≈ 2.0^2.0 * (log(2.0) + 1.0)
 end
 
 @testitem "differentiate - Gaussian pulse to 2nd derivative" begin
     import FiniteElementContainers.Expressions: ScalarExpressionFunction, differentiate
     # g(t) = a * exp(-(t - tc)^2 / (2 τ^2))
-    # g'   = -((t-tc)/τ^2) * g
-    # g''  = ((t-tc)^2/τ^4 - 1/τ^2) * g
     g   = ScalarExpressionFunction{Float64}(
         "a * exp(-(t - tc)^2 / (2 * tau^2))", ["a", "tc", "tau", "t"])
-    gp  = differentiate(g, "t")
-    gpp = differentiate(gp, "t")
+    gp  = differentiate(g, 4)
+    gpp = differentiate(gp, 4)
     a, tc, τ = 1.0e-3, 2.5e-4, 5.0e-5
     for t in (0.0, 1.0e-4, 2.5e-4, 4.0e-4, 5.0e-4)
         η  = t - tc
@@ -302,17 +293,16 @@ end
         gpp_= (η^2 / τ^4 - 1 / τ^2) * g_
         @test g([a, tc, τ, t])   ≈ g_   rtol=1e-12
         @test gp([a, tc, τ, t])  ≈ gp_  rtol=1e-12
-        @test gpp([a, tc, τ, t]) ≈ gpp_ rtol=1e-12
+        @test gpp([a, tc, τ, t]) ≈ gpp_ rtol=1e-9
     end
 end
 
-@testitem "differentiate - spatial derivatives (traveling wave IC)" begin
+@testitem "differentiate - spatial (traveling-wave IC)" begin
     import FiniteElementContainers.Expressions: ScalarExpressionFunction, differentiate
-    # u₀(z) = a * exp(-z^2 / (2 s^2))
-    # u₀'(z) = -(z/s^2) * u₀
+    # u₀(z) = a * exp(-z^2 / (2 s^2));  ∂u/∂z = -(z/s²) u
     u0  = ScalarExpressionFunction{Float64}(
         "a * exp(-z^2 / (2 * s^2))", ["a", "s", "z"])
-    duz = differentiate(u0, "z")
+    duz = differentiate(u0, 3)
     a, s = 0.01, 0.02
     for z in (-0.04, -0.01, 0.0, 0.01, 0.04)
         u_  = a * exp(-z^2 / (2 * s^2))
@@ -322,8 +312,11 @@ end
     end
 end
 
-@testitem "differentiate - error on unknown variable" begin
+@testitem "differentiate - var-name overload + error on unknown" begin
     import FiniteElementContainers.Expressions: ScalarExpressionFunction, differentiate
-    f = ScalarExpressionFunction{Float64}("x", ["x"])
-    @test_throws AssertionError differentiate(f, "y")
+    f = ScalarExpressionFunction{Float64}("x + 2*t", ["x", "t"])
+    fp = differentiate(f, ["x", "t"], "t")
+    @test fp([0.5, 0.3]) ≈ 2.0
+    @test_throws AssertionError differentiate(f, ["x", "t"], "y")
+    @test_throws AssertionError differentiate(f, 5)
 end