diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..5631645
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,5 @@
+__pycache__
+.pytest__cache
+/htmlcov/
+.coverage
+.tox/
\ No newline at end of file
diff --git a/coverage-report.pdf b/coverage-report.pdf
new file mode 100644
index 0000000..570abfe
Binary files /dev/null and b/coverage-report.pdf differ
diff --git a/diffusion2d.py b/diffusion2d.py
index 51a07f2..3dd21d5 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 float"
+        assert isinstance(h, float), "h must be float"
+        assert isinstance(dx, float), "dx must be float"
+        assert isinstance(dy, float), "dy must be 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 float"
+        assert isinstance(T_cold, float), "T_cold must be float"
+        assert isinstance(T_hot, float), "T_hot must be float"
+
         self.D = d
         self.T_cold = T_cold
         self.T_hot = T_hot
@@ -53,7 +62,7 @@ def initialize_physical_parameters(self, d=4., T_cold=300, T_hot=700):
         # Computing a stable time step
         dx2, dy2 = self.dx * self.dx, self.dy * self.dy
         self.dt = dx2 * dy2 / (2 * self.D * (dx2 + dy2))
-
+        
         print("dt = {}".format(self.dt))
 
     def set_initial_condition(self):
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000..6cdf7cd
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,3 @@
+numpy
+matplotlib
+pytest
\ No newline at end of file
diff --git a/shell.nix b/shell.nix
new file mode 100644
index 0000000..70bdbdf
--- /dev/null
+++ b/shell.nix
@@ -0,0 +1,12 @@
+{ pkgs ? import <nixpkgs> {} }:
+pkgs.mkShell {
+  packages = with pkgs; [
+    python313
+    python313Packages.tox
+    stdenv.cc.cc.lib  # provides libstdc++.so.6
+  ];
+
+  shellHook = ''
+    export LD_LIBRARY_PATH=${pkgs.stdenv.cc.cc.lib}/lib:$LD_LIBRARY_PATH
+  '';
+}
\ No newline at end of file
diff --git a/tests/integration/test_diffusion2d.py b/tests/integration/test_diffusion2d.py
index fd026b4..164100d 100644
--- a/tests/integration/test_diffusion2d.py
+++ b/tests/integration/test_diffusion2d.py
@@ -2,18 +2,60 @@
 Tests for functionality checks in class SolveDiffusion2D
 """
 
+import pytest
+import numpy as np
 from diffusion2d import SolveDiffusion2D
 
+@pytest.fixture
+def solver():
+    solver = SolveDiffusion2D()
+    return solver
 
-def test_initialize_physical_parameters():
+def test_initialize_physical_parameters(solver):
     """
     Checks function SolveDiffusion2D.initialize_domain
     """
-    solver = SolveDiffusion2D()
+    w = 20.
+    h = 15.
+    dx = 0.2
+    dy = 0.3
+    d = 5.
+    T_cold = 200.
+    T_hot = 500.
+
+    solver.initialize_domain(w,h,dx,dy)
+    solver.initialize_physical_parameters(d,T_cold,T_hot)
 
+    expected_dt = 0.00276923076923077
 
-def test_set_initial_condition():
+    assert solver.dt == pytest.approx(expected_dt), "dt should be approximately {expected_dt}"
+
+
+def test_set_initial_condition(solver):
     """
     Checks function SolveDiffusion2D.get_initial_function
     """
-    solver = SolveDiffusion2D()
+    w = 20.
+    h = 15.
+    dx = 0.2
+    dy = 0.3
+    d = 5.
+    T_cold = 200.
+    T_hot = 500.
+
+    solver.initialize_domain(w,h,dx,dy)
+    solver.initialize_physical_parameters(d,T_cold,T_hot)
+    computed_u = solver.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
+
+    assert np.array_equal(expected_u, computed_u), "Initial condition array is not equal to the expected"
+
diff --git a/tests/unit/test_diffusion2d_functions.py b/tests/unit/test_diffusion2d_functions.py
index c4277ff..e0ebdab 100644
--- a/tests/unit/test_diffusion2d_functions.py
+++ b/tests/unit/test_diffusion2d_functions.py
@@ -2,6 +2,8 @@
 Tests for functions in class SolveDiffusion2D
 """
 
+import numpy as np
+import pytest
 from diffusion2d import SolveDiffusion2D
 
 
@@ -9,7 +11,17 @@ def test_initialize_domain():
     """
     Check function SolveDiffusion2D.initialize_domain
     """
+    w = 20.
+    h = 15.
+    dx = 0.2
+    dy = 0.3
+    expected_nx  = 100
+    expected_ny = 50
     solver = SolveDiffusion2D()
+    solver.initialize_domain(w,h,dx,dy)
+    assert solver.nx == pytest.approx(expected_nx), "nx should be 100"
+    assert solver.ny == pytest.approx(expected_ny), "ny should be 50"
+
 
 
 def test_initialize_physical_parameters():
@@ -17,10 +29,41 @@ def test_initialize_physical_parameters():
     Checks function SolveDiffusion2D.initialize_domain
     """
     solver = SolveDiffusion2D()
+    solver.dx = 0.2
+    solver.dy = 0.3
+    d = 5.
+    T_cold = 200.
+    T_hot =  500.
 
+    expected_dt = 0.00276923076923077
+    solver.initialize_physical_parameters(d, T_cold, T_hot)
+    assert solver.dt == pytest.approx(expected_dt), "dt should be approximately {expected_dt}"
+    assert solver.T_cold == pytest.approx(T_cold), "T_cold should be 200"
+    assert solver.T_hot == pytest.approx(T_hot), "T_hot should be 500"
 
 def test_set_initial_condition():
     """
     Checks function SolveDiffusion2D.get_initial_function
     """
     solver = SolveDiffusion2D()
+    solver.dx = 0.2
+    solver.dy = 0.3
+    solver.nx = 100
+    solver.ny = 50
+    solver.T_cold = 200.
+    solver.T_hot = 500.
+    solver.dt = 0.00276923076923077
+    solver.D = 5.
+    
+    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
+
+    computed_u = solver.set_initial_condition()
+    assert np.array_equal(expected_u, computed_u), "Initial condition array is not equal to the expected"
diff --git a/tox.toml b/tox.toml
new file mode 100644
index 0000000..8c95dcc
--- /dev/null
+++ b/tox.toml
@@ -0,0 +1,7 @@
+# tox.toml
+requires = ["tox>=4.0"]
+env_list = ["testing"]
+
+[env.testing]
+deps = ["-rrequirements.txt"]
+commands = [["python", "-m", "pytest"]]
\ No newline at end of file