diff --git a/Snakefile b/Snakefile
index 8ae8e4a649121c292a5552890cc4c62f88ccf10a..b25379db932a97033958d8462e79f02e674ee118 100644
--- a/Snakefile
+++ b/Snakefile
@@ -21,14 +21,21 @@ TODO: Discuss name for quadrature mesh (now: grid)
 TODO: Contemplate using lambdify for basis
 TODO: Contemplate allowing vector input for ICs
 TODO: Discuss how wavelet details should be plotted
+
 TODO: Ask whether we are dealing with inviscid Burgers only
 TODO: Ask why we limit the initial time step
 TODO: Ask whether cfl number should be absolute value
 TODO: Ask why cells 1 and -2 are adjusted for Dirichlet boundary
+TODO: Ask why height adjustment and stretch factor are considered outside of
+    implicit solver function
+TODO: Ask about workings of implicit solver for Burgers
+TODO: Ask whether Burgers has compatible ICs other than Sine and
+    DiscontinuousConstant
 
 Urgent:
 TODO: Add Burgers class completely
 TODO: Clean up exact solution for Burgers
+TODO: Adapt code to accommodate boundary decorator -> Done
 TODO: Add subclasses of Burgers
 TODO: Ignore ghost domains by flagging all cells with it for extreme outlier
 TODO: Add functions to create each object from dict
diff --git a/scripts/approximate_solution.py b/scripts/approximate_solution.py
index ec3669045939471c3bc23d1a6d3cb97c2c43af39..dbf9277289bb43b9ee5945ce8c613128231ddf8c 100644
--- a/scripts/approximate_solution.py
+++ b/scripts/approximate_solution.py
@@ -45,7 +45,8 @@ def main() -> None:
         mesh = Mesh(num_cells=params.pop('num_mesh_cells', 64),
                     left_bound=params.pop('left_bound', -1),
                     right_bound=params.pop('right_bound', 1),
-                    num_ghost_cells=num_ghost_cells)
+                    num_ghost_cells=num_ghost_cells,
+                    boundary_config={'type': 'periodic'})
 
         # Initialize basis
         basis = OrthonormalLegendre(params.pop('polynomial_degree', 2))
diff --git a/scripts/tcd/ANN_Data_Generator.py b/scripts/tcd/ANN_Data_Generator.py
index 1aef27ce705284c194bac112fd2e20282b89e9c2..a3900861f86d8fe9ca6d1cd817080ef7e34d15e7 100644
--- a/scripts/tcd/ANN_Data_Generator.py
+++ b/scripts/tcd/ANN_Data_Generator.py
@@ -40,7 +40,8 @@ class TrainingDataGenerator:
         self._quadrature_list = [Gauss({'num_nodes': pol_deg+1})
                                  for pol_deg in range(7)]
         self._mesh_list = [Mesh(left_bound=-1, right_bound=1,
-                                num_cells=2**exp, num_ghost_cells=0)
+                                num_cells=2**exp, num_ghost_cells=0,
+                                boundary_config={'type': 'periodic'})
                            for exp in range(5, 12)]
 
     def build_training_data(self, init_cond_list, num_samples, balance=0.5,
diff --git a/scripts/tcd/Boundary_Condition.py b/scripts/tcd/Boundary_Condition.py
index 1cde4ea28cb6c6f839a0564c6da446c11bdc1348..cb9f61082c50295e4cc4f5f52ae53aa616cc76d6 100644
--- a/scripts/tcd/Boundary_Condition.py
+++ b/scripts/tcd/Boundary_Condition.py
@@ -88,7 +88,7 @@ def enforce_boundary(num_ghost_cells=None):
                     projection=projection, config=config)
             elif self._mesh.boundary_config['type'] == 'dirichlet':
                 return dirichlet_boundary(projection=projection,
-                                          )
+                                          config=self._mesh.boundary_config)
             else:
                 raise Exception('Not implemented!')
         return boundary
@@ -98,7 +98,7 @@ def enforce_boundary(num_ghost_cells=None):
 def enforce_boxplot_boundaries(func):
     def boxplot_boundary(self, *args, **kwargs):
         bounds = np.array(func(self, *args, **kwargs))
-        if self._mesh.boundary['type'] == 'periodic':
+        if self._mesh.boundary_config['type'] == 'periodic':
             return tuple(periodic_boundary(projection=bounds,
                                            config=self._mesh.boundary_config))
         else:
@@ -109,7 +109,7 @@ def enforce_boxplot_boundaries(func):
 def enforce_fold_boundaries(func):
     def fold_boundaries(self, *args, **kwargs):
         folds = np.array(func(self, *args, **kwargs))
-        if self._mesh.boundary['type'] == 'periodic':
+        if self._mesh.boundary_config['type'] == 'periodic':
             return folds % self._mesh.num_cells
         else:
             raise Exception('Not implemented!')
diff --git a/scripts/tcd/DG_Approximation.py b/scripts/tcd/DG_Approximation.py
index d3e044b4d84e9e912ad2d77dde7106e167584915..02dfab79adecb509b15f5ec09a1a32096fca198f 100644
--- a/scripts/tcd/DG_Approximation.py
+++ b/scripts/tcd/DG_Approximation.py
@@ -86,7 +86,8 @@ class DGScheme:
         time_history = []
         while current_time < self._equation.final_time:
             cfl_number, time_step = self._equation.update_time_step(
-                current_time=current_time, time_step=time_step)
+                projection=projection, current_time=current_time,
+                time_step=time_step)
 
             # Update projection
             projection, troubled_cells = self._update_scheme.step(projection,
diff --git a/scripts/tcd/Mesh.py b/scripts/tcd/Mesh.py
index 2129342547296368b9a404da523ab27084d8e5dc..866d03a482a70c04174787040893615b9cd95f3e 100644
--- a/scripts/tcd/Mesh.py
+++ b/scripts/tcd/Mesh.py
@@ -84,8 +84,10 @@ class Mesh:
         self._left_bound = left_bound
         self._right_bound = right_bound
 
+        boundary_config['num_ghost_cells'] = num_ghost_cells
         self._boundary_config = boundary_config
         boundary_type = self._boundary_config.get('type', 'periodic')
+
         if boundary_type not in BOUNDARY_TYPES:
             raise ValueError(f'The boundary type "{boundary_type}" '
                              f'is not implemented')
@@ -170,4 +172,5 @@ class Mesh:
         # Return new mesh instance
         return Mesh(left_bound=point - stencil_len/2 * mesh_spacing,
                     right_bound=point + stencil_len/2 * mesh_spacing,
-                    num_cells=stencil_len, num_ghost_cells=0)
+                    num_cells=stencil_len, num_ghost_cells=0,
+                    boundary_config=self.boundary_config)
diff --git a/scripts/tcd/Plotting.py b/scripts/tcd/Plotting.py
index 84891296b2d7eabf4a35f4150647e7b86f8a77c0..598065d3026dbaf93e382eb9ffee63a033a14686 100644
--- a/scripts/tcd/Plotting.py
+++ b/scripts/tcd/Plotting.py
@@ -428,7 +428,8 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
         coarse_mesh = Mesh(num_cells=mesh.num_cells//2,
                            left_bound=mesh.bounds[0],
                            right_bound=mesh.bounds[1],
-                           num_ghost_cells=0)
+                           num_ghost_cells=0,
+                           boundary_config=mesh.boundary_config)
 
         # Plot exact and approximate solutions for coarse mesh
         coarse_grid, coarse_exact = equation.solve_exactly(coarse_mesh)