diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..297b916
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,11 @@
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.pyc
+
+# OS files
+.DS_Store
+
+# Testing artifacts
+.coverage
+htmlcov/
+.tox/
\ No newline at end of file
diff --git a/coverage-report.pdf b/coverage-report.pdf
new file mode 100644
index 0000000..7408f11
Binary files /dev/null and b/coverage-report.pdf differ
diff --git a/diffusion2d.py b/diffusion2d.py
index 51a07f2..7db10be 100644
--- a/diffusion2d.py
+++ b/diffusion2d.py
@@ -38,6 +38,11 @@ def __init__(self):
         self.dt = None
 
     def initialize_domain(self, w=10., h=10., dx=0.1, dy=0.1):
+        assert isinstance(w, float), "w must be a float"
+        assert isinstance(h, float), "h must be a float"
+        assert isinstance(dx, float), "dx must be a float"
+        assert isinstance(dy, float), "dy must be a float"
+
         self.w = w
         self.h = h
         self.dx = dx
@@ -45,7 +50,11 @@ 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 isinstance(d, float), "d must be a float"
+        assert isinstance(T_cold, float), "T_cold must be a float"
+        assert isinstance(T_hot, float), "T_hot must be a 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..2eae7b2 100644
--- a/tests/integration/test_diffusion2d.py
+++ b/tests/integration/test_diffusion2d.py
@@ -2,18 +2,49 @@
 Tests for functionality checks in class SolveDiffusion2D
 """
 
+"""
+Integration tests for SolveDiffusion2D
+"""
+
+import pytest
+import numpy as np
+import numpy.testing as npt
+
 from diffusion2d import SolveDiffusion2D
 
 
 def test_initialize_physical_parameters():
-    """
-    Checks function SolveDiffusion2D.initialize_domain
-    """
+    """Initialize domain and physical parameters, then verify dt."""
     solver = SolveDiffusion2D()
+    solver.initialize_domain(w=10.0, h=10.0, dx=0.5, dy=0.5)
+    solver.initialize_physical_parameters(d=4.0, T_cold=100.0, T_hot=500.0)
+
+    dx2 = solver.dx * solver.dx
+    dy2 = solver.dy * solver.dy
+    expected = (dx2 * dy2) / (2 * solver.D * (dx2 + dy2))
+
+    assert solver.dt == pytest.approx(expected)
 
 
-def test_set_initial_condition():
-    """
-    Checks function SolveDiffusion2D.get_initial_function
-    """
+def test_set_initial_conditions():
+    """Initialize domain and physical parameters, then verify initial condition array."""
     solver = SolveDiffusion2D()
+    solver.initialize_domain(w=10.0, h=10.0, dx=0.5, dy=0.5)
+    solver.initialize_physical_parameters(d=4.0, T_cold=100.0, T_hot=500.0)
+
+    # call method and store result
+    u = solver.set_initial_condition()
+
+    # build expected array the same way as in set_initial_condition
+    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
+
+    npt.assert_array_equal(u, expected_u)
+
+
diff --git a/tests/unit/test_diffusion2d_functions.py b/tests/unit/test_diffusion2d_functions.py
index c4277ff..8e8509a 100644
--- a/tests/unit/test_diffusion2d_functions.py
+++ b/tests/unit/test_diffusion2d_functions.py
@@ -2,6 +2,9 @@
 Tests for functions in class SolveDiffusion2D
 """
 
+import pytest
+import unittest
+
 from diffusion2d import SolveDiffusion2D
 
 
@@ -10,6 +13,9 @@ def test_initialize_domain():
     Check function SolveDiffusion2D.initialize_domain
     """
     solver = SolveDiffusion2D()
+    solver.initialize_domain(w=20.0, h=10.0, dx=0.5, dy=0.5)
+    assert solver.nx == 40
+    assert solver.ny == 20
 
 
 def test_initialize_physical_parameters():
@@ -17,6 +23,15 @@ def test_initialize_physical_parameters():
     Checks function SolveDiffusion2D.initialize_domain
     """
     solver = SolveDiffusion2D()
+    solver.dx = 0.5
+    solver.dy = 0.5
+    solver.initialize_physical_parameters(d=4.0, T_cold=100.0, T_hot=500.0)
+
+    expected = (solver.dx * solver.dx) * (solver.dy * solver.dy) / (
+        2 * 4.0 * (solver.dx * solver.dx + solver.dy * solver.dy)
+    )
+
+    assert solver.dt == pytest.approx(expected)
 
 
 def test_set_initial_condition():
@@ -24,3 +39,46 @@ def test_set_initial_condition():
     Checks function SolveDiffusion2D.get_initial_function
     """
     solver = SolveDiffusion2D()
+    solver.nx = 40
+    solver.ny = 20
+    solver.dx = 0.5
+    solver.dy = 0.5
+    solver.T_cold = 100.0
+    solver.T_hot = 500.0
+    solver.u = solver.set_initial_condition()
+
+    assert hasattr(solver, "u")
+    assert solver.u is not None
+
+
+# class TestDiffusion2D(unittest.TestCase):
+#     def setUp(self):
+#         self.solver = SolveDiffusion2D()
+
+#     def test_initialize_domain(self):
+#         self.solver.initialize_domain(w=20.0, h=10.0, dx=0.5, dy=0.5)
+#         self.assertEqual(self.solver.nx, 40)
+#         self.assertEqual(self.solver.ny, 20)
+
+#     def test_initialize_physical_parameters(self):
+#         self.solver.dx = 0.5
+#         self.solver.dy = 0.5
+#         self.solver.initialize_physical_parameters(d=4.0, T_cold=100.0, T_hot=500.0)
+
+#         expected = (
+#             (self.solver.dx * self.solver.dx) * (self.solver.dy * self.solver.dy)
+#             / (2 * 4.0 * (self.solver.dx * self.solver.dx + self.solver.dy * self.solver.dy))
+#         )
+
+#         self.assertAlmostEqual(self.solver.dt, expected, places=7)
+
+#     def test_set_initial_condition(self):
+#         self.solver.nx = 40
+#         self.solver.ny = 20
+#         self.solver.dx = 0.5
+#         self.solver.dy = 0.5
+#         self.solver.T_cold = 100.0
+#         self.solver.T_hot = 500.0
+#         self.solver.u = self.solver.set_initial_condition()
+
+#         self.assertIsNotNone(self.solver.u)
diff --git a/tox.toml b/tox.toml
new file mode 100644
index 0000000..2f08654
--- /dev/null
+++ b/tox.toml
@@ -0,0 +1,26 @@
+[tox]
+env_list = ["unit", "integration", "report"]
+skipsdist = true
+
+[testenv]
+deps = ["-rrequirements.txt"]
+
+["testenv:unit"]
+commands = [
+    ["coverage", "run", "-a", "-m", "unittest", "discover", "tests/unit"]
+]
+set_env = { PYTHONPATH = "{toxinidir}" }
+
+["testenv:integration"]
+commands = [
+    ["coverage", "run", "-a", "-m", "pytest", "tests/integration"]
+]
+set_env = { PYTHONPATH = "{toxinidir}" }
+
+["testenv:report"]
+deps = ["coverage"]
+skip_install = true
+commands = [
+    ["coverage", "report"],
+    ["coverage", "html"]
+]
\ No newline at end of file