diff --git a/src/methods/jacobiMethodScript.js b/src/methods/jacobiMethodScript.js
index 1a3e3df..bebe702 100644
--- a/src/methods/jacobiMethodScript.js
+++ b/src/methods/jacobiMethodScript.js
@@ -22,6 +22,36 @@
  */
 export function jacobiMethod(A, b, x0, maxIterations = 100, tolerance = 1e-7) {
   const n = A.length; // Size of the square matrix
+
+  // Sanity checks for input dimensions
+  if (!Array.isArray(A) || n === 0) {
+    throw new Error("Matrix A must be a non-empty array");
+  }
+
+  // Verify A is square
+  for (let i = 0; i < n; i++) {
+    if (!Array.isArray(A[i]) || A[i].length !== n) {
+      throw new Error(`Matrix A must be square. Row ${i} has length ${A[i].length}, expected ${n}`);
+    }
+  }
+
+  // Verify b is a vector of correct length
+  if (!Array.isArray(b) || b.length !== n) {
+    throw new Error(`Vector b must have length ${n}, got ${b.length}`);
+  }
+
+  // Verify x0 is a vector of correct length
+  if (!Array.isArray(x0) || x0.length !== n) {
+    throw new Error(`Initial guess x0 must have length ${n}, got ${x0.length}`);
+  }
+
+  // Verify no zero diagonal elements (required for Jacobi method)
+  for (let i = 0; i < n; i++) {
+    if (A[i][i] === 0) {
+      throw new Error(`Diagonal element A[${i}][${i}] is zero; Jacobi method requires non-zero diagonal elements`);
+    }
+  }
+
   let x = [...x0]; // Current solution (starts with initial guess)
   let xNew = new Array(n); // Next iteration's solution
 
diff --git a/tests/regression/HeatConduction1DWall/REGRESSION.md b/tests/regression/HeatConduction1DWall/REGRESSION.md
new file mode 100644
index 0000000..579ffac
--- /dev/null
+++ b/tests/regression/HeatConduction1DWall/REGRESSION.md
@@ -0,0 +1,61 @@
+# Regression Test — HeatConduction1DWall
+
+## Purpose
+
+This test guards the numerical output of the 1D heat-conduction-through-a-wall example
+against unintended changes to the solver, assembler, or mesh-generation logic.
+
+It replicates exactly the problem set up in
+[`HeatConduction1DWall.html`](../../../examples/solidHeatTransferScript/HeatConduction1DWall/HeatConduction1DWall.html)
+and asserts a known good value.
+
+## Problem setup
+
+| Parameter | Value |
+|-----------|-------|
+| Domain | 1D, 0 – 0.15 m |
+| Mesh | 10 linear elements |
+| Boundary 0 (x = 0) | Convection, h = 1, T∞ = 25 °C |
+| Boundary 1 (x = 0.15) | Constant temperature, T = 5 °C |
+| Solver | LU decomposition (`lusolve`) |
+
+## Expected value
+
+| Quantity | Value |
+|----------|-------|
+| Temperature at node 0 (x = 0) | **10.29412 °C** |
+
+Tolerance used in the assertion: `1e-4`.
+
+## How to run
+
+From the repository root:
+
+```bash
+node tests/regression/HeatConduction1DWall/regression.test.js
+```
+
+A passing run prints:
+
+```
+PASS: T(x=0) = 10.29412 (expected 10.29412)
+```
+
+A failing run prints a `FAIL:` message and exits with code 1.
+
+The `test` script in `package.json` also runs this file, so `npm test` works too.
+
+## After modifying the code
+
+| Situation | Action |
+|-----------|--------|
+| Bug fix that should not change results | Run the test — it must still pass. |
+| Intentional algorithm change (new element type, new integration rule, etc.) | Re-derive the expected value, update `EXPECTED_T0` in `regression.test.js`, and document the reason here. |
+| New boundary condition API | Update both the test and the reference HTML example together. |
+| Adding a new solver method | Add a separate assertion block for the new method; keep the `lusolve` block untouched as the baseline. |
+
+## Change log
+
+| Date | Change | New expected value |
+|------|--------|--------------------|
+| 2026-06-14 | Initial regression baseline | 10.29412 |
diff --git a/tests/regression/HeatConduction1DWall/regression.test.js b/tests/regression/HeatConduction1DWall/regression.test.js
new file mode 100644
index 0000000..fbd5f2b
--- /dev/null
+++ b/tests/regression/HeatConduction1DWall/regression.test.js
@@ -0,0 +1,61 @@
+/**
+ * Regression test for HeatConduction1DWall
+ *
+ * Replicates the exact setup from HeatConduction1DWall.html and asserts
+ * that the temperature at node 0 (convection boundary) remains 10.29412.
+ *
+ * Run: node tests/regression/HeatConduction1DWall/regression.test.js
+ */
+
+import * as mathjs from "mathjs";
+import { FEAScriptModel } from "../../../src/FEAScript.js";
+import { basicLog } from "../../../src/utilities/loggingScript.js";
+
+// FEAScript.js references `math` as a global (loaded via CDN in browser).
+// Set it here before any solve() call.
+globalThis.math = mathjs;
+
+const EXPECTED_T0 = 10.29412;
+const TOLERANCE = 1e-4;
+
+function runSimulation() {
+  const model = new FEAScriptModel();
+
+  basicLog("")
+  basicLog("================================")
+  basicLog("Starting test in solid heat transfer 1D wall...")
+
+  model.setSolverConfig("solidHeatTransferScript");
+  model.setMeshConfig({
+    meshDimension: "1D",
+    elementOrder: "linear",
+    numElementsX: 10,
+    maxX: 0.15,
+  });
+
+  model.addBoundaryCondition("0", ["convection", 1, 25]);
+  model.addBoundaryCondition("1", ["constantTemp", 5]);
+  model.setSolverMethod("lusolve");
+
+  return model.solve();
+}
+
+function assert(condition, message) {
+  if (!condition) {
+    console.error(`FAIL: ${message}`);
+    process.exit(1);
+  }
+}
+
+const { solutionVector } = runSimulation();
+
+// solutionVector from math.lusolve is a nested array: [[T0], [T1], ...]
+const T0 = Array.isArray(solutionVector[0]) ? solutionVector[0][0] : solutionVector[0];
+
+assert(
+  Math.abs(T0 - EXPECTED_T0) < TOLERANCE,
+  `Temperature at node 0: expected ${EXPECTED_T0}, got ${T0} (tolerance ${TOLERANCE})`
+);
+
+console.log(`PASS: T(x=0) = ${T0.toFixed(5)} (expected ${EXPECTED_T0})`);
+basicLog("================================")
diff --git a/tests/regression/HeatConduction2DFin/REGRESSION.md b/tests/regression/HeatConduction2DFin/REGRESSION.md
new file mode 100644
index 0000000..5bea510
--- /dev/null
+++ b/tests/regression/HeatConduction2DFin/REGRESSION.md
@@ -0,0 +1,67 @@
+# Regression Test — HeatConduction2DFin
+
+## Purpose
+
+This test guards the numerical output of the 2D heat-conduction-in-a-fin example
+against unintended changes to the solver, assembler, or mesh-generation logic.
+
+It replicates exactly the problem set up in
+[`HeatConduction2DFin.html`](../../../examples/solidHeatTransferScript/HeatConduction2DFin/HeatConduction2DFin.html)
+and asserts a known good value at a representative interior point.
+
+## Problem setup
+
+| Parameter | Value |
+|-----------|-------|
+| Domain | 2D, x ∈ [0, 4] m, y ∈ [0, 2] m |
+| Mesh | 8 × 4 quadratic elements |
+| Boundary 0 (bottom, y = 0) | Constant temperature, T = 200 °C |
+| Boundary 1 (left, x = 0) | Symmetry (zero flux) |
+| Boundary 2 (top, y = 2) | Convection, h = 1, T∞ = 20 °C |
+| Boundary 3 (right, x = 4) | Constant temperature, T = 200 °C |
+| Solver | LU decomposition (`lusolve`) |
+
+## Expected value
+
+| Quantity | Value |
+|----------|-------|
+| Temperature at node (x = 0, y = 2) | **81.31873 °C** |
+
+This point sits at the top-left corner of the fin — on the symmetry boundary and
+the convection boundary — and is sensitive to both the heat transfer coefficient
+and the thermal gradient across the domain.
+
+Tolerance used in the assertion: `1e-4`.
+
+## How to run
+
+From the repository root:
+
+```bash
+node tests/regression/HeatConduction2DFin/regression.test.js
+```
+
+A passing run prints:
+
+```
+PASS: T(x=0, y=2) = 81.31873 (expected 81.31873)
+```
+
+A failing run prints a `FAIL:` message and exits with code 1.
+
+Running `npm test` executes all regression tests, including this one.
+
+## After modifying the code
+
+| Situation | Action |
+|-----------|--------|
+| Bug fix that should not change results | Run the test — it must still pass. |
+| Intentional algorithm change (new element type, new integration rule, etc.) | Re-derive the expected value, update `EXPECTED_T` in `regression.test.js`, and document the reason here. |
+| New boundary condition API | Update both the test and the reference HTML example together. |
+| Mesh refinement study | Add a separate assertion block for the refined mesh; keep the current block as the coarse-mesh baseline. |
+
+## Change log
+
+| Date | Change | New expected value |
+|------|--------|--------------------|
+| 2026-06-14 | Initial regression baseline | 81.31873 |
diff --git a/tests/regression/HeatConduction2DFin/regression.test.js b/tests/regression/HeatConduction2DFin/regression.test.js
new file mode 100644
index 0000000..b57d48e
--- /dev/null
+++ b/tests/regression/HeatConduction2DFin/regression.test.js
@@ -0,0 +1,77 @@
+/**
+ * Regression test for HeatConduction2DFin
+ *
+ * Replicates the exact setup from HeatConduction2DFin.html and asserts
+ * that the temperature at (x=0, y=2) remains 81.31873.
+ *
+ * Run: node tests/regression/HeatConduction2DFin/regression.test.js
+ */
+
+import * as mathjs from "mathjs";
+import { FEAScriptModel } from "../../../src/FEAScript.js";
+import { basicLog } from "../../../src/utilities/loggingScript.js";
+
+basicLog("")
+basicLog("================================")
+basicLog("Starting test in solid heat transfer 2D fin...")
+
+// FEAScript.js references `math` as a global (loaded via CDN in browser).
+// Set it here before any solve() call.
+globalThis.math = mathjs;
+
+const EXPECTED_X = 0;
+const EXPECTED_Y = 2;
+const EXPECTED_T = 81.31873;
+const TOLERANCE = 1e-4;
+
+function runSimulation() {
+  const model = new FEAScriptModel();
+
+  model.setSolverConfig("solidHeatTransferScript");
+  model.setMeshConfig({
+    meshDimension: "2D",
+    elementOrder: "quadratic",
+    numElementsX: 8,
+    numElementsY: 4,
+    maxX: 4,
+    maxY: 2,
+  });
+
+  model.addBoundaryCondition("0", ["constantTemp", 200]);
+  model.addBoundaryCondition("1", ["symmetry"]);
+  model.addBoundaryCondition("2", ["convection", 1, 20]);
+  model.addBoundaryCondition("3", ["constantTemp", 200]);
+  model.setSolverMethod("lusolve");
+
+  return model.solve();
+}
+
+function assert(condition, message) {
+  if (!condition) {
+    console.error(`FAIL: ${message}`);
+    process.exit(1);
+  }
+}
+
+const { solutionVector, nodesCoordinates } = runSimulation();
+const { nodesXCoordinates, nodesYCoordinates } = nodesCoordinates;
+
+// Locate the node at (x=0, y=2) by searching coordinates.
+// This is robust against changes in mesh ordering conventions.
+const nodeIndex = nodesXCoordinates.findIndex(
+  (x, i) => Math.abs(x - EXPECTED_X) < 1e-10 && Math.abs(nodesYCoordinates[i] - EXPECTED_Y) < 1e-10
+);
+
+assert(nodeIndex !== -1, `No node found at (x=${EXPECTED_X}, y=${EXPECTED_Y})`);
+
+// solutionVector from math.lusolve is a nested array: [[T0], [T1], ...]
+const T = Array.isArray(solutionVector[nodeIndex]) ? solutionVector[nodeIndex][0] : solutionVector[nodeIndex];
+
+assert(
+  Math.abs(T - EXPECTED_T) < TOLERANCE,
+  `Temperature at (x=${EXPECTED_X}, y=${EXPECTED_Y}): expected ${EXPECTED_T}, got ${T} (tolerance ${TOLERANCE})`
+);
+
+console.log(`PASS: T(x=${EXPECTED_X}, y=${EXPECTED_Y}) = ${T.toFixed(5)} (expected ${EXPECTED_T})`);
+
+basicLog("================================")
\ No newline at end of file
diff --git a/tests/run-all-tests.js b/tests/run-all-tests.js
new file mode 100644
index 0000000..8901fcc
--- /dev/null
+++ b/tests/run-all-tests.js
@@ -0,0 +1,53 @@
+/**
+ * Test runner — discovers and executes every *.test.js file under tests/.
+ * Add a new test file anywhere in this tree and it runs automatically.
+ *
+ * Usage: node tests/run-all-tests.js
+ */
+
+import { readdirSync, statSync } from "fs";
+import { join, relative } from "path";
+import { spawnSync } from "child_process";
+import { fileURLToPath } from "url";
+
+const __dirname = fileURLToPath(new URL(".", import.meta.url));
+
+function collectTestFiles(dir) {
+  const entries = readdirSync(dir);
+  const files = [];
+  for (const entry of entries) {
+    const fullPath = join(dir, entry);
+    if (statSync(fullPath).isDirectory()) {
+      files.push(...collectTestFiles(fullPath));
+    } else if (entry.endsWith(".test.js")) {
+      files.push(fullPath);
+    }
+  }
+  return files;
+}
+
+const testFiles = collectTestFiles(__dirname);
+
+if (testFiles.length === 0) {
+  console.log("No test files found.");
+  process.exit(0);
+}
+
+console.log(`Found ${testFiles.length} test file(s).\n`);
+
+let passed = 0;
+let failed = 0;
+
+for (const file of testFiles) {
+  const label = relative(__dirname, file);
+  const result = spawnSync(process.execPath, [file], { stdio: "inherit" });
+  if (result.status === 0) {
+    passed++;
+  } else {
+    console.error(`\nFAILED: ${label}\n`);
+    failed++;
+  }
+}
+
+console.log(`\n${passed} passed, ${failed} failed.`);
+process.exit(failed > 0 ? 1 : 0);
diff --git a/tests/unit/README.md b/tests/unit/README.md
new file mode 100644
index 0000000..dee4e93
--- /dev/null
+++ b/tests/unit/README.md
@@ -0,0 +1,9 @@
+# Unit Tests
+
+This folder will contain unit tests for individual FEAScript modules (solvers, assemblers, utilities, etc.).
+
+Each test file should target a single module and be runnable with:
+
+```bash
+node tests/unit/<test-file>.js
+```