diff --git a/coverage-report.pdf b/coverage-report.pdf
new file mode 100644
index 0000000..20adb06
Binary files /dev/null and b/coverage-report.pdf differ
diff --git a/diffusion2d.py b/diffusion2d.py
index 51a07f2..c17221c 100644
--- a/diffusion2d.py
+++ b/diffusion2d.py
@@ -44,8 +44,12 @@ def initialize_domain(self, w=10., h=10., dx=0.1, dy=0.1):
         self.dy = dy
         self.nx = int(w / dx)
         self.ny = int(h / dy)
+        assert isinstance(self.w, float), "type(w) must be float"
+        assert isinstance(self.h, float), "type(h) must be float"
+        assert isinstance(self.dx, float), "type(dx) must be float"
+        assert isinstance(self.dy, float), "type(dy) must be float"
 
-    def initialize_physical_parameters(self, d=4., T_cold=300, T_hot=700):
+    def initialize_physical_parameters(self, d=4., T_cold=300., T_hot=700.):
         self.D = d
         self.T_cold = T_cold
         self.T_hot = T_hot
@@ -56,6 +60,10 @@ def initialize_physical_parameters(self, d=4., T_cold=300, T_hot=700):
 
         print("dt = {}".format(self.dt))
 
+        assert isinstance(self.D, float), "type(D) must be float"
+        assert isinstance(self.T_cold, float), "type(T_cold) must be float"
+        assert isinstance(self.T_hot, float), "type(T_hot) must be float"
+
     def set_initial_condition(self):
         u = self.T_cold * np.ones((self.nx, self.ny))
 
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..a6c4ede
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,3 @@
+pytest
+numpy
+matplotlib
diff --git a/tests/integration/test_diffusion2d.py b/tests/integration/test_diffusion2d.py
index fd026b4..a3ebcfc 100644
--- a/tests/integration/test_diffusion2d.py
+++ b/tests/integration/test_diffusion2d.py
@@ -3,17 +3,40 @@
 """
 
 from diffusion2d import SolveDiffusion2D
-
+import numpy as np
 
 def test_initialize_physical_parameters():
     """
     Checks function SolveDiffusion2D.initialize_domain
     """
     solver = SolveDiffusion2D()
+    solver.initialize_domain(dx=0.2, dy=0.6)
 
+    d = 2.
+    T_cold = 150.
+    T_hot = 350.
+    dx2, dy2 = solver.dx * solver.dx, solver.dy * solver.dy
+    dt = dx2 * dy2 / (2 * d * (dx2 + dy2))
+    solver.initialize_physical_parameters(d=d, T_cold=T_cold, T_hot=T_hot)
+    assert dt == solver.dt, "wrong value of dt"
 
 def test_set_initial_condition():
     """
     Checks function SolveDiffusion2D.get_initial_function
     """
     solver = SolveDiffusion2D()
+    solver.initialize_domain(dx=0.2, dy=0.6)
+    solver.initialize_physical_parameters(d=2., T_cold=150., T_hot=350.)
+    u = 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[i, j] = solver.T_hot
+
+    ground_truth = u.copy()
+
+    assert np.array_equal(ground_truth, solver.set_initial_condition()), "initial condition not set correctly"
\ No newline at end of file
diff --git a/tests/unit/test_diffusion2d_functions.py b/tests/unit/test_diffusion2d_functions.py
index c4277ff..fa4f4ec 100644
--- a/tests/unit/test_diffusion2d_functions.py
+++ b/tests/unit/test_diffusion2d_functions.py
@@ -3,13 +3,23 @@
 """
 
 from diffusion2d import SolveDiffusion2D
-
+import numpy as np
 
 def test_initialize_domain():
     """
     Check function SolveDiffusion2D.initialize_domain
     """
     solver = SolveDiffusion2D()
+    w = 5.
+    h = 20.
+    dx = 0.2
+    dy = 0.6
+    nx = int(w / dx)
+    ny = int(h / dy)
+
+    solver.initialize_domain(w=w, h=h, dx=dx, dy=dy)
+    assert nx == solver.nx, "wrong value of nx"
+    assert ny == solver.ny, "wrong value of ny"
 
 
 def test_initialize_physical_parameters():
@@ -17,6 +27,17 @@ def test_initialize_physical_parameters():
     Checks function SolveDiffusion2D.initialize_domain
     """
     solver = SolveDiffusion2D()
+    
+    solver.dx = 0.2
+    solver.dy = 0.6
+
+    d = 2.
+    T_cold = 150.
+    T_hot = 350.
+    dx2, dy2 = solver.dx * solver.dx, solver.dy * solver.dy
+    dt = dx2 * dy2 / (2 * d * (dx2 + dy2))
+    solver.initialize_physical_parameters(d=d, T_cold=T_cold, T_hot=T_hot)
+    assert dt == solver.dt, "wrong value of dt"
 
 
 def test_set_initial_condition():
@@ -24,3 +45,30 @@ def test_set_initial_condition():
     Checks function SolveDiffusion2D.get_initial_function
     """
     solver = SolveDiffusion2D()
+    solver.w = 5.
+    solver.h = 20.
+    solver.dx = 0.2
+    solver.dy = 0.6
+    solver.nx = int(solver.w / solver.dx)
+    solver.ny = int(solver.h / solver.dy)
+    solver.dx = 0.2
+    solver.dy = 0.6
+
+    solver.d = 2.
+    solver.T_cold = 150.
+    solver.T_hot = 350.
+
+    u = 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[i, j] = solver.T_hot
+
+    ground_truth = u.copy()
+
+    assert np.array_equal(ground_truth, solver.set_initial_condition()), "initial condition not set correctly"
+
diff --git a/tox.toml b/tox.toml
new file mode 100644
index 0000000..f7bf328
--- /dev/null
+++ b/tox.toml
@@ -0,0 +1,13 @@
+requires = ["tox>=4"]
+env_list = ["pytest_testing"]
+
+[env.pytest_testing]
+description = "Run pytest"
+setenv = { PYTHONPATH = ".", MPLBACKEND = "Agg" } # had to add this due to my python configuration on mac
+deps = [
+    "pytest>=8",
+    "-rrequirements.txt"
+]
+commands = [
+    ["pytest", "-q"]
+]
\ No newline at end of file