diff --git a/.gitignore b/.gitignore
new file mode 100644
index 00000000..f48d3305
--- /dev/null
+++ b/.gitignore
@@ -0,0 +1,209 @@
+# Created by https://www.toptal.com/developers/gitignore/api/macos,python
+# Edit at https://www.toptal.com/developers/gitignore?templates=macos,python
+
+### macOS ###
+# General
+.DS_Store
+.AppleDouble
+.LSOverride
+
+# Icon must end with two \r
+Icon
+
+
+# Thumbnails
+._*
+
+# Files that might appear in the root of a volume
+.DocumentRevisions-V100
+.fseventsd
+.Spotlight-V100
+.TemporaryItems
+.Trashes
+.VolumeIcon.icns
+.com.apple.timemachine.donotpresent
+
+# Directories potentially created on remote AFP share
+.AppleDB
+.AppleDesktop
+Network Trash Folder
+Temporary Items
+.apdisk
+
+### macOS Patch ###
+# iCloud generated files
+*.icloud
+
+### Python ###
+# Byte-compiled / optimized / DLL files
+__pycache__/
+*.py[cod]
+*$py.class
+
+# C extensions
+*.so
+
+# Distribution / packaging
+.Python
+build/
+develop-eggs/
+dist/
+downloads/
+eggs/
+.eggs/
+lib/
+lib64/
+parts/
+sdist/
+var/
+wheels/
+share/python-wheels/
+*.egg-info/
+.installed.cfg
+*.egg
+MANIFEST
+
+# PyInstaller
+#  Usually these files are written by a python script from a template
+#  before PyInstaller builds the exe, so as to inject date/other infos into it.
+*.manifest
+*.spec
+
+# Installer logs
+pip-log.txt
+pip-delete-this-directory.txt
+
+# Unit test / coverage reports
+htmlcov/
+.tox/
+.nox/
+.coverage
+.coverage.*
+.cache
+nosetests.xml
+coverage.xml
+*.cover
+*.py,cover
+.hypothesis/
+.pytest_cache/
+cover/
+
+# Translations
+*.mo
+*.pot
+
+# Django stuff:
+*.log
+local_settings.py
+db.sqlite3
+db.sqlite3-journal
+
+# Flask stuff:
+instance/
+.webassets-cache
+
+# Scrapy stuff:
+.scrapy
+
+# Sphinx documentation
+docs/_build/
+
+# PyBuilder
+.pybuilder/
+target/
+
+# Jupyter Notebook
+.ipynb_checkpoints
+
+# IPython
+profile_default/
+ipython_config.py
+
+# pyenv
+#   For a library or package, you might want to ignore these files since the code is
+#   intended to run in multiple environments; otherwise, check them in:
+# .python-version
+
+# pipenv
+#   According to pypa/pipenv#598, it is recommended to include Pipfile.lock in version control.
+#   However, in case of collaboration, if having platform-specific dependencies or dependencies
+#   having no cross-platform support, pipenv may install dependencies that don't work, or not
+#   install all needed dependencies.
+#Pipfile.lock
+
+# poetry
+#   Similar to Pipfile.lock, it is generally recommended to include poetry.lock in version control.
+#   This is especially recommended for binary packages to ensure reproducibility, and is more
+#   commonly ignored for libraries.
+#   https://python-poetry.org/docs/basic-usage/#commit-your-poetrylock-file-to-version-control
+#poetry.lock
+
+# pdm
+#   Similar to Pipfile.lock, it is generally recommended to include pdm.lock in version control.
+#pdm.lock
+#   pdm stores project-wide configurations in .pdm.toml, but it is recommended to not include it
+#   in version control.
+#   https://pdm.fming.dev/#use-with-ide
+.pdm.toml
+
+# PEP 582; used by e.g. github.com/David-OConnor/pyflow and github.com/pdm-project/pdm
+__pypackages__/
+
+# Celery stuff
+celerybeat-schedule
+celerybeat.pid
+
+# SageMath parsed files
+*.sage.py
+
+# Environments
+.env
+.venv
+env/
+venv/
+ENV/
+env.bak/
+venv.bak/
+
+# Spyder project settings
+.spyderproject
+.spyproject
+
+# Rope project settings
+.ropeproject
+
+# mkdocs documentation
+/site
+
+# mypy
+.mypy_cache/
+.dmypy.json
+dmypy.json
+
+# Pyre type checker
+.pyre/
+
+# pytype static type analyzer
+.pytype/
+
+# Cython debug symbols
+cython_debug/
+
+# PyCharm
+#  JetBrains specific template is maintained in a separate JetBrains.gitignore that can
+#  be found at https://github.com/github/gitignore/blob/main/Global/JetBrains.gitignore
+#  and can be added to the global gitignore or merged into this file.  For a more nuclear
+#  option (not recommended) you can uncomment the following to ignore the entire idea folder.
+#.idea/
+
+### Python Patch ###
+# Poetry local configuration file - https://python-poetry.org/docs/configuration/#local-configuration
+poetry.toml
+
+# ruff
+.ruff_cache/
+
+# LSP config files
+pyrightconfig.json
+
+# End of https://www.toptal.com/developers/gitignore/api/macos,python
\ No newline at end of file
diff --git a/README.md b/README.md
index fe548354..a8ae5fdb 100644
--- a/README.md
+++ b/README.md
@@ -8,12 +8,39 @@ The code used in this exercise is based on [Chapter 7 of the book "Learning Scie
 
 ## 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
+pip3 install chenla_diffusion2D
+```
+
 ### Required dependencies
 
+* Check if your system has Python version >= 3.6 and update it if it is older than 3.6.
+* Install pip, build, and Twine.
+* Install NumPy and Matplotlib with pip.
+
 ## Running this package
 
+You could run the package via `diffusion2d.solve()` after import. This will generate 4 plots with the default parameters `dx=0.1, dy=0.1, D=4.`, where `dx` and `dy` are intervals in x-, y- directions (mm), and `D` is the thermal diffusivity of steel (mm^2/s).
+
+```python
+from chenla_diffusion2D import diffusion2d
+
+diffusion2d.solve()
+```
+
+You could also assign the values of `dx`, `dy` and `D` by yourself.
+
+```python
+from chenla_diffusion2D import diffusion2d
+
+diffusion2d.solve(dx=0.5, dy=0.5, D=16)
+```
+
 ## Citing
diff --git a/__init__.py b/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/chenla_diffusion2D/__init__.py b/chenla_diffusion2D/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/chenla_diffusion2D/diffusion2d.py b/chenla_diffusion2D/diffusion2d.py
new file mode 100644
index 00000000..ec8cc0f7
--- /dev/null
+++ b/chenla_diffusion2D/diffusion2d.py
@@ -0,0 +1,77 @@
+"""
+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 .output import create_plot, output_plots
+
+def solve(dx=0.1, dy=0.1, D=4.):
+    """
+    dx, dy: intervals in x-, y- directions, mm
+    D: Thermal diffusivity of steel, mm^2/s
+    """
+    # 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:
+            fig_counter += 1
+            im = create_plot(fig, u.copy(), n*dt*1000, 220+fig_counter, T_cold, T_hot)
+
+    # Plot output figures
+    output_plots(fig, im)
diff --git a/chenla_diffusion2D/output.py b/chenla_diffusion2D/output.py
new file mode 100644
index 00000000..c2033e3e
--- /dev/null
+++ b/chenla_diffusion2D/output.py
@@ -0,0 +1,20 @@
+import numpy as np
+import matplotlib.pyplot as plt
+
+def create_plot(fig, data, timestamp, subplot_pos, vmin, vmax):
+    """This is the function to create one plot for a particular time stamp.
+    """
+    ax = fig.add_subplot(subplot_pos)
+    im = ax.imshow(data, cmap=plt.get_cmap('hot'), vmin=vmin, vmax=vmax)  # image for color bar axes
+    ax.set_axis_off()
+    ax.set_title('{:.1f} ms'.format(timestamp))
+    return im
+
+def output_plots(fig, im):
+    """This is the function to output all the four plots as one figure.
+    """
+    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/diffusion2d.py b/diffusion2d.py
deleted file mode 100644
index c0c6083a..00000000
--- a/diffusion2d.py
+++ /dev/null
@@ -1,81 +0,0 @@
-"""
-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
-
-# plate size, mm
-w = h = 10.
-# intervals in x-, y- directions, mm
-dx = dy = 0.1
-# Thermal diffusivity of steel, mm^2/s
-D = 4.
-
-# 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:
-        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))
-
-# Plot output figures
-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()
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 00000000..f656d835
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,34 @@
+[build-system]
+requires = ["setuptools"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "chenla_diffusion2D"
+authors = [
+    {name = "LingChia Chen", email = "st191429@stud.uni-stuttgart.de"},
+]
+description = "This is a package to show plots of the solutions of a diffusion2D problem."
+readme = "README.md"
+requires-python = ">=3.6"
+keywords = ["diffusion2D", "plots"]
+license = {text = "BSD-3-Clause"}
+classifiers = [
+    "License :: OSI Approved :: BSD License",
+    "Development Status :: 2 - Pre-Alpha",
+    "Programming Language :: Python :: 3",
+]
+dependencies = [
+    "numpy",
+    "matplotlib",
+]
+version = "0.0.3"
+
+[project.urls]
+Homepage = "https://github.com/lingachen/diffusion2D/tree/chenla-dev"
+
+[project.optional-dependencies]
+pdf = ["ReportLab>=1.2", "RXP"]
+rest = ["docutils>=0.3", "pack ==1.1, ==1.3"]
+
+[project.scripts]
+chenla_diffusion2d = "chenla_diffusion2D.diffusion2d:solve"
\ No newline at end of file