diff --git a/coverage-report.pdf b/coverage-report.pdf
new file mode 100644
index 0000000..3c96873
Binary files /dev/null and b/coverage-report.pdf differ
diff --git a/diffusion2d.py b/diffusion2d.py
index 51a07f2..ef46d10 100644
--- a/diffusion2d.py
+++ b/diffusion2d.py
@@ -38,14 +38,17 @@ def __init__(self):
         self.dt = None
 
     def initialize_domain(self, w=10., h=10., dx=0.1, dy=0.1):
+
         self.w = w
         self.h = h
         self.dx = dx
         self.dy = dy
         self.nx = int(w / dx)
         self.ny = int(h / dy)
+        assert isinstance(self.w, float) and isinstance(self.h, float) and isinstance(self.dx, float) and isinstance(self.dy, float)
+        assert isinstance(self.nx, int) and isinstance(self.ny, int)
 
-    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 +59,8 @@ def initialize_physical_parameters(self, d=4., T_cold=300, T_hot=700):
 
         print("dt = {}".format(self.dt))
 
+        assert isinstance(self.D, float) and isinstance(self.T_cold, float) and isinstance(self.T_hot, float) and isinstance(self.dt, 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..7552c4e
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,3 @@
+numpy
+pytest
+matplotlib
\ No newline at end of file
diff --git a/tests/integration/test_diffusion2d.py b/tests/integration/test_diffusion2d.py
index fd026b4..ae6a5be 100644
--- a/tests/integration/test_diffusion2d.py
+++ b/tests/integration/test_diffusion2d.py
@@ -3,6 +3,7 @@
 """
 
 from diffusion2d import SolveDiffusion2D
+import numpy as np
 
 
 def test_initialize_physical_parameters():
@@ -11,9 +12,53 @@ def test_initialize_physical_parameters():
     """
     solver = SolveDiffusion2D()
 
+    w = 12.0
+    h = 8.0
+    dx = 0.2
+    dy = 0.2
+    d = 5.
+    T_cold = 290.0
+    T_hot = 750.0
+
+    dx2, dy2 = dx * dx, dy * dy
+    dt = dx2 * dy2 / (2 * d * (dx2 + dy2))
+
+    solver.initialize_domain(dx=dx, dy=dy)
+    solver.initialize_physical_parameters(d=d, T_cold=T_cold, T_hot=T_hot)
+
+    assert solver.D == d
+    assert solver.T_cold == T_cold
+    assert solver.T_hot == T_hot
+    assert solver.dt == dt
+
 
 def test_set_initial_condition():
     """
     Checks function SolveDiffusion2D.get_initial_function
     """
     solver = SolveDiffusion2D()
+
+    w = 12.0
+    h = 8.0
+    d = 5.
+    dx = 0.2
+    dy = 0.2
+    T_cold = 290.0
+    T_hot = 750.0
+    nx = int(w / dx)
+    ny = int(h / dy)
+
+    u = T_cold * np.ones((nx, ny))
+    r, cx, cy = 2, 5, 5
+    r2 = r ** 2
+    for i in range(nx):
+        for j in range(ny):
+            p2 = (i * dx - cx) ** 2 + (j * dy - cy) ** 2
+            if p2 < r2:
+                u[i, j] = T_hot
+
+    solver.initialize_domain(w=w, h=h, dx=dx, dy=dy)
+    solver.initialize_physical_parameters(d=d, T_cold=T_cold, T_hot=T_hot)
+    u_return = solver.set_initial_condition()
+
+    assert np.array_equal(u, u_return)
diff --git a/tests/unit/test_diffusion2d_functions.py b/tests/unit/test_diffusion2d_functions.py
index c4277ff..ca369d9 100644
--- a/tests/unit/test_diffusion2d_functions.py
+++ b/tests/unit/test_diffusion2d_functions.py
@@ -3,7 +3,7 @@
 """
 
 from diffusion2d import SolveDiffusion2D
-
+import numpy as np
 
 def test_initialize_domain():
     """
@@ -11,6 +11,23 @@ def test_initialize_domain():
     """
     solver = SolveDiffusion2D()
 
+    w = 12.0
+    h = 8.0
+    dx = 0.2
+    dy = 0.2
+
+    nx = int(w / dx)
+    ny = int(h / dy)
+
+    solver.initialize_domain(w=w, h=h, dx=dx, dy=dy)
+
+    assert solver.w == w
+    assert solver.h == h
+    assert solver.dx == dx
+    assert solver.dy == dy
+    assert solver.nx == nx
+    assert solver.ny == ny
+
 
 def test_initialize_physical_parameters():
     """
@@ -18,9 +35,52 @@ def test_initialize_physical_parameters():
     """
     solver = SolveDiffusion2D()
 
+    dx = 0.2
+    dy = 0.2
+    d = 5.
+    T_cold = 290.0
+    T_hot = 750.0
+
+    dx2, dy2 = dx * dx, dy * dy
+    dt = dx2 * dy2 / (2 * d * (dx2 + dy2))
+
+    solver.dx = dx
+    solver.dy = dy
+    solver.initialize_physical_parameters(d=d, T_cold=T_cold, T_hot=T_hot)
+
+    assert solver.D == d
+    assert solver.T_cold == T_cold
+    assert solver.T_hot == T_hot
+    assert solver.dt == dt
 
 def test_set_initial_condition():
     """
     Checks function SolveDiffusion2D.get_initial_function
     """
     solver = SolveDiffusion2D()
+
+    dx = 0.2
+    dy = 0.2
+    T_cold = 290.0
+    T_hot = 750.0
+    nx = int(12.0 / dx)
+    ny = int(8.0 / dy)
+
+    u = T_cold * np.ones((nx, ny))
+    r, cx, cy = 2, 5, 5
+    r2 = r ** 2
+    for i in range(nx):
+        for j in range(ny):
+            p2 = (i * dx - cx) ** 2 + (j * dy - cy) ** 2
+            if p2 < r2:
+                u[i, j] = T_hot
+
+    solver.nx = nx
+    solver.ny = ny
+    solver.dx = dx
+    solver.dy = dy
+    solver.T_cold = T_cold
+    solver.T_hot = T_hot
+    u_return = solver.set_initial_condition()
+
+    assert np.array_equal(u, u_return)
\ No newline at end of file
diff --git a/tox.toml b/tox.toml
new file mode 100644
index 0000000..340f584
--- /dev/null
+++ b/tox.toml
@@ -0,0 +1,12 @@
+requires = ["tox>=4"]
+env_list = ["py314"]
+
+[env_run_base]
+description = "run tests"
+deps = ["-r requirements.txt"]
+
+set_env = {PYTHONPATH = "{tox_root}" }
+
+commands = [
+    ["pytest"]
+]