diff --git a/Snakefile b/Snakefile
index e91491fd49c59216591d4b230bdc810524fde3b0..e5ddb4e8195e4e2be82eb7ff57640517697ad9ea 100644
--- a/Snakefile
+++ b/Snakefile
@@ -22,16 +22,20 @@ TODO: Contemplate using lambdify for basis
 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
+
 Urgent:
 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
 TODO: Move height adjustment and stretch factor into implicit solver
     function -> Done
+TODO: Enable choice of equation -> Done
+TODO: Rework burgers_volume_integral()
+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: Enable choice of equation
-TODO: Add subclasses of Burgers (rarefaction and shock case)
 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
diff --git a/config.yaml b/config.yaml
index 89bdfe0ff3de5190021e79e2df5dcd1e69a16381..83eceaa28d71e5c331fe604669f082f08a8fb23a 100644
--- a/config.yaml
+++ b/config.yaml
@@ -1,4 +1,4 @@
-data_dir: 'Oct04'
+data_dir: 'May12'
 random_seed: 1234
 
 # Parameter for Approximation with Troubled Cell Detection
@@ -44,6 +44,9 @@ Approximation:
       detector_config:
           fold_len: 16
       init_cond: 'FourPeakWave'
+    Burgers:
+      equation: 'Burgers'
+      detector: 'Theoretical'
 
 # Parameter for Training Data Generation
 ANN_Data:
diff --git a/scripts/approximate_solution.py b/scripts/approximate_solution.py
index dbf9277289bb43b9ee5945ce8c613128231ddf8c..b7a06df2e25713050dfd659b8f1e390bd75c43ad 100644
--- a/scripts/approximate_solution.py
+++ b/scripts/approximate_solution.py
@@ -6,6 +6,7 @@
 import sys
 import time
 
+from tcd import Equation
 from tcd import Initial_Condition
 from tcd import Limiter
 from tcd import Quadrature
@@ -14,7 +15,7 @@ from tcd import Update_Scheme
 from tcd.Basis_Function import OrthonormalLegendre
 from tcd.DG_Approximation import DGScheme
 from tcd.Mesh import Mesh
-from tcd.Equation import LinearAdvection
+# from tcd.Equation import LinearAdvection
 
 
 def main() -> None:
@@ -84,15 +85,18 @@ def main() -> None:
                              % init_name)
         init_cond = getattr(Initial_Condition, init_name)(config=init_config)
 
-        # Initialize linear advection equation
+        # Initialize equation after checking its legitimacy
         wave_speed = params.pop('wave_speed', 1)
         final_time = params.pop('final_time', 1)
         cfl_number = params.pop('cfl_number', 0.2)
-        equation = LinearAdvection(final_time=final_time,
-                                   cfl_number=cfl_number,
-                                   basis=basis, init_cond=init_cond,
-                                   quadrature=quadrature,
-                                   mesh=mesh, wave_speed=wave_speed)
+
+        equation_name = params.pop('equation', 'LinearAdvection')
+        if not hasattr(Equation, equation_name):
+            raise ValueError('Invalid equation: "%s"' % equation_name)
+        equation = getattr(Equation, equation_name)(
+            final_time=final_time, cfl_number=cfl_number, basis=basis,
+            init_cond=init_cond, quadrature=quadrature, mesh=mesh,
+            wave_speed=wave_speed)
 
         # Initialize update scheme after checking its legitimacy
         scheme_name = params.pop('update_scheme', 'SSPRK3')