diff --git a/coverage-report.pdf b/coverage-report.pdf
new file mode 100644
index 0000000..bab5ace
Binary files /dev/null and b/coverage-report.pdf differ
diff --git a/diffusion2d.py b/diffusion2d.py
index 51a07f2..7d6ff53 100644
--- a/diffusion2d.py
+++ b/diffusion2d.py
@@ -38,6 +38,7 @@ def __init__(self):
         self.dt = None
 
     def initialize_domain(self, w=10., h=10., dx=0.1, dy=0.1):
+        assert type(w) is float and type(h) is float and type(dx) is float and type(dy) is float
         self.w = w
         self.h = h
         self.dx = dx
@@ -45,7 +46,8 @@ 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., T_hot=700.):
+        assert type(d) is float and type(T_cold) is float and type(T_hot) is 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..8f87fb4
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,4 @@
+numpy
+matplotlib
+pytest
+coverage
\ No newline at end of file
diff --git a/tests/integration/test_diffusion2d.py b/tests/integration/test_diffusion2d.py
index fd026b4..d2d875c 100644
--- a/tests/integration/test_diffusion2d.py
+++ b/tests/integration/test_diffusion2d.py
@@ -3,17 +3,66 @@
 """
 
 from diffusion2d import SolveDiffusion2D
+import pytest
+import numpy as np
 
 
-def test_initialize_physical_parameters():
+@pytest.fixture
+def solver():
+    solver = SolveDiffusion2D()
+    return solver
+
+
+def test_initialize_physical_parameters(solver):
     """
     Checks function SolveDiffusion2D.initialize_domain
     """
-    solver = SolveDiffusion2D()
+    input_w = 100.
+    input_h = 200.
+    input_dx = 0.1
+    input_dy = 0.1
+
+    expected_T_cold = 100.
+    expected_T_hot = 1000.
+    expected_D = 5.
+
+    expected_dt = 5 * 10e-5
+
+    solver.initialize_domain(input_w, input_h, input_dx, input_dy)
+    solver.initialize_physical_parameters(expected_D, expected_T_cold, expected_T_hot)
 
+    assert solver.dt == pytest.approx(expected_dt, rel=1e-12, abs=0.0)
 
-def test_set_initial_condition():
+
+def test_set_initial_condition(solver):
     """
     Checks function SolveDiffusion2D.get_initial_function
     """
-    solver = SolveDiffusion2D()
+    expected_nx = 1000
+    expected_ny = 2000
+
+    input_w = 100.
+    input_h = 200.
+    input_dx = 0.1
+    input_dy = 0.1
+
+    expected_T_cold = 100.
+    expected_T_hot = 1000.
+    expected_D = 5.
+    
+    solver.initialize_domain(input_w, input_h, input_dx, input_dy)
+    solver.initialize_physical_parameters(expected_D, expected_T_cold, expected_T_hot)
+    calculated_u = solver.set_initial_condition()
+
+    expected_u = solver.T_cold * np.ones((expected_nx, expected_ny))
+
+    r, cx, cy = 2, 5, 5
+    r2 = r ** 2
+    for i in range(expected_nx):
+        for j in range(expected_ny):
+            p2 = (i * input_dx - cx) ** 2 + (j * input_dy - cy) ** 2
+            if p2 < r2:
+                expected_u[i, j] = solver.T_hot
+
+    assert np.array_equal(expected_u, calculated_u)
+    
\ No newline at end of file
diff --git a/tests/unit/test_diffusion2d_functions.py b/tests/unit/test_diffusion2d_functions.py
index c4277ff..90e5251 100644
--- a/tests/unit/test_diffusion2d_functions.py
+++ b/tests/unit/test_diffusion2d_functions.py
@@ -3,24 +3,73 @@
 """
 
 from diffusion2d import SolveDiffusion2D
+import pytest
+import numpy as np
 
 
-def test_initialize_domain():
+@pytest.fixture
+def solver():
+    solver = SolveDiffusion2D()
+    return solver
+
+def test_initialize_domain(solver: SolveDiffusion2D):
     """
     Check function SolveDiffusion2D.initialize_domain
     """
-    solver = SolveDiffusion2D()
+    expected_nx = 50
+    expected_ny = 100
 
+    input_w = 100.
+    input_h = 200.
+    input_dx = 2.
+    input_dy = 2.
+    
+    solver.initialize_domain(input_w, input_h, input_dx, input_dy)
+    assert solver.nx == expected_nx 
+    assert solver.ny == expected_ny
+    
 
-def test_initialize_physical_parameters():
+def test_initialize_physical_parameters(solver: SolveDiffusion2D):
     """
     Checks function SolveDiffusion2D.initialize_domain
     """
-    solver = SolveDiffusion2D()
+    expected_T_cold = 100.
+    expected_T_hot = 1000.
+    expected_D = 5.
+
+    expected_dt = 5 * 10e-5
+
+    solver.dx = 0.1
+    solver.dy = 0.1
 
+    solver.initialize_physical_parameters(expected_D, expected_T_cold, expected_T_hot)
 
-def test_set_initial_condition():
+    assert solver.T_cold == expected_T_cold 
+    assert solver.T_hot == expected_T_hot
+    assert solver.D == expected_D
+    assert solver.dt == pytest.approx(expected_dt, rel=1e-12, abs=0.0)
+
+
+def test_set_initial_condition(solver: SolveDiffusion2D):
     """
     Checks function SolveDiffusion2D.get_initial_function
     """
-    solver = SolveDiffusion2D()
+    solver.nx = 200
+    solver.ny = 200
+    solver.dx = 0.1
+    solver.dy = 0.1
+    solver.T_cold = 50
+    solver.T_hot = 200
+
+    expected_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:
+                expected_u[i, j] = solver.T_hot
+
+    calculated_u = solver.set_initial_condition()
+    assert np.array_equal(expected_u, calculated_u)
diff --git a/tox.toml b/tox.toml
new file mode 100644
index 0000000..493bd67
--- /dev/null
+++ b/tox.toml
@@ -0,0 +1,7 @@
+requires = ["tox>=4"]
+env_list = ["pytest_testing"]
+
+[env.pytest_testing]
+description = "Run pytest"
+deps = ["-rrequirements.txt"]
+commands = [["python", "-m", "pytest"]]