diff --git a/Snakefile b/Snakefile
index 7a11ef2edfba3a894d60f8d13a57715f0ffe18ac..6e791853db3a62120e66b9b48b22bd067cfc3357 100644
--- a/Snakefile
+++ b/Snakefile
@@ -23,7 +23,7 @@ TODO: Contemplate allowing vector input for ICs
 TODO: Discuss how wavelet details should be plotted
 
 Urgent:
-TODO: Rename stiffness matrix to volume integral matrix
+TODO: Rename stiffness matrix to volume integral matrix -> Done
 TODO: Rename boundary matrix to flux matrix
 TODO: Change num_ghost_cells to be 1 for calc and 0 for other
 TODO: Move boundary condition to Mesh class
diff --git a/scripts/tcd/Equation.py b/scripts/tcd/Equation.py
index e0e7c21cfcc19befeb0ee50dff89d14b70908485..afe266b7763de5449eef8a3ba3dcdd071e5995d7 100644
--- a/scripts/tcd/Equation.py
+++ b/scripts/tcd/Equation.py
@@ -179,8 +179,8 @@ class LinearAdvection(Equation):
 
     Attributes
     ----------
-    stiffness_matrix : ndarray
-        Stiffness matrix.
+    volume_integral_matrix : ndarray
+        Volume integral matrix.
     boundary_matrix : ndarray
         Boundary matrix.
 
@@ -206,7 +206,7 @@ class LinearAdvection(Equation):
         degree_matrix = np.matmul(root_vector[:, np.newaxis],
                                   root_vector[:, np.newaxis].T)
 
-        # Set stiffness matrix
+        # Set volume integral matrix
         matrix = np.ones(matrix_shape)
         if self._wave_speed > 0:
             matrix[np.fromfunction(
@@ -214,7 +214,7 @@ class LinearAdvection(Equation):
         else:
             matrix[np.fromfunction(
                 lambda i, j: (j > i) & ((i+j) % 2 == 1), matrix_shape)] = -1.0
-        self._stiffness_matrix = matrix * degree_matrix
+        self._volume_integral_matrix = matrix * degree_matrix
 
         # Set boundary matrix
         matrix = np.fromfunction(lambda i, j: (-1.0)**i, matrix_shape) \
@@ -326,7 +326,7 @@ class LinearAdvection(Equation):
                 2 * (self._boundary_matrix @
                      projection[:, self._mesh.num_ghost_cells-1:
                                 -self._mesh.num_ghost_cells-1] +
-                     self._stiffness_matrix @
+                     self._volume_integral_matrix @
                      projection[:, self._mesh.num_ghost_cells:
                                 -self._mesh.num_ghost_cells])
         else:
@@ -334,7 +334,7 @@ class LinearAdvection(Equation):
                             -self._mesh.num_ghost_cells] = \
                 2 * (self._boundary_matrix @
                      projection[:, self._mesh.num_ghost_cells+1:] +
-                     self._stiffness_matrix @
+                     self._volume_integral_matrix @
                      projection[:, self._mesh.num_ghost_cells:
                                 -self._mesh.num_ghost_cells])