diff --git a/.jules/bolt.md b/.jules/bolt.md
new file mode 100644
index 00000000..425ac507
--- /dev/null
+++ b/.jules/bolt.md
@@ -0,0 +1,3 @@
+## 2025-03-05 - [Grid evaluation vectorization in topoplot]
+**Learning:** Nested Python loops over coordinate meshes (`xq`, `yq`) are extremely slow when executing NumPy calculations inside, especially for mathematical formulas like biharmonic spline interpolation. Evaluating elements individually causes substantial overhead. Vectorizing coordinate meshes using complex numbers and broadcasting enables fast operations at scale.
+**Action:** Always refactor meshgrid-based spatial evaluations to use complex 1D/2D broadcasting and matrix multiplication with `@`.
diff --git a/src/eegprep/functions/sigprocfunc/topoplot.py b/src/eegprep/functions/sigprocfunc/topoplot.py
index f6bdb838..e973f313 100644
--- a/src/eegprep/functions/sigprocfunc/topoplot.py
+++ b/src/eegprep/functions/sigprocfunc/topoplot.py
@@ -45,19 +45,21 @@ def griddata_v4(x, y, v, xq, yq):
         # If still singular, use pseudoinverse as last resort
         weights = np.linalg.pinv(g_reg) @ v
 
-    # Initialize output array
-    m, n = xq.shape
-    vq = np.zeros_like(xq)
-
-    # Evaluate at requested points
-    xy = xy[:, None]  # Make it column vector for broadcasting
-    for i in range(m):
-        for j in range(n):
-            d = np.abs(xq[i, j] + 1j * yq[i, j] - xy.ravel())
-            with np.errstate(divide='ignore', invalid='ignore'):
-                g = (d**2) * (np.log(d) - 1)  # Green's function
-            g[d == 0] = 0  # Handle Green's function at zero
-            vq[i, j] = np.dot(g, weights)
+    # Vectorized evaluation at query points:
+    # 1. Flatten the coordinates to form a 1D complex array of query locations
+    xy_query_flat = (xq + 1j * yq).ravel()
+    xy_flat = xy.ravel()
+
+    # 2. Compute pairwise distances using broadcasting
+    d_q = np.abs(xy_query_flat[:, None] - xy_flat[None, :])
+
+    # 3. Compute Green's function values
+    with np.errstate(divide='ignore', invalid='ignore'):
+        g_q = (d_q**2) * (np.log(d_q) - 1)
+    g_q[d_q == 0] = 0
+
+    # 4. Perform matrix multiplication and reshape back to original grid shape
+    vq = (g_q @ weights).reshape(xq.shape)
 
     return vq