diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..4f80bf4
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+.venv
+__pycache__
+htmlcov
+.pytest_cache
+.coverage
\ No newline at end of file
diff --git a/coverage-report.pdf b/coverage-report.pdf
new file mode 100644
index 0000000..1e553e9
Binary files /dev/null and b/coverage-report.pdf differ
diff --git a/diffusion2d.py b/diffusion2d.py
index 51a07f2..8919154 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(isinstance(w, float) and isinstance(h, float) and isinstance(dx, float) and isinstance(dy, float))
         self.w = w
         self.h = h
         self.dx = dx
@@ -45,7 +46,9 @@ 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) and isinstance(T_cold, float) and 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..e072a28
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,3 @@
+matplotlib
+numpy
+pytest
\ No newline at end of file
diff --git a/tests/integration/test_diffusion2d.py b/tests/integration/test_diffusion2d.py
index fd026b4..a0647f0 100644
--- a/tests/integration/test_diffusion2d.py
+++ b/tests/integration/test_diffusion2d.py
@@ -4,6 +4,8 @@
 
 from diffusion2d import SolveDiffusion2D
 
+import pytest 
+import numpy as np
 
 def test_initialize_physical_parameters():
     """
@@ -11,9 +13,41 @@ def test_initialize_physical_parameters():
     """
     solver = SolveDiffusion2D()
 
+    solver.initialize_domain()
+    solver.initialize_physical_parameters(5.0, 300.0, 700.0)
+
+    solver.initialize_physical_parameters(5.0, 300.0, 700.0)
+    assert solver.dt == pytest.approx(0.0005, rel=1e-4)
+
+    solver.initialize_physical_parameters(5.0, 500.0, 800.0)
+    assert solver.dt == pytest.approx(0.0005, rel=1e-4)
+
+    with pytest.raises(AssertionError):
+        solver.initialize_physical_parameters(5, 500, 800)
+
 
 def test_set_initial_condition():
     """
     Checks function SolveDiffusion2D.get_initial_function
     """
     solver = SolveDiffusion2D()
+
+    solver.initialize_domain(8.0, 8.0, 1.0, 1.0)
+    solver.initialize_physical_parameters(4.0, 0.0, 1.0)
+
+    # 8x8 domain with circle at (5,5)
+    expected_field = np.array(
+        [[0., 0., 0., 0., 0., 0., 0., 0.],
+         [0., 0., 0., 0., 0., 0., 0., 0.],
+         [0., 0., 0., 0., 0., 0., 0., 0.],
+         [0., 0., 0., 0., 0., 0., 0., 0.],
+         [0., 0., 0., 0., 1., 1., 1., 0.],
+         [0., 0., 0., 0., 1., 1., 1., 0.],
+         [0., 0., 0., 0., 1., 1., 1., 0.],
+         [0., 0., 0., 0., 0., 0., 0., 0.]]
+    )
+
+    actual_field = solver.set_initial_condition()
+
+    error_norm = np.linalg.norm(expected_field - actual_field)
+    assert 0 == pytest.approx(error_norm, abs=1e-14)
diff --git a/tests/unit/test_diffusion2d_functions.py b/tests/unit/test_diffusion2d_functions.py
index c4277ff..24217e1 100644
--- a/tests/unit/test_diffusion2d_functions.py
+++ b/tests/unit/test_diffusion2d_functions.py
@@ -4,6 +4,8 @@
 
 from diffusion2d import SolveDiffusion2D
 
+import pytest 
+import numpy as np
 
 def test_initialize_domain():
     """
@@ -11,6 +13,14 @@ def test_initialize_domain():
     """
     solver = SolveDiffusion2D()
 
+    solver.initialize_domain(1.5, 2.0, 0.2, 0.2)
+
+    assert solver.nx == 7
+    assert solver.ny == 10
+
+    with pytest.raises(AssertionError):
+        solver.initialize_domain(1, 2, 0.2, 0.2)
+
 
 def test_initialize_physical_parameters():
     """
@@ -18,9 +28,47 @@ def test_initialize_physical_parameters():
     """
     solver = SolveDiffusion2D()
 
+    solver.dx = 0.1
+    solver.dy = 0.1
+    solver.nx = 100
+    solver.ny = 100
+
+    solver.initialize_physical_parameters(5.0, 300.0, 700.0)
+    assert solver.dt == pytest.approx(0.0005, rel=1e-4)
+
+    solver.initialize_physical_parameters(5.0, 500.0, 800.0)
+    assert solver.dt == pytest.approx(0.0005, rel=1e-4)
+
+    with pytest.raises(AssertionError):
+        solver.initialize_physical_parameters(5, 500, 800)
+
 
 def test_set_initial_condition():
     """
     Checks function SolveDiffusion2D.get_initial_function
     """
     solver = SolveDiffusion2D()
+
+    solver.T_cold = 0.0
+    solver.T_hot = 1.0
+    solver.dx = 1.0
+    solver.dy = 1.0
+    solver.nx = 8
+    solver.ny = 8
+
+    # 8x8 domain with circle at (5,5)
+    expected_field = np.array(
+        [[0., 0., 0., 0., 0., 0., 0., 0.],
+         [0., 0., 0., 0., 0., 0., 0., 0.],
+         [0., 0., 0., 0., 0., 0., 0., 0.],
+         [0., 0., 0., 0., 0., 0., 0., 0.],
+         [0., 0., 0., 0., 1., 1., 1., 0.],
+         [0., 0., 0., 0., 1., 1., 1., 0.],
+         [0., 0., 0., 0., 1., 1., 1., 0.],
+         [0., 0., 0., 0., 0., 0., 0., 0.]]
+    )
+
+    actual_field = solver.set_initial_condition()
+
+    error_norm = np.linalg.norm(expected_field - actual_field)
+    assert 0 == pytest.approx(error_norm, abs=1e-14)
diff --git a/tox.toml b/tox.toml
new file mode 100644
index 0000000..b1d0b4f
--- /dev/null
+++ b/tox.toml
@@ -0,0 +1,8 @@
+[tox]
+envlist = ["test"]
+requires = ["tox>=4"]
+
+[env.test]
+description = "Run test suite with pytest."
+deps = "-r requirements.txt"
+commands = [["python","-m","pytest"]]