diff --git a/Basis_Function.py b/Basis_Function.py
index 14315ad7341c96d5e5379955eaf542dc98972b76..5e69d31b0d896d1df9c6c612200b6b85fae3e363 100644
--- a/Basis_Function.py
+++ b/Basis_Function.py
@@ -5,7 +5,9 @@
 
 """
 from functools import cache
+from typing import Tuple
 import numpy as np
+from numpy import ndarray
 from sympy import Symbol, integrate
 
 from projection_utils import calculate_approximate_solution
@@ -44,7 +46,7 @@ class Basis:
 
     """
 
-    def __init__(self, polynomial_degree):
+    def __init__(self, polynomial_degree: int) -> None:
         """Initialize Basis.
 
         Parameters
@@ -56,17 +58,17 @@ class Basis:
         self._polynomial_degree = polynomial_degree
 
     @property
-    def polynomial_degree(self):
+    def polynomial_degree(self) -> int:
         """Return polynomial degree."""
         return self._polynomial_degree
 
     @property
     @cache
-    def basis(self):
+    def basis(self) -> ndarray:
         """Return basis vector."""
         return self._build_basis_vector(x)
 
-    def _build_basis_vector(self, eval_point):
+    def _build_basis_vector(self, eval_point: float) -> ndarray:
         """Construct basis vector.
 
         Parameters
@@ -80,15 +82,15 @@ class Basis:
             Vector containing basis evaluated at evaluation point.
 
         """
-        return []
+        pass
 
     @property
     @cache
-    def wavelet(self):
+    def wavelet(self) -> ndarray:
         """Return wavelet vector."""
         return self._build_wavelet_vector(z)
 
-    def _build_wavelet_vector(self, eval_point):
+    def _build_wavelet_vector(self, eval_point: float) -> ndarray:
         """Construct wavelet vector.
 
         Parameters
@@ -102,15 +104,15 @@ class Basis:
             Vector containing wavelet evaluated at evaluation point.
 
         """
-        return []
+        pass
 
     @property
     @cache
-    def inverse_mass_matrix(self):
+    def inverse_mass_matrix(self) -> ndarray:
         """Return inverse mass matrix."""
         return self._build_inverse_mass_matrix()
 
-    def _build_inverse_mass_matrix(self):
+    def _build_inverse_mass_matrix(self) -> ndarray:
         """Construct inverse mass matrix.
 
         Returns
@@ -123,18 +125,18 @@ class Basis:
 
     @property
     @cache
-    def basis_projection(self):
+    def basis_projection(self) -> Tuple[ndarray, ndarray]:
         """Return basis projection."""
-        return []
+        return np.array([]), np.array([])
 
     @property
     @cache
-    def wavelet_projection(self):
+    def wavelet_projection(self) -> Tuple[ndarray, ndarray]:
         """Return wavelet projection."""
-        return []
+        return np.array([]), np.array([])
 
-    def calculate_cell_average(self, projection, stencil_length,
-                               add_reconstructions=True):
+    def calculate_cell_average(self, projection: ndarray, stencil_length: int,
+                               add_reconstructions: bool = True) -> ndarray:
         """Calculate cell averages for a given projection.
 
         Calculate the cell averages of all cells in a projection.
@@ -174,7 +176,8 @@ class Basis:
                                 cell_averages[:, middle_idx+1:]))))
         return np.array(list(map(np.float64, cell_averages)))
 
-    def _calculate_reconstructions(self, projection):
+    def _calculate_reconstructions(self, projection: ndarray) \
+            -> Tuple[ndarray, ndarray]:
         """Calculate left and right reconstructions for a given projection.
 
         Parameters
@@ -202,7 +205,7 @@ class Basis:
 class Legendre(Basis):
     """Class for Legendre basis."""
 
-    def _build_basis_vector(self, eval_point):
+    def _build_basis_vector(self, eval_point: float) -> ndarray:
         """Construct basis vector.
 
         Parameters
@@ -218,7 +221,7 @@ class Legendre(Basis):
         """
         return self._calculate_legendre_vector(eval_point)
 
-    def _calculate_legendre_vector(self, eval_point):
+    def _calculate_legendre_vector(self, eval_point: float) -> ndarray:
         """Construct Legendre vector.
 
         Parameters
@@ -244,13 +247,13 @@ class Legendre(Basis):
                     poly = (2.0*degree - 1)/degree * eval_point * vector[-1] \
                            - (degree-1)/degree * vector[-2]
                     vector.append(poly)
-        return vector
+        return np.array(vector)
 
 
 class OrthonormalLegendre(Legendre):
     """Class for orthonormal Legendre basis."""
 
-    def _build_basis_vector(self, eval_point):
+    def _build_basis_vector(self, eval_point: float) -> ndarray:
         """Construct basis vector.
 
         Parameters
@@ -265,10 +268,10 @@ class OrthonormalLegendre(Legendre):
 
         """
         leg_vector = self._calculate_legendre_vector(eval_point)
-        return [leg_vector[degree] * np.sqrt(degree+0.5)
-                for degree in range(self._polynomial_degree+1)]
+        return np.array([leg_vector[degree] * np.sqrt(degree+0.5)
+                         for degree in range(self._polynomial_degree+1)])
 
-    def _build_wavelet_vector(self, eval_point):
+    def _build_wavelet_vector(self, eval_point: float) -> ndarray:
         """Construct wavelet vector.
 
         Parameters
@@ -289,51 +292,51 @@ class OrthonormalLegendre(Legendre):
         degree = self._polynomial_degree
 
         if degree == 0:
-            return [np.sqrt(0.5) + eval_point*0]
+            return np.array([np.sqrt(0.5) + eval_point*0])
         if degree == 1:
-            return [np.sqrt(1.5) * (-1 + 2*eval_point),
-                    np.sqrt(0.5) * (-2 + 3*eval_point)]
+            return np.array([np.sqrt(1.5) * (-1 + 2*eval_point),
+                             np.sqrt(0.5) * (-2 + 3*eval_point)])
         if degree == 2:
-            return [1/3 * np.sqrt(0.5) *
-                    (1 - 24*eval_point + 30*(eval_point**2)),
-                    1/2 * np.sqrt(1.5) *
-                    (3 - 16*eval_point + 15*(eval_point**2)),
-                    1/3 * np.sqrt(2.5) *
-                    (4 - 15*eval_point + 12*(eval_point**2))]
+            return np.array([1/3 * np.sqrt(0.5) *
+                             (1 - 24*eval_point + 30*(eval_point**2)),
+                             1/2 * np.sqrt(1.5) *
+                             (3 - 16*eval_point + 15*(eval_point**2)),
+                             1/3 * np.sqrt(2.5) *
+                             (4 - 15*eval_point + 12*(eval_point**2))])
         if degree == 3:
-            return [np.sqrt(15/34) *
-                    (1 + 4*eval_point - 30*(eval_point**2)
-                     + 28*(eval_point**3)),
-                    np.sqrt(1/42) *
-                    (-4 + 105*eval_point - 300*(eval_point**2)
-                     + 210*(eval_point**3)),
-                    1/2 * np.sqrt(35/34) *
-                    (-5 + 48*eval_point - 105*(eval_point**2)
-                     + 64*(eval_point**3)),
-                    1/2 * np.sqrt(5/42) *
-                    (-16 + 105*eval_point - 192*(eval_point**2)
-                     + 105*(eval_point**3))]
+            return np.array([np.sqrt(15/34) *
+                             (1 + 4*eval_point - 30*(eval_point**2)
+                              + 28*(eval_point**3)),
+                             np.sqrt(1/42) *
+                             (-4 + 105*eval_point - 300*(eval_point**2)
+                              + 210*(eval_point**3)),
+                             1/2 * np.sqrt(35/34) *
+                             (-5 + 48*eval_point - 105*(eval_point**2)
+                             + 64*(eval_point**3)),
+                             1/2 * np.sqrt(5/42) *
+                             (-16 + 105*eval_point - 192*(eval_point**2)
+                             + 105*(eval_point**3))])
         if degree == 4:
-            return [np.sqrt(1/186) *
-                    (1 + 30*eval_point + 210*(eval_point**2)
-                     - 840*(eval_point**3) + 630*(eval_point**4)),
-                    0.5 * np.sqrt(1/38) *
-                    (-5 - 144*eval_point + 1155*(eval_point**2)
-                     - 2240*(eval_point**3) + 1260*(eval_point**4)),
-                    np.sqrt(35/14694) *
-                    (22 - 735*eval_point + 3504*(eval_point**2)
-                     - 5460*(eval_point**3) + 2700*(eval_point**4)),
-                    1/8 * np.sqrt(21/38) *
-                    (35 - 512*eval_point + 1890*(eval_point**2)
-                     - 2560*(eval_point**3) + 1155*(eval_point**4)),
-                    0.5 * np.sqrt(7/158) *
-                    (32 - 315*eval_point + 960*(eval_point**2)
-                     - 1155*(eval_point**3) + 480*(eval_point**4))]
+            return np.array([np.sqrt(1/186) *
+                             (1 + 30*eval_point + 210*(eval_point**2)
+                              - 840*(eval_point**3) + 630*(eval_point**4)),
+                             0.5 * np.sqrt(1/38) *
+                             (-5 - 144*eval_point + 1155*(eval_point**2)
+                             - 2240*(eval_point**3) + 1260*(eval_point**4)),
+                             np.sqrt(35/14694) *
+                             (22 - 735*eval_point + 3504*(eval_point**2)
+                             - 5460*(eval_point**3) + 2700*(eval_point**4)),
+                             1/8 * np.sqrt(21/38) *
+                             (35 - 512*eval_point + 1890*(eval_point**2)
+                             - 2560*(eval_point**3) + 1155*(eval_point**4)),
+                             0.5 * np.sqrt(7/158) *
+                             (32 - 315*eval_point + 960*(eval_point**2)
+                             - 1155*(eval_point**3) + 480*(eval_point**4))])
 
         raise ValueError('Invalid value: Alpert\'s wavelet is only available \
                          up to degree 4 for this application')
 
-    def _build_inverse_mass_matrix(self):
+    def _build_inverse_mass_matrix(self) -> ndarray:
         """Construct inverse mass matrix.
 
         Returns
@@ -355,7 +358,7 @@ class OrthonormalLegendre(Legendre):
 
     @property
     @cache
-    def basis_projections(self):
+    def basis_projections(self) -> Tuple[ndarray, ndarray]:
         """Return basis projection.
 
         Construct matrices containing the integrals of the
@@ -374,7 +377,8 @@ class OrthonormalLegendre(Legendre):
         right_basis_projection = self._build_basis_matrix(z, 0.5 * (z + 1))
         return left_basis_projection, right_basis_projection
 
-    def _build_basis_matrix(self, first_param, second_param):
+    def _build_basis_matrix(self, first_param: float, second_param: float) \
+            -> ndarray:
         """Construct a basis matrix.
 
         Parameters
@@ -399,11 +403,11 @@ class OrthonormalLegendre(Legendre):
                                   (z, -1, 1))
                 row.append(np.float64(entry))
             matrix.append(row)
-        return matrix
+        return np.array(matrix)
 
     @property
     @cache
-    def multiwavelet_projections(self):
+    def multiwavelet_projections(self) -> Tuple[ndarray, ndarray]:
         """Return wavelet projection.
 
         Construct matrices containing the integrals of the
@@ -424,8 +428,9 @@ class OrthonormalLegendre(Legendre):
             z, 0.5*(z+1), False)
         return left_wavelet_projection, right_wavelet_projection
 
-    def _build_multiwavelet_matrix(self, first_param, second_param,
-                                   is_left_matrix):
+    def _build_multiwavelet_matrix(self, first_param: float,
+                                   second_param: float, is_left_matrix: bool) \
+            -> ndarray:
         """Construct a multiwavelet matrix.
 
         Parameters
@@ -455,10 +460,10 @@ class OrthonormalLegendre(Legendre):
                     entry = entry * (-1)**(j + self._polynomial_degree + 1)
                 row.append(np.float64(entry))
             matrix.append(row)
-        return matrix
+        return np.array(matrix)
 
-    def calculate_cell_average(self, projection, stencil_length,
-                               add_reconstructions=True):
+    def calculate_cell_average(self, projection: ndarray, stencil_length: int,
+                               add_reconstructions: bool = True) -> ndarray:
         """Calculate cell averages for a given projection.
 
         Calculate the cell averages of all cells in a projection.
@@ -502,7 +507,8 @@ class OrthonormalLegendre(Legendre):
                                 cell_averages[:, middle_idx+1:]))))
         return np.array(list(map(np.float64, cell_averages)))
 
-    def _calculate_reconstructions(self, projection):
+    def _calculate_reconstructions(self, projection: ndarray) \
+            -> Tuple[list, list]:
         """Calculate left and right reconstructions for a given projection.
 
         Parameters