diff --git a/README.md b/README.md
index fe54835..056000e 100644
--- a/README.md
+++ b/README.md
@@ -7,13 +7,28 @@ Please follow the instructions in [pypi_exercise.md](https://github.com/Simulati
 The code used in this exercise is based on [Chapter 7 of the book "Learning Scientific Programming with Python"](https://scipython.com/book/chapter-7-matplotlib/examples/the-two-dimensional-diffusion-equation/).
 
 ## Project description
+This code solves the diffusion equation in 2D over a square domain which is at a certain temperature and a circular disc at the center which is at a higher temperature. This code solves the diffusion equation using the Finite Difference Method. The thermal diffusivity and initial conditions of the system can be changed by the user. The code produces four plots at various timepoints of the simulation. The diffusion process can be clearly observed in these plots.
 
 ## Installing the package
 
 ### Using pip3 to install from PyPI
 
+```bash
+pip install -i https://test.pypi.org/simple/ luosg-diffusion2d==0.0.1
+```
+
 ### Required dependencies
 
+- **Python** >= 3.6
+- **NumPy**
+- **Matplotlib**
+
 ## Running this package
+Following Commands could be run in python Shell:
+```python
+>>> from luosg_diffusion2d import diffusion2d
+>>> diffusion2d.solve()
+```
 
 ## Citing
+Instructions in [pypi_exercise.md](https://github.com/Simulation-Software-Engineering/Lecture-Material/blob/main/03_building_and_packaging/pypi_exercise.md).
\ No newline at end of file
diff --git a/diff_plot_1.png b/diff_plot_1.png
new file mode 100644
index 0000000..67145b1
Binary files /dev/null and b/diff_plot_1.png differ
diff --git a/diff_plot_2.png b/diff_plot_2.png
new file mode 100644
index 0000000..67145b1
Binary files /dev/null and b/diff_plot_2.png differ
diff --git a/dist/luosg_diffusion2d-0.0.1-py3-none-any.whl b/dist/luosg_diffusion2d-0.0.1-py3-none-any.whl
new file mode 100644
index 0000000..9ee3dc1
Binary files /dev/null and b/dist/luosg_diffusion2d-0.0.1-py3-none-any.whl differ
diff --git a/dist/luosg_diffusion2d-0.0.1.tar.gz b/dist/luosg_diffusion2d-0.0.1.tar.gz
new file mode 100644
index 0000000..3e1c4b1
Binary files /dev/null and b/dist/luosg_diffusion2d-0.0.1.tar.gz differ
diff --git a/luosg_diffusion2d/__init__.py b/luosg_diffusion2d/__init__.py
new file mode 100644
index 0000000..e69de29
diff --git a/luosg_diffusion2d/__pycache__/__init__.cpython-38.pyc b/luosg_diffusion2d/__pycache__/__init__.cpython-38.pyc
new file mode 100644
index 0000000..19bb453
Binary files /dev/null and b/luosg_diffusion2d/__pycache__/__init__.cpython-38.pyc differ
diff --git a/luosg_diffusion2d/__pycache__/diffusion2d.cpython-38.pyc b/luosg_diffusion2d/__pycache__/diffusion2d.cpython-38.pyc
new file mode 100644
index 0000000..c4cd8bf
Binary files /dev/null and b/luosg_diffusion2d/__pycache__/diffusion2d.cpython-38.pyc differ
diff --git a/luosg_diffusion2d/__pycache__/output.cpython-38.pyc b/luosg_diffusion2d/__pycache__/output.cpython-38.pyc
new file mode 100644
index 0000000..8a6f9e7
Binary files /dev/null and b/luosg_diffusion2d/__pycache__/output.cpython-38.pyc differ
diff --git a/luosg_diffusion2d/diffusion2d.py b/luosg_diffusion2d/diffusion2d.py
new file mode 100644
index 0000000..ff99918
--- /dev/null
+++ b/luosg_diffusion2d/diffusion2d.py
@@ -0,0 +1,73 @@
+"""
+Solving the two-dimensional diffusion equation
+
+Example acquired from https://scipython.com/book/chapter-7-matplotlib/examples/the-two-dimensional-diffusion-equation/
+"""
+
+import numpy as np
+import matplotlib.pyplot as plt
+from luosg_diffusion2d.output import create_plot, output_plots
+
+
+def solve(dx = 0.1, dy = 0.1, D = 4):
+    # plate size, mm
+    w = h = 10.
+
+
+    # Initial cold temperature of square domain
+    T_cold = 300
+
+    # Initial hot temperature of circular disc at the center
+    T_hot = 700
+
+    # Number of discrete mesh points in X and Y directions
+    nx, ny = int(w / dx), int(h / dy)
+
+    # Computing a stable time step
+    dx2, dy2 = dx * dx, dy * dy
+    dt = dx2 * dy2 / (2 * D * (dx2 + dy2))
+
+    print("dt = {}".format(dt))
+
+    u0 = T_cold * np.ones((nx, ny))
+    u = u0.copy()
+
+    # Initial conditions - circle of radius r centred at (cx,cy) (mm
+    r = min(h, w) / 4.0
+    cx = w / 2.0
+    cy = h / 2.0
+    r2 = r ** 2
+    for i in range(nx):
+        for j in range(ny):
+            p2 = (i * dx - cx) ** 2 + (j * dy - cy) ** 2
+            if p2 < r2:
+                u0[i, j] = T_hot
+
+
+    def do_timestep(u_nm1, u, D, dt, dx2, dy2):
+        # Propagate with forward-difference in time, central-difference in space
+        u[1:-1, 1:-1] = u_nm1[1:-1, 1:-1] + D * dt * (
+            (u_nm1[2:, 1:-1] - 2 * u_nm1[1:-1, 1:-1] + u_nm1[:-2, 1:-1]) / dx2
+            + (u_nm1[1:-1, 2:] - 2 * u_nm1[1:-1, 1:-1] + u_nm1[1:-1, :-2]) / dy2)
+
+        u_nm1 = u.copy()
+        return u_nm1, u
+
+
+    # Number of timesteps
+    nsteps = 101
+    # Output 4 figures at these timesteps
+    n_output = [0, 10, 50, 100]
+    fig_counter = 0
+    fig = plt.figure()
+
+    # Time loop
+    for n in range(nsteps):
+        u0, u = do_timestep(u0, u, D, dt, dx2, dy2)
+
+        # Create figure
+        if n in n_output:
+            im, fig_counter = create_plot(T_cold, T_hot, dt, u, fig_counter, fig, n)
+
+    # Plot output figures
+    output_plots(fig, im)
\ No newline at end of file
diff --git a/luosg_diffusion2d/output.py b/luosg_diffusion2d/output.py
new file mode 100644
index 0000000..571c160
--- /dev/null
+++ b/luosg_diffusion2d/output.py
@@ -0,0 +1,16 @@
+import matplotlib.pyplot as plt
+
+def create_plot(T_cold, T_hot, dt, u, fig_counter, fig, n): 
+    fig_counter += 1
+    ax = fig.add_subplot(220 + fig_counter)
+    im = ax.imshow(u.copy(), cmap=plt.get_cmap('hot'), vmin=T_cold, vmax=T_hot)  # image for color bar axes
+    ax.set_axis_off()
+    ax.set_title('{:.1f} ms'.format(n * dt * 1000))
+    return im, fig_counter
+
+def output_plots(fig, im):
+    fig.subplots_adjust(right=0.85)
+    cbar_ax = fig.add_axes([0.9, 0.15, 0.03, 0.7])
+    cbar_ax.set_xlabel('$T$ / K', labelpad=20)
+    fig.colorbar(im, cax=cbar_ax)
+    plt.show()
\ No newline at end of file
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000..3d4f7b2
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,28 @@
+[build-system]
+requires = ["setuptools", "wheel"]
+
+[project]
+name = "luosg_diffusion2d"
+version = "0.0.1"
+authors = [
+	{name="Luo, Shiting", email="schitingluo@gmail.com"}
+]
+description = "2D diffusion equation simulation"
+readme = "README.md"
+requires-python = ">=3.6"
+license = { file = "LICENSE" }
+keywords = ["sse", "2D diffusion"]
+classifiers = [
+    "Programming Language :: Python :: 3",
+]
+dependencies = [
+    "numpy",
+    "matplotlib",
+]
+
+[project.urls]
+Homepage = "https://github.com/Simulation-Software-Engineering/diffusion2D"
+
+
+[project.entry-points."simulation"]
+solve = "luosg_diffusion2d:solve"
\ No newline at end of file
diff --git a/setup.py b/setup.py
new file mode 100644
index 0000000..e342256
--- /dev/null
+++ b/setup.py
@@ -0,0 +1,4 @@
+from setuptools import setup
+
+if __name__ == "__main__":
+    setup()
\ No newline at end of file
diff --git a/test.py b/test.py
new file mode 100644
index 0000000..45974a6
--- /dev/null
+++ b/test.py
@@ -0,0 +1,3 @@
+from luosg_diffusion2d.diffusion2d import solve
+
+solve(dx = 0.1, dy = 0.1, D = 4)
\ No newline at end of file