diff --git a/Snakefile b/Snakefile
index ea3fb9b83c650d2bb010cde3e0976ee6a6e090f7..1b5ddc0c2ec70c050a68850623cdf60c21774bf6 100644
--- a/Snakefile
+++ b/Snakefile
@@ -26,6 +26,7 @@ TODO: Ask why implicit solver is always over the interval [0,1] with 200 cells
 TODO: Ask about reasoning behind height_adjustment and stretch_factor
 
 Urgent:
+TODO: Move code to other Git
 TODO: Mark Burgers' equation as inviscid -> Done
 TODO: Make sure CFL number is absolute value -> Done
 TODO: Make sure only ghost cells are limited for Dirichlet boundary -> Done
@@ -39,6 +40,7 @@ TODO: Vectorize 'burgers_volume_integral()'
 TODO: Vectorize 'burgers_local_Lax_Friedrich()'
 TODO: Vectorize 'burgers_boundary_matrix()'
 TODO: Vectorize '_calculate_boundary_points()'
+TODO: Vectorize 'update_right_hand_side()' for Burgers -> Done
 TODO: Try to extract calculation of matrices for the volume integral and flux
     for Burgers
 TODO: Rework/refactor implicit solver
diff --git a/scripts/tcd/Equation.py b/scripts/tcd/Equation.py
index 783dd22e8bad4c1b05c6dee9cbec5978f89461e4..c518c91ad650c91584f773a182e97689112071bc 100644
--- a/scripts/tcd/Equation.py
+++ b/scripts/tcd/Equation.py
@@ -614,44 +614,9 @@ class Burgers(Equation):
             Matrix of right-hand side.
 
         """
-        pass
         volume_integral_matrix = self._burgers_volume_integral(projection)
         boundary_matrix = self._burgers_boundary_matrix(projection)
-
-        # Initialize vector and set first entry to accommodate for ghost cell
-        right_hand_side = [0]
-
-        for j in range(self._mesh.num_cells):
-            right_hand_side.append(2*(volume_integral_matrix[:, j + 1] +
-                                      boundary_matrix[:, j + 1]))
-
-        # Set ghost cells to respective value
-        # (Periodic, Updated to Dirichlet in enforce_boundary_condition
-        # for DiscontinuousConstant problem)
-        right_hand_side[0] = right_hand_side[self._mesh.num_cells]
-        right_hand_side.append(right_hand_side[1])
-
-        return np.transpose(right_hand_side)
-        # right_hand_side = np.zeros_like(projection)
-        # if self._wave_speed > 0:
-        #     right_hand_side[:, self._mesh.num_ghost_cells:
-        #                     -self._mesh.num_ghost_cells] = \
-        #         2 * (self._flux_matrix @
-        #              projection[:, self._mesh.num_ghost_cells-1:
-        #                         -self._mesh.num_ghost_cells-1] +
-        #              self._volume_integral_matrix @
-        #              projection[:, self._mesh.num_ghost_cells:
-        #                         -self._mesh.num_ghost_cells])
-        # else:
-        #     right_hand_side[:, self._mesh.num_ghost_cells:
-        #                     -self._mesh.num_ghost_cells] = \
-        #         2 * (self._flux_matrix @
-        #              projection[:, self._mesh.num_ghost_cells+1:] +
-        #              self._volume_integral_matrix @
-        #              projection[:, self._mesh.num_ghost_cells:
-        #                         -self._mesh.num_ghost_cells])
-        #
-        # return right_hand_side
+        return 2*(volume_integral_matrix+boundary_matrix)
 
     def _burgers_volume_integral(self, projection):
         basis_vector = self._basis.basis