diff --git a/src/eegprep/functions/sigprocfunc/topoplot.py b/src/eegprep/functions/sigprocfunc/topoplot.py
index a63aff4c..6ed8003b 100644
--- a/src/eegprep/functions/sigprocfunc/topoplot.py
+++ b/src/eegprep/functions/sigprocfunc/topoplot.py
@@ -45,19 +45,17 @@ 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)
+    # Evaluate at requested points (vectorized)
+    # q is (M, N), xy is (L,) -> d_q is (M, N, L)
+    q = xq + 1j * yq
+    d_q = np.abs(q[:, :, np.newaxis] - xy[np.newaxis, np.newaxis, :])
+
+    with np.errstate(divide='ignore', invalid='ignore'):
+        g_q = (d_q**2) * (np.log(d_q) - 1)  # Green's function
+    g_q[d_q == 0] = 0  # Handle Green's function at zero
+
+    # Weights is (L,), g_q is (M, N, L) -> vq is (M, N)
+    vq = g_q @ weights
 
     return vq