@@ -28,7 +28,19 @@ export type TapeOp =
   | { type: "tan"; a: number }
   | { type: "exp"; a: number }
   | { type: "log"; a: number }
-  | { type: "sqrt"; a: number };
+  | { type: "sqrt"; a: number }
+  // ── Vector ops: SIMD-style operations on contiguous blocks ──
+  // These occupy `size` consecutive tape slots starting at the op's index.
+  | { type: "vec_var"; baseName: string; size: number }
+  | { type: "vec_const"; vals: number[]; size: number }
+  | { type: "vec_add"; a: number; b: number; size: number }
+  | { type: "vec_sub"; a: number; b: number; size: number }
+  | { type: "vec_mul"; a: number; b: number; size: number }
+  | { type: "vec_neg"; a: number; size: number }
+  // Extracts scalar t[a + offset] from a vector block starting at tape index `a`
+  | { type: "vec_subscript"; a: number; offset: number }
+  // No-op placeholder for vector padding slots (evaluators skip these)
+  | { type: "nop" };
 
 function formatCDouble(v: number): string {
   if (!isFinite(v)) return v === Infinity ? "INFINITY" : v === -Infinity ? "(-INFINITY)" : "NAN";
@@ -187,6 +199,69 @@ export class StaticTapeBuilder {
     return [this.walk(expr)];
   }
 
+  /**
+   * Push a vector op that reserves `size` consecutive tape slots.
+   * The op itself occupies slot `idx`, and `size-1` placeholder slots follow.
+   * Returns the starting index of the vector block.
+   */
+  public pushVecOp(op: TapeOp & { size: number }): number {
+    const idx = this.ops.length;
+    this.ops.push(op);
+    // Reserve remaining size-1 slots with nop placeholders
+    for (let i = 1; i < op.size; i++) {
+      this.ops.push({ type: "nop" }); // Placeholder; evaluators skip these
+    }
+    return idx;
+  }
+
+  /**
+   * Walk an array expression as a single vectorized block.
+   * If `expr` is a `ModelicaArray` of named variables, emits `vec_var` etc.
+   * Returns the starting tape index and the block size.
+   *
+   * For array binary ops like `{a,b} + {c,d}`, emits `vec_add`.
+   * Falls back to element-wise scalar tapes for complex expressions.
+   */
+  public walkArrayVectorized(expr: ModelicaExpression): { startIdx: number; size: number } {
+    if (!(expr instanceof ModelicaArray)) {
+      // Scalar: single tape slot
+      const idx = this.walk(expr);
+      return { startIdx: idx, size: 1 };
+    }
+
+    const flatElems = [...expr.flatElements];
+    const size = flatElems.length;
+    if (size === 0) return { startIdx: this.pushOp({ type: "const", val: 0 }), size: 1 };
+
+    // Check if all elements are simple named variables with a common base
+    const allNamed = flatElems.every((e) => e instanceof ModelicaNameExpression);
+    if (allNamed && size > 1) {
+      // Extract common base name (e.g., "x[1]", "x[2]" → base = "x")
+      const firstName = (flatElems[0] as ModelicaNameExpression).name;
+      const bracketPos = firstName.indexOf("[");
+      const baseName = bracketPos >= 0 ? firstName.substring(0, bracketPos) : firstName;
+      const startIdx = this.pushVecOp({ type: "vec_var", baseName, size });
+      return { startIdx, size };
+    }
+
+    // Check if all elements are constants
+    const allConst = flatElems.every((e) => e instanceof ModelicaRealLiteral || e instanceof ModelicaIntegerLiteral);
+    if (allConst && size > 1) {
+      const vals = flatElems.map((e) =>
+        e instanceof ModelicaRealLiteral ? e.value : e instanceof ModelicaIntegerLiteral ? e.value : 0,
+      );
+      const startIdx = this.pushVecOp({ type: "vec_const", vals, size });
+      return { startIdx, size };
+    }
+
+    // Fallback: walk each element as scalar, pack into contiguous block
+    const startIdx = this.ops.length;
+    for (const elem of flatElems) {
+      this.walk(elem);
+    }
+    return { startIdx, size };
+  }
+
   /**
    * Emit the C-code evaluating the expressions step-by-step.
    * `varResolver` maps the Modelica variable name into a valid C-code getter string
@@ -264,6 +339,32 @@ export class StaticTapeBuilder {
         case "sqrt":
           rhs = `sqrt(t[${op.a}])`;
           break;
+        // ── Vector ops: emit for-loops ──
+        case "vec_var":
+          lines.push(`  for (int _k = 0; _k < ${op.size}; _k++) t[${i}+_k] = ${varResolver(`${op.baseName}[_k+1]`)};`);
+          continue; // Skip default t[i] = rhs line
+        case "vec_const":
+          for (let k = 0; k < op.size; k++) {
+            lines.push(`  t[${i + k}] = ${formatCDouble(op.vals[k] ?? 0)};`);
+          }
+          continue;
+        case "vec_add":
+          lines.push(`  for (int _k = 0; _k < ${op.size}; _k++) t[${i}+_k] = t[${op.a}+_k] + t[${op.b}+_k];`);
+          continue;
+        case "vec_sub":
+          lines.push(`  for (int _k = 0; _k < ${op.size}; _k++) t[${i}+_k] = t[${op.a}+_k] - t[${op.b}+_k];`);
+          continue;
+        case "vec_mul":
+          lines.push(`  for (int _k = 0; _k < ${op.size}; _k++) t[${i}+_k] = t[${op.a}+_k] * t[${op.b}+_k];`);
+          continue;
+        case "vec_neg":
+          lines.push(`  for (int _k = 0; _k < ${op.size}; _k++) t[${i}+_k] = -t[${op.a}+_k];`);
+          continue;
+        case "vec_subscript":
+          rhs = `t[${op.a + op.offset}]`;
+          break;
+        case "nop":
+          continue; // Skip padding slots
       }
       lines.push(`  t[${i}] = ${rhs};`);
     }
 
@@ -78,6 +78,42 @@ export function evaluateTapeForward(ops: TapeOp[], varValues: Map<string, number
       case "sqrt":
         t[i] = Math.sqrt(t[op.a] ?? 0);
         break;
+      // ── Vector ops ──
+      case "vec_var":
+        for (let k = 0; k < op.size; k++) {
+          t[i + k] = varValues.get(`${op.baseName}[${k + 1}]`) ?? 0;
+        }
+        break;
+      case "vec_const":
+        for (let k = 0; k < op.size; k++) {
+          t[i + k] = op.vals[k] ?? 0;
+        }
+        break;
+      case "vec_add":
+        for (let k = 0; k < op.size; k++) {
+          t[i + k] = (t[op.a + k] ?? 0) + (t[op.b + k] ?? 0);
+        }
+        break;
+      case "vec_sub":
+        for (let k = 0; k < op.size; k++) {
+          t[i + k] = (t[op.a + k] ?? 0) - (t[op.b + k] ?? 0);
+        }
+        break;
+      case "vec_mul":
+        for (let k = 0; k < op.size; k++) {
+          t[i + k] = (t[op.a + k] ?? 0) * (t[op.b + k] ?? 0);
+        }
+        break;
+      case "vec_neg":
+        for (let k = 0; k < op.size; k++) {
+          t[i + k] = -(t[op.a + k] ?? 0);
+        }
+        break;
+      case "vec_subscript":
+        t[i] = t[op.a + op.offset] ?? 0;
+        break;
+      case "nop":
+        break;
     }
   }
   return t;
@@ -140,6 +176,38 @@ export function evaluateTapeReverse(ops: TapeOp[], t: Float64Array, outputIndex:
       case "sqrt":
         dt[op.a] = (dt[op.a] ?? 0) + dti / (2 * (t[i] ?? 1));
         break;
+      // ── Vector ops reverse ──
+      case "vec_add":
+        for (let k = 0; k < op.size; k++) {
+          const dk = dt[i + k] ?? 0;
+          dt[op.a + k] = (dt[op.a + k] ?? 0) + dk;
+          dt[op.b + k] = (dt[op.b + k] ?? 0) + dk;
+        }
+        break;
+      case "vec_sub":
+        for (let k = 0; k < op.size; k++) {
+          const dk = dt[i + k] ?? 0;
+          dt[op.a + k] = (dt[op.a + k] ?? 0) + dk;
+          dt[op.b + k] = (dt[op.b + k] ?? 0) - dk;
+        }
+        break;
+      case "vec_mul":
+        for (let k = 0; k < op.size; k++) {
+          const dk = dt[i + k] ?? 0;
+          dt[op.a + k] = (dt[op.a + k] ?? 0) + dk * (t[op.b + k] ?? 0);
+          dt[op.b + k] = (dt[op.b + k] ?? 0) + dk * (t[op.a + k] ?? 0);
+        }
+        break;
+      case "vec_neg":
+        for (let k = 0; k < op.size; k++) {
+          dt[op.a + k] = (dt[op.a + k] ?? 0) - (dt[i + k] ?? 0);
+        }
+        break;
+      case "vec_subscript":
+        dt[op.a + op.offset] = (dt[op.a + op.offset] ?? 0) + dti;
+        break;
+      case "nop":
+        break;
     }
   }
 
@@ -149,6 +217,11 @@ export function evaluateTapeReverse(ops: TapeOp[], t: Float64Array, outputIndex:
     const op = ops[i]!; // eslint-disable-line @typescript-eslint/no-non-null-assertion
     if (op.type === "var") {
       gradients.set(op.name, (gradients.get(op.name) ?? 0) + (dt[i] ?? 0));
+    } else if (op.type === "vec_var") {
+      for (let k = 0; k < op.size; k++) {
+        const name = `${op.baseName}[${k + 1}]`;
+        gradients.set(name, (gradients.get(name) ?? 0) + (dt[i + k] ?? 0));
+      }
     }
   }
   return gradients;
 
@@ -395,3 +395,28 @@ export function buildSbbFromDAE(
     constraintTapes,
   };
 }
+
+/**
+ * Expand array variable bounds in a DomainBox.
+ * Given a box with `x → [lo, hi]` and a variable `x` with arrayDimensions = [3],
+ * produces `x[1] → [lo,hi]`, `x[2] → [lo,hi]`, `x[3] → [lo,hi]`.
+ *
+ * @param box       Initial domain box (may contain array-valued variable names)
+ * @param dae       DAE for variable metadata (arrayDimensions)
+ * @returns Expanded domain box with per-element entries
+ */
+export function expandArrayBounds(box: DomainBox, dae: ModelicaDAE): DomainBox {
+  const expanded: DomainBox = new Map();
+  for (const [name, interval] of box) {
+    const v = dae.variables.find((dv) => dv.name === name);
+    if (v?.arrayDimensions && v.arrayDimensions.length > 0) {
+      const size = v.arrayDimensions.reduce((a: number, b: number) => a * b, 1);
+      for (let i = 0; i < size; i++) {
+        expanded.set(`${name}[${i + 1}]`, new Interval(interval.lo, interval.hi));
+      }
+    } else {
+      expanded.set(name, interval);
+    }
+  }
+  return expanded;
+}
@@ -17,6 +17,8 @@
 import { StaticTapeBuilder } from "./ad-codegen.js";
 import type { ModelicaDAE, ModelicaExpression } from "./dae.js";
 import {
+  ModelicaArray,
+  ModelicaArrayEquation,
   ModelicaBinaryExpression,
   ModelicaBooleanLiteral,
   ModelicaFunctionCallEquation,
@@ -618,12 +620,24 @@ function generateInitializeSolve(id: string, dae: ModelicaDAE, result: FmuResult
     return lines;
   }
 
-  // Collect initial equations that are simple equations (LHS = RHS)
+  // Collect initial equations — unroll array equations into per-element scalar equations
   const initEqs: { lhs: ModelicaExpression; rhs: ModelicaExpression }[] = [];
   for (const eq of dae.initialEquations) {
     if ("expression1" in eq && "expression2" in eq) {
       const se = eq as { expression1: ModelicaExpression; expression2: ModelicaExpression };
-      initEqs.push({ lhs: se.expression1, rhs: se.expression2 });
+      if (eq instanceof ModelicaArrayEquation) {
+        // Unroll array equation into per-element scalar equations
+        const lhsElems = se.expression1 instanceof ModelicaArray ? [...se.expression1.flatElements] : [se.expression1];
+        const rhsElems = se.expression2 instanceof ModelicaArray ? [...se.expression2.flatElements] : [se.expression2];
+        const n = Math.max(lhsElems.length, rhsElems.length);
+        for (let i = 0; i < n; i++) {
+          const lhs = lhsElems[i] ?? lhsElems[0];
+          const rhs = rhsElems[i] ?? rhsElems[0];
+          if (lhs && rhs) initEqs.push({ lhs, rhs });
+        }
+      } else {
+        initEqs.push({ lhs: se.expression1, rhs: se.expression2 });
+      }
     }
   }
 
@@ -2193,6 +2207,13 @@ function extractDerName(expr: unknown): string | null {
 /** Recursively collect all variable names referenced in an expression. */
 function collectReferencedNames(expr: unknown, names: Set<string>): void {
   if (!expr || typeof expr !== "object") return;
+  // Handle ModelicaArray — recurse into each element
+  if (expr instanceof ModelicaArray) {
+    for (const elem of expr.flatElements) {
+      collectReferencedNames(elem, names);
+    }
+    return;
+  }
   if ("name" in expr) {
     const nameVal = (expr as { name: unknown }).name;
     if (typeof nameVal === "string") {
 
@@ -12,6 +12,8 @@
 import { StaticTapeBuilder, type TapeOp } from "./ad-codegen.js";
 import { evaluateTapeForward, evaluateTapeReverse } from "./ad-jacobian.js";
 import {
+  ModelicaArray,
+  ModelicaArrayEquation,
   ModelicaBooleanLiteral,
   type ModelicaDAE,
   type ModelicaEquation,
@@ -58,6 +60,13 @@ function collectVarNames(expr: ModelicaExpression, names: Set<string>): void {
     names.add(expr.name);
     return;
   }
+  // Recurse into array elements
+  if (expr instanceof ModelicaArray) {
+    for (const elem of expr.flatElements) {
+      collectVarNames(elem, names);
+    }
+    return;
+  }
   if ("name" in expr) names.add((expr as { name: string }).name);
   if ("operand" in expr) collectVarNames((expr as { operand: ModelicaExpression }).operand, names);
   if ("operand1" in expr) collectVarNames((expr as { operand1: ModelicaExpression }).operand1, names);
@@ -283,10 +292,38 @@ export function solveInitialEquations(
   // Classify equations
   const classified: ClassifiedEq[] = [];
   for (const eq of dae.initialEquations) {
+    if (eq instanceof ModelicaArrayEquation) {
+      // Unroll array initial equation into per-element implicit equations
+      const se = eq as { expression1: ModelicaExpression; expression2: ModelicaExpression };
+      const lhsElems = se.expression1 instanceof ModelicaArray ? [...se.expression1.flatElements] : [se.expression1];
+      const rhsElems = se.expression2 instanceof ModelicaArray ? [...se.expression2.flatElements] : [se.expression2];
+      const n = Math.max(lhsElems.length, rhsElems.length);
+      for (let i = 0; i < n; i++) {
+        const lhs = lhsElems[i] ?? lhsElems[0];
+        const rhs = rhsElems[i] ?? rhsElems[0];
+        if (lhs && rhs) {
+          // Create a synthetic simple equation for each element
+          const elemEq = { expression1: lhs, expression2: rhs } as unknown as ModelicaEquation;
+          const c = classifyInitialEquation(elemEq, unknowns);
+          if (c) classified.push(c);
+        }
+      }
+      continue;
+    }
     const c = classifyInitialEquation(eq, unknowns);
     if (c) classified.push(c);
   }
 
+  // Expand array variable unknowns: for vars with arrayDimensions, add subscripted names
+  for (const v of dae.variables) {
+    if (v.arrayDimensions && v.arrayDimensions.length > 0 && unknowns.has(v.name)) {
+      const size = v.arrayDimensions.reduce((a: number, b: number) => a * b, 1);
+      for (let i = 0; i < size; i++) {
+        unknowns.add(`${v.name}[${i + 1}]`);
+      }
+    }
+  }
+
   // Phase 1: Solve explicit equations first (direct assignment)
   const explicitEqs = classified.filter((c): c is ExplicitEq => c.type === "explicit");
   const implicitEqs = classified.filter((c): c is ImplicitEq => c.type === "implicit");