diff --git a/Coverage report.pdf b/Coverage report.pdf
new file mode 100644
index 0000000..9297440
Binary files /dev/null and b/Coverage report.pdf differ
diff --git a/diffusion2d.py b/diffusion2d.py
index 51a07f2..cc2b3d2 100644
--- a/diffusion2d.py
+++ b/diffusion2d.py
@@ -38,6 +38,10 @@ def __init__(self):
         self.dt = None
 
     def initialize_domain(self, w=10., h=10., dx=0.1, dy=0.1):
+        assert isinstance(w, float)
+        assert isinstance(h, float)
+        assert isinstance(dx, float)
+        assert isinstance(dy, float)
         self.w = w
         self.h = h
         self.dx = dx
@@ -45,7 +49,10 @@ def initialize_domain(self, w=10., h=10., dx=0.1, dy=0.1):
         self.nx = int(w / dx)
         self.ny = int(h / dy)
 
-    def initialize_physical_parameters(self, d=4., T_cold=300, T_hot=700):
+    def initialize_physical_parameters(self, d=4., T_cold=300.0, T_hot=700.0):
+        assert isinstance(d, float)
+        assert isinstance(T_cold, float)
+        assert isinstance(T_hot, float)
         self.D = d
         self.T_cold = T_cold
         self.T_hot = T_hot
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..c8368a0
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,9 @@
+# Testing & Tools
+pytest>=7.0.0
+pytest-cov>=4.0.0
+tox>=4.0.0
+coverage>=7.0.0
+
+# Data Science / Plotting
+numpy
+matplotlib
\ No newline at end of file
diff --git a/tests/integration/test_diffusion2d.py b/tests/integration/test_diffusion2d.py
index fd026b4..9f144b3 100644
--- a/tests/integration/test_diffusion2d.py
+++ b/tests/integration/test_diffusion2d.py
@@ -1,19 +1,48 @@
 """
 Tests for functionality checks in class SolveDiffusion2D
 """
-
+import numpy as np
+import pytest
 from diffusion2d import SolveDiffusion2D
 
+@pytest.fixture
+def initial_values():
+    return {
+        "w": 5.0, "h": 5.0, "dx": 0.2, "dy": 0.2, "d": 3.5, "T_cold": 290.0, "T_hot": 750.0}
 
-def test_initialize_physical_parameters():
+def test_initialize_physical_parameters(initial_values):
     """
     Checks function SolveDiffusion2D.initialize_domain
     """
     solver = SolveDiffusion2D()
+    solver.initialize_domain(w=initial_values["w"], h=initial_values["h"], dx=initial_values["dx"], dy=initial_values["dy"])
+
+    d = initial_values["d"]
+
+    dx2, dy2 = solver.dx * solver.dx, solver.dy * solver.dy
+    dt = dx2 * dy2 / (2 * d * (dx2 + dy2))
 
+    solver.initialize_physical_parameters(d=initial_values["d"], T_cold=initial_values["T_cold"], T_hot=initial_values["T_hot"])
+    assert solver.dt == dt
 
-def test_set_initial_condition():
+def test_set_initial_condition(initial_values):
     """
     Checks function SolveDiffusion2D.get_initial_function
     """
     solver = SolveDiffusion2D()
+    solver.initialize_domain(w=initial_values["w"], h=initial_values["h"], dx=initial_values["dx"], dy=initial_values["dy"])
+    solver.initialize_physical_parameters(d=initial_values["d"], T_cold=initial_values["T_cold"], T_hot=initial_values["T_hot"])    
+    u_test = solver.T_cold * np.ones((solver.nx, solver.ny))
+    r, cx, cy = 2, 5, 5
+    r2 = r ** 2
+    for i in range(solver.nx):
+        for j in range(solver.ny):
+            p2 = (i * solver.dx - cx) ** 2 + (j * solver.dy - cy) ** 2
+            if p2 < r2:
+                u_test[i, j] = solver.T_hot
+
+    ground_truth = u_test.copy()
+
+    test_value = solver.set_initial_condition()
+
+    assert np.array_equal(ground_truth, test_value)
\ No newline at end of file
diff --git a/tests/unit/test_diffusion2d_functions.py b/tests/unit/test_diffusion2d_functions.py
index c4277ff..a3741bd 100644
--- a/tests/unit/test_diffusion2d_functions.py
+++ b/tests/unit/test_diffusion2d_functions.py
@@ -2,6 +2,7 @@
 Tests for functions in class SolveDiffusion2D
 """
 
+import pytest, numpy as np
 from diffusion2d import SolveDiffusion2D
 
 
@@ -10,6 +11,9 @@ def test_initialize_domain():
     Check function SolveDiffusion2D.initialize_domain
     """
     solver = SolveDiffusion2D()
+    solver.initialize_domain(w=5.0, h=5.0, dx=0.2, dy=0.2)
+    assert solver.nx == 25
+    assert solver.ny == 25
 
 
 def test_initialize_physical_parameters():
@@ -17,6 +21,11 @@ def test_initialize_physical_parameters():
     Checks function SolveDiffusion2D.initialize_domain
     """
     solver = SolveDiffusion2D()
+    solver.dx = 0.2
+    solver.dy = 0.2
+    solver.D = 3.5
+    solver.initialize_physical_parameters(d=3.5, T_cold=290.0, T_hot=750.0)
+    assert solver.dt == pytest.approx(0.002857142857142857, rel=1e-6)
 
 
 def test_set_initial_condition():
@@ -24,3 +33,24 @@ def test_set_initial_condition():
     Checks function SolveDiffusion2D.get_initial_function
     """
     solver = SolveDiffusion2D()
+    solver.nx = 100
+    solver.ny = 100
+    solver.dx = 0.1
+    solver.dy = 0.1
+    solver.T_cold = 300.0
+    solver.T_hot = 750.0
+
+    u = solver.set_initial_condition()
+
+
+    u_test = solver.T_cold * np.ones((solver.nx, solver.ny))
+    r, cx, cy = 2, 5, 5
+    r2 = r ** 2
+    for i in range(solver.nx):
+        for j in range(solver.ny):
+            p2 = (i * solver.dx - cx) ** 2 + (j * solver.dy - cy) ** 2
+            if p2 < r2:
+                u_test[i, j] = solver.T_hot
+
+    ground_truth = u_test.copy()
+    assert np.array_equal(ground_truth, u)
\ No newline at end of file
diff --git a/tox.toml b/tox.toml
new file mode 100644
index 0000000..050ccf9
--- /dev/null
+++ b/tox.toml
@@ -0,0 +1,12 @@
+[tox]
+env_list = ["pytest"]
+
+[env_common]
+# Common settings to avoid repetition
+base_python = ["python3.9"]
+deps = ["-rrequirements.txt"]
+
+[env.pytest_env]
+description = "Run tests using pytest"
+deps = ["-rrequirements.txt"]
+commands = [["pytest"]]
\ No newline at end of file