diff --git a/Basis_Function.py b/Basis_Function.py
index 7efd115cb0970fa9969248ab3fdcc499495909d7..63da57912fabd0f19c55fd66c9cd743ae8648697 100644
--- a/Basis_Function.py
+++ b/Basis_Function.py
@@ -3,11 +3,12 @@
 @author: Laura C. Kühle
 
 TODO: Contemplate whether calculating projections during initialization can
-    save time
+    save time -> Done (solved due to property attributes)
 
 """
 import numpy as np
 from sympy import Symbol, integrate
+from functools import cache
 
 from projection_utils import calculate_approximate_solution
 
@@ -16,7 +17,7 @@ z = Symbol('z')
 
 
 class Basis:
-    """Class for basis vector.
+    """Class for polynomial basis.
 
     Attributes
     ----------
@@ -26,19 +27,11 @@ class Basis:
         Array of wavelet.
     inv_mass: ndarray
         Inverse mass matrix.
+    basis_projection: ndarray
+
 
     Methods
     -------
-    get_basis_vector()
-        Returns basis vector.
-    get_wavelet_vector
-        Returns wavelet vector.
-    get_inverse_mass_matrix()
-        Returns inverse mass matrix.
-    get_basis_projections()
-        Returns basis projections.
-    get_wavelet_projections()
-        Returns wavelet projections.
     calculate_cell_average(projection, stencil_length, add_reconstructions)
         Calculate cell averages for a given projection.
 
@@ -53,13 +46,12 @@ class Basis:
 
         """
         self._polynomial_degree = polynomial_degree
-        self._basis = self._build_basis_vector(x)
-        self._wavelet = self._build_wavelet_vector(z)
-        self._inv_mass = self._build_inverse_mass_matrix()
 
-    def get_basis_vector(self):
-        """Returns basis vector."""
-        return self._basis
+    @property
+    @cache
+    def basis(self):
+        """Return basis vector."""
+        return self._build_basis_vector(x)
 
     def _build_basis_vector(self, eval_point):
         """Constructs basis vector.
@@ -77,9 +69,11 @@ class Basis:
         """
         return []
 
-    def get_wavelet_vector(self):
-        """Returns wavelet vector."""
-        return self._wavelet
+    @property
+    @cache
+    def wavelet(self):
+        """Return basis vector."""
+        return self._build_wavelet_vector(z)
 
     def _build_wavelet_vector(self, eval_point):
         """Constructs wavelet vector.
@@ -97,9 +91,11 @@ class Basis:
         """
         return []
 
-    def get_inverse_mass_matrix(self):
-        """Returns inverse mass matrix."""
-        return self._inv_mass
+    @property
+    @cache
+    def inverse_mass_matrix(self):
+        """Return inverse mass matrix."""
+        return self._build_inverse_mass_matrix()
 
     def _build_inverse_mass_matrix(self):
         """Constructs inverse mass matrix.
@@ -112,13 +108,17 @@ class Basis:
         """
         pass
 
-    def get_basis_projections(self):
-        """Returns basis projection."""
-        pass
+    @property
+    @cache
+    def basis_projection(self):
+        """Return basis projection."""
+        return []
 
-    def get_multiwavelet_projections(self):
-        """Returns wavelet projection."""
-        pass
+    @property
+    @cache
+    def wavelet_projection(self):
+        """Return wavelet projection."""
+        return []
 
     def calculate_cell_average(self, projection, stencil_length,
                                add_reconstructions=True):
@@ -146,7 +146,7 @@ class Basis:
 
         """
         cell_averages = calculate_approximate_solution(
-            projection, np.array([0]), 0, self._basis)
+            projection, np.array([0]), 0, self.basis)
 
         if add_reconstructions:
             middle_idx = stencil_length // 2
@@ -178,9 +178,11 @@ class Basis:
 
         """
         left_reconstructions = calculate_approximate_solution(
-            projection, np.array([-1]), self._polynomial_degree, self._basis)
+            projection, np.array([-1]), self._polynomial_degree,
+            self.basis)
         right_reconstructions = calculate_approximate_solution(
-            projection, np.array([1]), self._polynomial_degree, self._basis)
+            projection, np.array([1]), self._polynomial_degree,
+            self.basis)
         return left_reconstructions, right_reconstructions
 
 
@@ -232,16 +234,7 @@ class Legendre(Basis):
 
 
 class OrthonormalLegendre(Legendre):
-    """Class for orthonormal Legendre basis.
-
-    Methods
-    -------
-    get_basis_projection()
-        Returns basis projection.
-    get_wavelet_projection()
-        Returns wavelet projection.
-
-    """
+    """Class for orthonormal Legendre basis."""
     def _build_basis_vector(self, eval_point):
         """Constructs basis vector.
 
@@ -337,7 +330,9 @@ class OrthonormalLegendre(Legendre):
             mass_matrix.append(new_row)
         return np.array(mass_matrix)
 
-    def get_basis_projections(self):
+    @property
+    @cache
+    def basis_projections(self):
         """Returns basis projection.
 
         Returns
@@ -371,14 +366,16 @@ class OrthonormalLegendre(Legendre):
         for i in range(self._polynomial_degree + 1):
             row = []
             for j in range(self._polynomial_degree + 1):
-                entry = integrate(self._basis[i].subs(x, first_param)
-                                  * self._basis[j].subs(x, second_param),
+                entry = integrate(self.basis[i].subs(x, first_param)
+                                  * self.basis[j].subs(x, second_param),
                                   (z, -1, 1))
                 row.append(np.float64(entry))
             matrix.append(row)
         return matrix
 
-    def get_multiwavelet_projections(self):
+    @property
+    @cache
+    def multiwavelet_projections(self):
         """Returns wavelet projection.
 
         Returns
@@ -419,8 +416,8 @@ class OrthonormalLegendre(Legendre):
         for i in range(self._polynomial_degree+1):
             row = []
             for j in range(self._polynomial_degree+1):
-                entry = integrate(self._basis[i].subs(x, first_param)
-                                  * self._wavelet[j].subs(z, second_param),
+                entry = integrate(self.basis[i].subs(x, first_param)
+                                  * self.wavelet[j].subs(z, second_param),
                                   (z, -1, 1))
                 if is_left_matrix:
                     entry = entry * (-1)**(j + self._polynomial_degree + 1)
diff --git a/DG_Approximation.py b/DG_Approximation.py
index 11a25c55d9d57cfd4f3785961e80d55b37ec3fdf..ef8931b2e3d75a39528852bd0cb671cd16eb247f 100644
--- a/DG_Approximation.py
+++ b/DG_Approximation.py
@@ -6,13 +6,16 @@ Discussion:
 TODO: Contemplate saving 5-CV split and evaluating models separately
 TODO: Contemplate separating cell average and reconstruction calculations
     completely
+TODO: Contemplate removing Methods section from class docstring
 
 Urgent:
 TODO: Investigate protected vs. public variables -> Done
-TODO: Investigate decorators and properties
-TODO: Investigate how to make variables read-only
-TODO: Make basis variables public (if feasible)
+TODO: Investigate decorators and properties -> Done
+TODO: Investigate how to make variables read-only -> Done
+TODO: Replace getter with property attributes for basis -> Done
 TODO: Contain polynomial degree in basis
+TODO: Improve docstring for Basis class
+TODO: Enforce abstract classes/methods (abc.ABC, abc.abstractmethod)
 TODO: Introduce Mesh class
     (mesh, cell_len, num_grid_cells, bounds, num_ghost_cells, etc.)
 TODO: Check whether ghost cells are handled/set correctly
@@ -29,7 +32,6 @@ Critical, but not urgent:
 TODO: Investigate g-mesh(?)
 TODO: Force input_size for each ANN model to be stencil length
 TODO: Use full path for ANN model state
-TODO: Enforce abstract classes/methods (abc.ABC, abc.abstractmethod)
 TODO: Extract object initialization from DGScheme
 TODO: Use cfl_number for updating, not just time
 
@@ -340,7 +342,7 @@ def do_initial_projection(initial_condition, basis, quadrature,
     ----------
     initial_condition : InitialCondition object
         Initial condition used for calculation.
-    basis: Vector object
+    basis: Basis object
         Basis used for calculation.
     quadrature: Quadrature object
         Quadrature used for evaluation.
@@ -365,7 +367,7 @@ def do_initial_projection(initial_condition, basis, quadrature,
     """
     # Initialize matrix and set first entry to accommodate for ghost cell
     output_matrix = [0]
-    basis_vector = basis.get_basis_vector()
+    basis_vector = basis.basis
 
     cell_len = (right_bound-left_bound)/num_grid_cells
     for cell in range(num_grid_cells):
@@ -382,7 +384,7 @@ def do_initial_projection(initial_condition, basis, quadrature,
                 for point in range(quadrature.get_num_points()))))
 
         new_row = np.array(new_row)
-        output_matrix.append(basis.get_inverse_mass_matrix() @ new_row)
+        output_matrix.append(basis.inverse_mass_matrix @ new_row)
 
     # Set ghost cells to respective value
     output_matrix[0] = output_matrix[num_grid_cells]
diff --git a/Plotting.py b/Plotting.py
index dd6a58205ad9bcb26323061956145d7f143ee658..cfff3bace7357637dea9057be64bd765efe02b6b 100644
--- a/Plotting.py
+++ b/Plotting.py
@@ -337,7 +337,7 @@ def plot_approximation_results(data_file: str, directory: str, plot_name: str,
         Path to directory in which plots will be saved.
     plot_name : str
         Name of plot.
-    basis: Vector object
+    basis: Basis object
         Basis used for calculation.
     quadrature: Quadrature object
         Quadrature used for evaluation.
@@ -444,7 +444,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
         init_cond)
     approx = calculate_approximate_solution(
         projection[:, 1:-1], quadrature.get_eval_points(),
-        polynomial_degree, basis.get_basis_vector())
+        polynomial_degree, basis.basis)
 
     # Plot multiwavelet solution (fine and coarse grid)
     if coarse_projection is not None:
@@ -460,7 +460,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
             init_cond)
         coarse_approx = calculate_approximate_solution(
             coarse_projection, quadrature.get_eval_points(),
-            polynomial_degree, basis.get_basis_vector())
+            polynomial_degree, basis.basis)
         plot_solution_and_approx(
             coarse_grid, coarse_exact, coarse_approx, colors['coarse_exact'],
             colors['coarse_approx'])
@@ -468,8 +468,7 @@ def plot_results(projection: ndarray, troubled_cell_history: list,
         # Plot multiwavelet details
         num_coarse_grid_cells = num_grid_cells//2
         plot_details(projection[:, 1:-1], mesh[2:-2], coarse_projection,
-                     basis.get_basis_vector(),
-                     basis.get_wavelet_vector(), multiwavelet_coeffs,
+                     basis.basis, basis.wavelet, multiwavelet_coeffs,
                      num_coarse_grid_cells,
                      polynomial_degree)
 
diff --git a/Troubled_Cell_Detector.py b/Troubled_Cell_Detector.py
index 00447cf219e0c3906973564fbc9159c83803c84a..6db101be3b8ef1d7f5e7f9d10c96b3a3767bd4a6 100644
--- a/Troubled_Cell_Detector.py
+++ b/Troubled_Cell_Detector.py
@@ -250,7 +250,7 @@ class WaveletDetector(TroubledCellDetector):
         # Set additional necessary parameter
         self._num_coarse_grid_cells = self._num_grid_cells//2
         self._wavelet_projection_left, self._wavelet_projection_right \
-            = self._basis.get_multiwavelet_projections()
+            = self._basis.multiwavelet_projections
 
     def get_cells(self, projection):
         """Calculates troubled cells in a given projection.
@@ -325,7 +325,7 @@ class WaveletDetector(TroubledCellDetector):
 
         """
         basis_projection_left, basis_projection_right\
-            = self._basis.get_basis_projections()
+            = self._basis.basis_projections
 
         # Remove ghost cells
         projection = projection[:, 1:-1]