diff --git a/Snakefile b/Snakefile
index 05ca0d46562b11052e129dff7698b96ce6e11be0..1312c34f22d784aa35f58446f8d483e105a8fff2 100644
--- a/Snakefile
+++ b/Snakefile
@@ -24,7 +24,7 @@ TODO: Contemplate using lambdify for basis
 TODO: Discuss how wavelet details should be plotted
 
 Urgent:
-TODO: Vectorize 'plot_details()'
+TODO: Vectorize 'plot_details()' -> Done
 TODO: Replace loops/list comprehension with vectorization if feasible
 TODO: Replace loops with list comprehension if feasible
 TODO: Rework ICs to allow vector input
diff --git a/scripts/tcd/Plotting.py b/scripts/tcd/Plotting.py
index 908292df7a815925c4275dfb74d8ddf032c77a69..14ffcf4de70284df5cb19efb76b0b21b016283cf 100644
--- a/scripts/tcd/Plotting.py
+++ b/scripts/tcd/Plotting.py
@@ -140,26 +140,22 @@ def plot_details(fine_projection: ndarray, fine_mesh: Mesh, basis: Basis,
         Matrix of multiwavelet coefficients.
 
     """
-    num_coarse_mesh_cells = len(coarse_projection[0])
-    averaged_projection = [[coarse_projection[degree][cell]
-                            * basis.basis[degree].subs(x, value)
-                            for cell in range(num_coarse_mesh_cells)
-                            for value in [-0.5, 0.5]]
-                           for degree in range(basis.polynomial_degree + 1)]
-
-    wavelet_projection = [[multiwavelet_coeffs[degree][cell]
-                           * basis.wavelet[degree].subs(z, 0.5) * value
-                           for cell in range(num_coarse_mesh_cells)
-                           for value in [(-1) ** (basis.polynomial_degree
-                                                  + degree + 1), 1]]
-                          for degree in range(basis.polynomial_degree + 1)]
-
-    projected_coarse = np.sum(averaged_projection, axis=0)
-    projected_fine = np.sum([fine_projection[degree]
-                             * basis.basis[degree].subs(x, 0)
-                             for degree in range(basis.polynomial_degree + 1)],
-                            axis=0)
-    projected_wavelet_coeffs = np.sum(wavelet_projection, axis=0)
+    basis_matrix = np.array(
+        [[basis.basis[degree].subs(x, point)
+          for point in [-0.5, 0.5]]
+         for degree in range(basis.polynomial_degree+1)])
+    projected_coarse = (basis_matrix.T @ coarse_projection).flatten('F')
+
+    wavelet_matrix = np.array(
+        [[basis.wavelet[degree].subs(z, 0.5) * value
+          for value in [(-1) ** (basis.polynomial_degree + degree + 1), 1]]
+         for degree in range(basis.polynomial_degree + 1)])
+    projected_wavelet_coeffs = (wavelet_matrix.T @
+                                multiwavelet_coeffs[:, ::2]).flatten('F')
+
+    basis_vector = [basis.basis[degree].subs(x, 0)
+                    for degree in range(basis.polynomial_degree + 1)]
+    projected_fine = basis_vector @ fine_projection
 
     plt.figure('coeff_details')
     plt.plot(fine_mesh.non_ghost_cells, projected_fine-projected_coarse, 'm-.')