diff --git a/scripts/tcd/Basis_Function.py b/scripts/tcd/Basis_Function.py
index 35701b09241055fa7e22e2d2149bcf44503df374..ef334a9f03d5cf178e591377214624fd41121bee 100644
--- a/scripts/tcd/Basis_Function.py
+++ b/scripts/tcd/Basis_Function.py
@@ -403,16 +403,15 @@ class OrthonormalLegendre(Legendre):
             Matrix containing the integral of basis products.
 
         """
-        matrix = []
-        for i in range(self._polynomial_degree + 1):
-            row = []
-            for j in range(self._polynomial_degree + 1):
-                entry = integrate(self.basis[i].subs(x, first_param)
-                                  * self.basis[j].subs(x, second_param),
-                                  (z, -1, 1))
-                row.append(np.float64(entry))
-            matrix.append(row)
-        return np.array(matrix)
+        basis_row = np.array([self.basis[idx].subs(x, first_param) for idx in
+                              range(self._polynomial_degree+1)])[: np.newaxis]
+        basis_col = np.array([self.basis[idx].subs(x, second_param) for idx in
+                              range(self._polynomial_degree+1)])[: np.newaxis]
+        basis_matrix = np.matmul(basis_row[:, np.newaxis],
+                                 basis_col[:, np.newaxis].T)
+        basis_matrix = np.float64(np.vectorize(
+            lambda y: integrate(y, (z, -1, 1)))(basis_matrix))
+        return basis_matrix
 
     @property
     @cache