diff --git a/tools/test_verify_lib.py b/tools/test_verify_lib.py
new file mode 100644
index 0000000..4292df5
--- /dev/null
+++ b/tools/test_verify_lib.py
@@ -0,0 +1,87 @@
+import unittest
+
+import numpy as np
+
+from verify_lib import compute_l1_error_fvm, simple_quad_2d, simple_quad_3d
+
+
+def cell_integral_power(bounds: list[tuple[float, float]]) -> float:
+    result = 1.0
+    for lower, upper in bounds:
+        result *= (upper**2 - lower**2) * 0.5
+    return result
+
+
+class ComputeL1ErrorFVMTest(unittest.TestCase):
+    def test_2d_uses_full_tensor_product_for_cell_integrals(self) -> None:
+        centers = np.array([0.25, 0.75])
+        xg, yg = np.meshgrid(centers, centers, indexing="ij")
+        x_num = np.array([xg.ravel(), yg.ravel()])
+
+        numerical = []
+        for x, y in zip(x_num[0], x_num[1]):
+            bounds = [
+                (x - 0.25, x + 0.25),
+                (y - 0.25, y + 0.25),
+            ]
+            numerical.append(cell_integral_power(bounds) / 0.25)
+
+        error = compute_l1_error_fvm(
+            x_num,
+            np.array(numerical),
+            lambda xy: xy[0] * xy[1],
+            2,
+        )
+
+        self.assertAlmostEqual(error, 0.0)
+
+    def test_3d_uses_full_tensor_product_for_cell_integrals(self) -> None:
+        centers = np.array([0.25, 0.75])
+        xg, yg, zg = np.meshgrid(centers, centers, centers, indexing="ij")
+        x_num = np.array([xg.ravel(), yg.ravel(), zg.ravel()])
+
+        numerical = []
+        for x, y, z in zip(x_num[0], x_num[1], x_num[2]):
+            bounds = [
+                (x - 0.25, x + 0.25),
+                (y - 0.25, y + 0.25),
+                (z - 0.25, z + 0.25),
+            ]
+            numerical.append(cell_integral_power(bounds) / 0.125)
+
+        error = compute_l1_error_fvm(
+            x_num,
+            np.array(numerical),
+            lambda xyz: xyz[0] * xyz[1] * xyz[2],
+            3,
+        )
+
+        self.assertAlmostEqual(error, 0.0)
+
+
+class SimpleQuadTensorProductTest(unittest.TestCase):
+    def test_2d_asymmetric_box_and_integrand(self) -> None:
+        got = simple_quad_2d(
+            lambda c: np.sin(c[0] + 2.0 * c[1]), 0.0, 1.0, 0.0, 2.0
+        )
+        expected = (np.sin(1.0) - np.sin(5.0) + np.sin(4.0)) * 0.5
+        self.assertAlmostEqual(got, expected)
+
+    def test_3d_asymmetric_box_and_integrand(self) -> None:
+        got = simple_quad_3d(
+            lambda c: np.sin(c[0] + 2.0 * c[1]) * np.exp(0.5 * c[2]),
+            0.0,
+            1.0,
+            0.0,
+            2.0,
+            0.0,
+            3.0,
+        )
+        expected = (np.sin(1.0) - np.sin(5.0) + np.sin(4.0)) * (
+            np.exp(1.5) - 1.0
+        )
+        self.assertAlmostEqual(got, expected)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/tools/verify_lib.py b/tools/verify_lib.py
index f54801f..821a87e 100644
--- a/tools/verify_lib.py
+++ b/tools/verify_lib.py
@@ -94,6 +94,61 @@ def transform(x):
     return np.sum(f(transform(points)) * weights) * 0.5 * (x1 - x0)
 
 
+def simple_quad_2d(
+    f: Callable, x0: float, x1: float, y0: float, y1: float, deg: int = 10
+) -> float:
+    """
+    Use a tensor-product Gauss-Legendre quadrature of degree deg
+    over the rectangle [x0, x1] x [y0, y1]
+    """
+
+    points, weights = np.polynomial.legendre.leggauss(deg)
+
+    xs = ((x1 - x0) * points + x1 + x0) * 0.5
+    ys = ((y1 - y0) * points + y1 + y0) * 0.5
+    xg, yg = np.meshgrid(xs, ys, indexing="ij")
+    wg = np.outer(weights, weights).ravel()
+
+    return (
+        np.sum(f([xg.ravel(), yg.ravel()]) * wg)
+        * 0.25
+        * (x1 - x0)
+        * (y1 - y0)
+    )
+
+
+def simple_quad_3d(
+    f: Callable,
+    x0: float,
+    x1: float,
+    y0: float,
+    y1: float,
+    z0: float,
+    z1: float,
+    deg: int = 10,
+) -> float:
+    """
+    Use a tensor-product Gauss-Legendre quadrature of degree deg
+    over the box [x0, x1] x [y0, y1] x [z0, z1]
+    """
+
+    points, weights = np.polynomial.legendre.leggauss(deg)
+
+    xs = ((x1 - x0) * points + x1 + x0) * 0.5
+    ys = ((y1 - y0) * points + y1 + y0) * 0.5
+    zs = ((z1 - z0) * points + z1 + z0) * 0.5
+    xg, yg, zg = np.meshgrid(xs, ys, zs, indexing="ij")
+    wg = np.outer(np.outer(weights, weights), weights).ravel()
+
+    return (
+        np.sum(f([xg.ravel(), yg.ravel(), zg.ravel()]) * wg)
+        * 0.125
+        * (x1 - x0)
+        * (y1 - y0)
+        * (z1 - z0)
+    )
+
+
 def get_dx(array: NDArray) -> float | None:
     """
     Get the minimum dx from array. Assume that the coordinates are ordered.
@@ -141,13 +196,7 @@ def compute_l1_error_fvm(
 
             error += abs(
                 dx * dy * numerical[i]
-                - simple_quad(
-                    lambda y: simple_quad(
-                        lambda x: analytical([x, y]), x0, x1
-                    ),
-                    y0,
-                    y1,
-                )
+                - simple_quad_2d(analytical, x0, x1, y0, y1)
             )
     else:
         dx = get_dx(x_num[0])
@@ -168,17 +217,7 @@ def compute_l1_error_fvm(
 
             error += abs(
                 dx * dy * dz * numerical[i]
-                - simple_quad(
-                    lambda z: simple_quad(
-                        lambda y: simple_quad(
-                            lambda x: analytical([x, y, z]), x0, x1
-                        ),
-                        y0,
-                        y1,
-                    ),
-                    z0,
-                    z1,
-                )
+                - simple_quad_3d(analytical, x0, x1, y0, y1, z0, z1)
             )
 
     return error