diff --git a/.gitignore b/.gitignore
new file mode 100644
index 0000000..2f1365b
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,22 @@
+cat > .gitignore << 'EOF'
+# Python
+__pycache__/
+*.py[cod]
+*$py.class
+*.so
+.Python
+venv/
+env/
+.tox/
+.pytest_cache/
+.coverage
+htmlcov/
+
+# OS
+.DS_Store
+Thumbs.db
+
+# IDE
+.vscode/
+.idea/
+EOF
diff --git a/coverage-report.pdf b/coverage-report.pdf
new file mode 100644
index 0000000..9f08855
Binary files /dev/null and b/coverage-report.pdf differ
diff --git a/diffusion2d.py b/diffusion2d.py
index 51a07f2..f5c65df 100644
--- a/diffusion2d.py
+++ b/diffusion2d.py
@@ -38,6 +38,12 @@ 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"
+        assert w > 0 and h > 0, "Dimensions must be positive"
+        assert dx > 0 and dy > 0, "Grid spacings must be positive"
         self.w = w
         self.h = h
         self.dx = dx
@@ -45,7 +51,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.0, T_hot=700.0):
+        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"
+        assert d > 0, "Diffusivity must be positive"
         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..2d0ef4d
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1,4 @@
+numpy>=1.24
+matplotlib>=3.7
+pytest>=7.0
+coverage>=6.0
diff --git a/tests/integration/test_diffusion2d.py b/tests/integration/test_diffusion2d.py
index fd026b4..3f78e94 100644
--- a/tests/integration/test_diffusion2d.py
+++ b/tests/integration/test_diffusion2d.py
@@ -1,19 +1,59 @@
 """
 Tests for functionality checks in class SolveDiffusion2D
 """
-
+import numpy as np
+import sys
+import os
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), '..')))  
 from diffusion2d import SolveDiffusion2D
 
-
 def test_initialize_physical_parameters():
     """
-    Checks function SolveDiffusion2D.initialize_domain
+    Checks function SolveDiffusion2D.initialize_physical_parameters 
+    Integration test: domain setup affects dt calculation
     """
     solver = SolveDiffusion2D()
-
+    
+    # Non-default values for integration test
+    w, h, dx, dy = 12.0, 8.0, 0.2, 0.1
+    D, T_cold, T_hot = 3.5, 280.0, 750.0
+    
+    solver.initialize_domain(w, h, dx, dy)
+    solver.initialize_physical_parameters(D, T_cold, T_hot)
+    
+    # Manual calculation of expected dt (depends on domain dx,dy)
+    dx2, dy2 = dx**2, dy**2
+    expected_dt = dx2 * dy2 / (2 * D * (dx2 + dy2))
+    
+    assert abs(solver.dt - expected_dt) < 1e-12
+    assert solver.D == D
+    assert solver.T_cold == T_cold
+    assert solver.T_hot == T_hot
 
 def test_set_initial_condition():
     """
-    Checks function SolveDiffusion2D.get_initial_function
+    Checks function SolveDiffusion2D.set_initial_condition
+    Integration test: domain + physics → correct initial field u
     """
     solver = SolveDiffusion2D()
+    
+    # Default-like params to match code's circle
+    solver.initialize_domain()
+    solver.initialize_physical_parameters()
+    
+    u = solver.set_initial_condition()
+    
+    # Replicate exact circle logic from set_initial_condition
+    r, cx, cy = 2, 5, 5
+    r2 = r**2
+    expected_u = solver.T_cold * np.ones((solver.nx, solver.ny))
+    
+    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
+    
+    # Exact match (no floating point issues)
+    assert np.array_equal(u, expected_u)
+    assert u.shape == (solver.nx, solver.ny)
diff --git a/tests/unit/test_diffusion2d_functions.py b/tests/unit/test_diffusion2d_functions.py
index c4277ff..792043a 100644
--- a/tests/unit/test_diffusion2d_functions.py
+++ b/tests/unit/test_diffusion2d_functions.py
@@ -1,26 +1,58 @@
-"""
-Tests for functions in class SolveDiffusion2D
-"""
+import pytest
+import sys
+import os
+import numpy as np
 
+sys.path.insert(0, os.path.abspath(os.path.join(os.path.dirname(__file__), "..")))
 from diffusion2d import SolveDiffusion2D
 
 
 def test_initialize_domain():
-    """
-    Check function SolveDiffusion2D.initialize_domain
-    """
     solver = SolveDiffusion2D()
+    w, h, dx, dy = 20.0, 15.0, 0.5, 0.3  # Non-default, w != h to catch bugs
+    expected_nx = int(w / dx)  # 40
+    expected_ny = int(h / dy)  # 50
+
+    solver.initialize_domain(w, h, dx, dy)
+
+    assert solver.nx == expected_nx
+    assert solver.ny == expected_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
-    """
     solver = SolveDiffusion2D()
+    solver.initialize_domain()  # Only for dx,dy (minimal dep, but isolated check)
+    d, T_cold, T_hot = 5.0, 250.0, 800.0
+    dx2, dy2 = solver.dx**2, solver.dy**2
+    expected_dt = dx2 * dy2 / (2 * d * (dx2 + dy2))  # ~6.25e-4
+
+    solver.initialize_physical_parameters(d, T_cold, T_hot)
+
+    assert solver.D == d
+    assert solver.T_cold == T_cold
+    assert solver.T_hot == T_hot
+    assert abs(solver.dt - expected_dt) < 1e-10  # Floating-point tolerance
 
 
 def test_set_initial_condition():
-    """
-    Checks function SolveDiffusion2D.get_initial_function
-    """
     solver = SolveDiffusion2D()
+    solver.initialize_domain(w=10.0, h=10.0, dx=0.1, dy=0.1)
+    solver.initialize_physical_parameters()
+
+    u = solver.set_initial_condition()
+
+    print("u shape:", u.shape)  # Should be (100,100)
+    print("T_cold:", solver.T_cold, type(solver.T_cold))
+    print("T_hot:", solver.T_hot, type(solver.T_hot))
+    center_i = int(5.0 / solver.dx)
+    center_j = int(5.0 / solver.dy)
+    print("center_i,j:", center_i, center_j)
+    print("u[center]:", u[center_i, center_j])
+    print("Hot count:", np.sum(u == solver.T_hot))
+    print("Sample u[49:51,49:51]:\n", u[49:51, 49:51])
+
+    # Keep existing asserts for now
diff --git a/tox.toml b/tox.toml
new file mode 100644
index 0000000..666a0f4
--- /dev/null
+++ b/tox.toml
@@ -0,0 +1,17 @@
+[tool.tox]
+envlist = ["py"]
+
+[tool.tox.env_run_base]
+deps = [
+    "-rrequirements.txt"
+]
+
+commands = [
+    "echo '=== PYTEST (verbose) ===",
+    "pytest tests/ -v -s --tb=short",
+    "echo '=== COVERAGE REPORT ===",
+    "coverage run --source=diffusion2d -m pytest tests/",
+    "coverage report",
+    "echo '=== HTML in htmlcov/ ===",
+    "coverage html"
+]
\ No newline at end of file