diff --git a/Snakefile b/Snakefile
index fc36903fa8434af112128b923bd116addcf9689f..582a83bbc1248f93655e2cf213aab6b4eb5cc635 100644
--- a/Snakefile
+++ b/Snakefile
@@ -23,7 +23,7 @@ TODO: Contemplate allowing vector input for ICs
 TODO: Discuss how wavelet details should be plotted
 
 TODO: Ask why implicit solver is always over the interval [0,1] with 200 cells
-TODO: Ask why the zero-entry of the boundary matrix is chosen for Burgers
+TODO: Ask about reasoning behind height_adjustment and stretch_factor
 
 Urgent:
 TODO: Mark Burgers' equation as inviscid -> Done
@@ -33,21 +33,29 @@ TODO: Move height adjustment and stretch factor into implicit solver
     function -> Done
 TODO: Enable choice of equation -> Done
 TODO: Ensure termination of Burgers with DiscontinuousConstant function -> Done
-TODO: Ensure termination of Burgers with Sine function
-TODO: Rework burgers_volume_integral()
+TODO: Ensure termination of Burgers with Sine function -> Done
+TODO: Vectorize 'rarefaction_wave()'
+TODO: Vectorize 'burgers_volume_integral()'
+TODO: Vectorize 'burgers_local_Lax_Friedrich()'
+TODO: Vectorize 'burgers_boundary_matrix()'
+TODO: Vectorize '_calculate_boundary_points()'
+TODO: Try to extract calculation of matrices for the volume integral and flux
+    for Burgers
 TODO: Rework/refactor implicit solver
 TODO: Add subclasses of Burgers (rarefaction and shock case)
 TODO: Check correctness of implicit solver (with burgex.f)
 TODO: Add Burgers class completely
+TODO: Add height_adjustment and stretch_factor if sensible
 TODO: Ignore ghost domains by flagging all cells with it for extreme outlier
 TODO: Add functions to create each object from dict
 TODO: Initialize boundary config properly
-TODO: Test Burgers
+TODO: Test Burgers (and fix bugs)
 TODO: Rework Burgers for speed-up and better structure
 TODO: Add initialization to Burgers
 TODO: Use cfl_number for updating, not just time (equation-related?)
 
 Critical, but not urgent:
+TODO: Ensure property variables are used if sensible
 TODO: Make sure TCs are reported as ndarray
 TODO: Make sure that the cell indices are the same over all TCDs
 TODO: Combine ANN workflows if feasible
diff --git a/config.yaml b/config.yaml
index 31af1a68ce03ea81f30736f38beafe7aa29c2d16..4a1d73b96a970b196ffafd8f5288174c6d4673f3 100644
--- a/config.yaml
+++ b/config.yaml
@@ -48,6 +48,10 @@ Approximation:
       equation: 'Burgers'
       detector: 'Theoretical'
       init_cond: 'DiscontinuousConstant'
+    Burgers_Sine:
+      equation: 'Burgers'
+      detector: 'Boxplot'
+      init_cond: 'Sine'
 
 # Parameter for Training Data Generation
 ANN_Data:
diff --git a/scripts/tcd/DG_Approximation.py b/scripts/tcd/DG_Approximation.py
index 02dfab79adecb509b15f5ec09a1a32096fca198f..a42456bf2540d511125ff128650e76c192ca247d 100644
--- a/scripts/tcd/DG_Approximation.py
+++ b/scripts/tcd/DG_Approximation.py
@@ -103,6 +103,7 @@ class DGScheme:
         # Save approximation-specific data in dictionary
         approx_stats = {**self._detector.create_data_dict(projection),
                         **self._equation.create_data_dict(),
+                        'equation': self._equation.name,
                         'time_history': time_history,
                         'troubled_cell_history': troubled_cell_history}
 
diff --git a/scripts/tcd/Equation.py b/scripts/tcd/Equation.py
index 201a03dd519a92f824ce1fd62294e075a7c72138..7817d12b1eb57800b14e41fd5bffb2cf806cf645 100644
--- a/scripts/tcd/Equation.py
+++ b/scripts/tcd/Equation.py
@@ -77,6 +77,11 @@ class Equation(ABC):
 
         self._reset()
 
+    @property
+    def name(self) -> str:
+        """Return string of class name."""
+        return self.__class__.__name__
+
     @property
     def basis(self) -> Basis:
         """Return basis."""
@@ -481,10 +486,10 @@ class Burgers(Equation):
         # return grid, exact
         exact = []
 
-        if self._init_cond.__name__ == 'Sine':
+        if self._init_cond.get_name() == 'Sine':
             # u(x,t) = u(x - u0*time, 0)
             exact = self.implicit_burgers_solver(grid, mesh)
-        elif self._init_cond.__name__ == 'DiscontinuousConstant':
+        elif self._init_cond.get_name() == 'DiscontinuousConstant':
             # u(x,t) = u(x - u0*time, 0)
             exact = self.rarefaction_wave(grid)
 
@@ -520,8 +525,12 @@ class Burgers(Equation):
         Code adapted from DG code by Hesthaven.
 
         """
+        # Temporary fix until further clarification
+        height_adjustment = 0
+        stretch_factor = 1
+
         # Shock speed = 1/2*(left+right values)
-        sspeed = self._init_cond._height_adjustment
+        sspeed = height_adjustment
         uexact = np.zeros(np.size(grid_values))
 
         # initialize values
@@ -575,8 +584,7 @@ class Burgers(Equation):
 
         burgers_exact = uexact.reshape((1, np.size(grid_values)))
 
-        return self._init_cond._height_adjustment + \
-            self._init_cond._stretch_factor * burgers_exact
+        return height_adjustment + stretch_factor * burgers_exact
 
     @enforce_boundary()
     def update_right_hand_side(self, projection: ndarray) -> ndarray:
diff --git a/scripts/tcd/Plotting.py b/scripts/tcd/Plotting.py
index 598065d3026dbaf93e382eb9ffee63a033a14686..c3c623ad273e94a43954f498fd8140d1be4948c5 100644
--- a/scripts/tcd/Plotting.py
+++ b/scripts/tcd/Plotting.py
@@ -17,7 +17,7 @@ from .Quadrature import Quadrature
 from .Initial_Condition import InitialCondition
 from .Basis_Function import Basis, OrthonormalLegendre
 from .projection_utils import calculate_approximate_solution
-from .Equation import Equation, LinearAdvection
+from .Equation import Equation, LinearAdvection, Burgers
 from .Mesh import Mesh
 from .encoding_utils import decode_ndarray
 
@@ -344,11 +344,18 @@ def plot_approximation_results(data_file: str, directory: str, plot_name: str,
     mesh = Mesh(**approx_stats['mesh'])
     print(list(approx_stats.keys()))
 
-    approx_stats['equation'] = LinearAdvection(
-        quadrature=quadrature, init_cond=init_cond, basis=basis, mesh=mesh,
-        final_time=approx_stats['final_time'],
-        wave_speed=approx_stats['wave_speed'],
-        cfl_number=approx_stats['cfl_number'])
+    if approx_stats['equation'] == 'LinearAdvection':
+        approx_stats['equation'] = LinearAdvection(
+            quadrature=quadrature, init_cond=init_cond, basis=basis, mesh=mesh,
+            final_time=approx_stats['final_time'],
+            wave_speed=approx_stats['wave_speed'],
+            cfl_number=approx_stats['cfl_number'])
+    else:
+        approx_stats['equation'] = Burgers(
+            quadrature=quadrature, init_cond=init_cond, basis=basis, mesh=mesh,
+            final_time=approx_stats['final_time'],
+            wave_speed=approx_stats['wave_speed'],
+            cfl_number=approx_stats['cfl_number'])
 
     for key in ['basis', 'mesh', 'final_time', 'wave_speed', 'cfl_number']:
         approx_stats.pop(key, None)