diff --git a/coverage-report.pdf b/coverage-report.pdf
new file mode 100644
index 0000000..e64c4c0
Binary files /dev/null and b/coverage-report.pdf differ
diff --git a/diffusion2d.py b/diffusion2d.py
index 51a07f2..17bc002 100644
--- a/diffusion2d.py
+++ b/diffusion2d.py
@@ -38,6 +38,10 @@ 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 +49,10 @@ 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..021aae7
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,4 @@
+numpy
+matplotlib
+pytest
+coverage
diff --git a/tests/integration/test_diffusion2d.py b/tests/integration/test_diffusion2d.py
index fd026b4..2780559 100644
--- a/tests/integration/test_diffusion2d.py
+++ b/tests/integration/test_diffusion2d.py
@@ -2,6 +2,7 @@
 Tests for functionality checks in class SolveDiffusion2D
 """
 
+import numpy as np
 from diffusion2d import SolveDiffusion2D
 
 
@@ -10,6 +11,35 @@ def test_initialize_physical_parameters():
     Checks function SolveDiffusion2D.initialize_domain
     """
     solver = SolveDiffusion2D()
+    
+    # Initialize domain with specific values
+    w = 12.0
+    h = 12.0
+    dx = 0.4
+    dy = 0.4
+    solver.initialize_domain(w, h, dx, dy)
+    
+    # Initialize physical parameters
+    d = 3.0
+    T_cold = 250.0
+    T_hot = 750.0
+    solver.initialize_physical_parameters(d, T_cold, T_hot)
+    
+    # Manually calculate expected dt
+    dx2 = dx * dx
+    dy2 = dy * dy
+    expected_dt = dx2 * dy2 / (2 * d * (dx2 + dy2))
+    
+    # Verify the dt calculation
+    assert abs(solver.dt - expected_dt) < 1e-10, \
+        f"Expected dt={expected_dt}, but got {solver.dt}"
+    
+    # Also verify that all parameters are set correctly
+    assert solver.D == d
+    assert solver.T_cold == T_cold
+    assert solver.T_hot == T_hot
+    assert solver.dx == dx
+    assert solver.dy == dy
 
 
 def test_set_initial_condition():
@@ -17,3 +47,43 @@ def test_set_initial_condition():
     Checks function SolveDiffusion2D.get_initial_function
     """
     solver = SolveDiffusion2D()
+    w = 10.0
+    h = 10.0
+    dx = 0.5
+    dy = 0.5
+    solver.initialize_domain(w, h, dx, dy)
+    
+    # Initialize physical parameters
+    d = 4.0
+    T_cold = 300.0
+    T_hot = 700.0
+    solver.initialize_physical_parameters(d, T_cold, T_hot)
+    
+    # Get initial condition
+    u = solver.set_initial_condition()
+    
+    # Verify shape matches domain
+    expected_shape = (solver.nx, solver.ny)
+    assert u.shape == expected_shape, \
+        f"Expected shape {expected_shape}, but got {u.shape}"
+    
+    # Manually compute expected initial condition
+    # Most values should be T_cold
+    u_expected = T_cold * np.ones((solver.nx, solver.ny))
+    
+    # Circular disc at center with T_hot
+    r, cx, cy = 2, 5, 5
+    r2 = r ** 2
+    for i in range(solver.nx):
+        for j in range(solver.ny):
+            p2 = (i * dx - cx) ** 2 + (j * dy - cy) ** 2
+            if p2 < r2:
+                u_expected[i, j] = T_hot
+    
+    # Check that the computed u matches expected
+    assert np.allclose(u, u_expected), \
+        "Initial condition array does not match expected values"
+    
+    # Verify all values are within [T_cold, T_hot]
+    assert np.all(u >= T_cold), f"Some values are below T_cold={T_cold}"
+    assert np.all(u <= T_hot), f"Some values are above T_hot={T_hot}"
diff --git a/tests/unit/test_diffusion2d_functions.py b/tests/unit/test_diffusion2d_functions.py
index c4277ff..36c019f 100644
--- a/tests/unit/test_diffusion2d_functions.py
+++ b/tests/unit/test_diffusion2d_functions.py
@@ -2,6 +2,7 @@
 Tests for functions in class SolveDiffusion2D
 """
 
+import numpy as np
 from diffusion2d import SolveDiffusion2D
 
 
@@ -10,17 +11,97 @@ def test_initialize_domain():
     Check function SolveDiffusion2D.initialize_domain
     """
     solver = SolveDiffusion2D()
+    
+    # Use non-default values to test properly
+    w = 20.0
+    h = 30.0
+    dx = 0.5
+    dy = 0.2
+    
+    # Call the function
+    solver.initialize_domain(w, h, dx, dy)
+    
+    # Manually calculate expected values
+    expected_nx = int(w / dx)  # int(20.0 / 0.5) = 40
+    expected_ny = int(h / dy)  # int(30.0 / 0.2) = 150
+    
+    # Assert that the values match
+    assert solver.nx == expected_nx, f"Expected nx={expected_nx}, but got {solver.nx}"
+    assert solver.ny == expected_ny, f"Expected ny={expected_ny}, but got {solver.ny}"
+    assert solver.w == w
+    assert solver.h == h
+    assert solver.dx == dx
+    assert solver.dy == dy
 
 
 def test_initialize_physical_parameters():
     """
-    Checks function SolveDiffusion2D.initialize_domain
+    Checks function SolveDiffusion2D.initialize_physical_parameters
     """
+    # Create solver and mock the required member variables (instead of calling initialize_domain)
     solver = SolveDiffusion2D()
+    
+    # Manually set the required member variables
+    dx = 0.3
+    dy = 0.3
+    solver.dx = dx
+    solver.dy = dy
+    
+    # Test parameters
+    d = 2.5
+    T_cold = 250.0
+    T_hot = 750.0
+    
+    # Call the function
+    solver.initialize_physical_parameters(d, T_cold, T_hot)
+    
+    # Manually calculate expected dt
+    dx2 = dx * dx
+    dy2 = dy * dy
+    expected_dt = dx2 * dy2 / (2 * d * (dx2 + dy2))
+    
+    # Assert values
+    assert solver.D == d
+    assert solver.T_cold == T_cold
+    assert solver.T_hot == T_hot
+    assert abs(solver.dt - expected_dt) < 1e-10, f"Expected dt={expected_dt}, but got {solver.dt}"
 
 
 def test_set_initial_condition():
     """
-    Checks function SolveDiffusion2D.get_initial_function
+    Checks function SolveDiffusion2D.set_initial_condition
     """
+    # Create solver and manually set required member variables (mocking initialize_domain and initialize_physical_parameters)
     solver = SolveDiffusion2D()
+    
+    # Mock the required member variables
+    w = 10.0
+    h = 10.0
+    dx = 0.5
+    dy = 0.5
+    solver.w = w
+    solver.h = h
+    solver.dx = dx
+    solver.dy = dy
+    solver.nx = int(w / dx)
+    solver.ny = int(h / dy)
+    
+    T_cold = 300.0
+    T_hot = 700.0
+    solver.T_cold = T_cold
+    solver.T_hot = T_hot
+    
+    # Call the function
+    u = solver.set_initial_condition()
+    
+    # Check properties of the returned array
+    assert u.shape == (solver.nx, solver.ny), f"Expected shape {(solver.nx, solver.ny)}, got {u.shape}"
+    
+    # Check that most values are T_cold
+    assert np.all(u >= T_cold), f"Some values are below T_cold={T_cold}"
+    assert np.all(u <= T_hot), f"Some values are above T_hot={T_hot}"
+    
+    # Check that center has hot temperature (circular disc at center)
+    cx, cy = 5, 5
+    r = 2
+    assert u[int(cx/dx), int(cy/dy)] == T_hot, "Center should have T_hot"
diff --git a/tox.ini b/tox.ini
new file mode 100644
index 0000000..fe5ab38
--- /dev/null
+++ b/tox.ini
@@ -0,0 +1,13 @@
+[tox]
+envlist = py39, py310, py311, py312
+
+[testenv]
+deps = 
+    -r{toxinidir}/requirements.txt
+setenv =
+    PYTHONPATH = {toxinidir}
+commands =
+    pytest {posargs:tests/}
+
+[pytest]
+testpaths = tests