diff --git a/examples/AllenCahn/plot_4000.py b/examples/AllenCahn/plot_4000.py
index c73f95a..9d58eb4 100644
--- a/examples/AllenCahn/plot_4000.py
+++ b/examples/AllenCahn/plot_4000.py
@@ -4,7 +4,7 @@
 Free Energy: Allen-Cahn
 ^^^^^^^^^^^^^^^^^^^^^^^
 
-In this example, we study the free energy potential of the Allen-Cahn equation. The model consists of a bar of length $2a$ with boundary conditions at $x=-a$ of $\phi=-1$ and $\phi=1$ at $x=a$.
+In this example, we study the free energy potential of the Allen-Cahn equation. The model consists of a bar of length $2a$ with boundary conditions at $x=-a$ of $\phi=-1$ and $\phi=1$ at $x=a$. The theoric framework is explained in (:ref:`theory_allen_cahn`).
 
 .. code-block::
 
@@ -72,8 +72,21 @@
 ###############################################################################
 # Parameters definition
 # ---------------------
-# A length scale parameter of $l = 10.0$ is defined. The phase field is saved in a VTU file.
-# All the results are stored in a folder named results_folder_name = name.
+# First, we define an input class, which contains all the parameters needed for the setup 
+# and results of the simulation.
+#
+# The first term, $l$, specifies the length scale parameter for the problem, with $l = 10.0$.
+#
+# The next two options, `save_solution_xdmf` and `save_solution_vtu`, determine the file format
+# used to save the phase-field results (.xdmf or .vtu), which can then be visualized using
+# tools like ParaView or pvista. Both parameters are boolean values (`True` or `False`).
+# In this case, we set `save_solution_vtu=True` to save the phase-field results in .vtu format.
+#
+# Lastly, `results_folder_name` specifies the name of the folder where all results
+# and log information will be saved. If the folder does not exist, `phasefieldx` will create it. 
+# However, if the folder already exists, any previous data in it will be removed, and a new 
+# empty folder will be created in its place.
+#
 Data = Input(l=10.0,
              save_solution_xdmf=False,
              save_solution_vtu=True,
@@ -84,9 +97,12 @@
 # Mesh Definition
 # ---------------
 # A 2a mm x 1 mm rectangular mesh with quadrilateral elements is created using Dolfinx.
+# The mesh is a structured grid with quadrilateral elements:
+# - `divx`, `divy`: Number of elements along the x and y axes .
+# - `lx`, `ly`: Physical domain dimensions in x and y.
 a = 50.0
-divx, divy, divz = 100, 1, 1
-lx, ly, lz = 100.0, 1.0, 1.0
+divx, divy = 100, 1
+lx, ly = 100.0, 1.0
 
 msh = dolfinx.mesh.create_rectangle(mpi4py.MPI.COMM_WORLD,
                                     [np.array([-a, -0.5]),
@@ -96,22 +112,37 @@
 
 
 ###############################################################################
-# The left and right parts are identified and used to impose
-# the boundary conditions.
+# Boundary Identification
+# -----------------------
+# Boundary conditions and forces are applied to specific regions of the domain:
+# - `left`: Identifies the $y=-a$ boundary.
+# - `right`: Identifies the $y=a$ boundary.
+# `fdim` is the dimension of boundary facets (1D for a 2D mesh).
 def left(x):
     return np.isclose(x[0], -a)
 
-
 def right(x):
     return np.isclose(x[0], a)
 
-
 fdim = msh.topology.dim - 1
 
+# Using the `bottom` and `top` functions, we locate the facets on the left and right sides of the mesh,
+# where $x = -a$ and $x = a$, respectively. The `locate_entities_boundary` function returns an array of facet
+# indices representing these identified boundary entities.
 left_facet_marker = dolfinx.mesh.locate_entities_boundary(msh, fdim, left)
 right_facet_marker = dolfinx.mesh.locate_entities_boundary(msh, fdim, right)
+
+# The `get_ds_bound_from_marker` function generates a measure for applying boundary conditions 
+# specifically to the facets identified by `left_facet_marker` and `right_facet_marker`, respectively. 
+# This measure is then assigned to `ds_left` and `ds_right`.
 ds_left = get_ds_bound_from_marker(left_facet_marker, msh, fdim)
 ds_right = get_ds_bound_from_marker(right_facet_marker, msh, fdim)
+
+# `ds_list` is an array that stores boundary condition measures along with names 
+# for each boundary, simplifying result-saving processes. Each entry in `ds_list` 
+# is formatted as `[ds_, "name"]`, where `ds_` represents the boundary condition measure, 
+# and `"name"` is a label used for saving. Here, `ds_left` and `ds_right` are labeled 
+# as `"left"` and `"right"`, respectively, to ensure clarity when saving results.
 ds_list = np.array([
                    [ds_left, "left"],
                    [ds_right, "right"]
@@ -127,9 +158,19 @@ def right(x):
 
 
 ###############################################################################
-# Boundary Conditions
-# -------------------
-# Apply boundary conditions:.
+# Boundary Condition Setup for Scalar Field $\phi$
+# ------------------------------------------------
+# We define and apply a Dirichlet boundary condition for the scalar field $\phi$
+# on the left side of the mesh, setting $phi = -1$ on this boundary and on the right side of the mesh $phi = 1$ on this boundary. This setup is 
+# for a simple, static linear problem, meaning the boundary conditions and loading 
+# are constant and do not change throughout the simulation.
+#
+# - `bc_phi` is a function that creates a Dirichlet boundary condition on a specified 
+#   facet of the mesh for the scalar field $\phi$.
+# - `bcs_list_phi` is a list that stores all the boundary conditions for $\phi$, 
+#   facilitating easy management and extension of conditions if needed.
+# - `update_boundary_conditions` and `update_loading` are set to `None` as they are 
+#   unused in this static case with constant boundary conditions and loading.
 bc_left = bc_phi(left_facet_marker, V_phi, fdim, value=-1.0)
 bc_right = bc_phi(right_facet_marker, V_phi, fdim, value=1.0)
 bcs_list_phi = [bc_left, bc_right]
@@ -138,9 +179,31 @@ def right(x):
 
 
 ###############################################################################
-# Call the solver
-# ---------------
-# We set a final time of t=1 and a time step of \Delta t=1.
+# Solver Call for a Static Linear Problem
+# ---------------------------------------
+# We define the parameters for a simple, static linear boundary value problem 
+# with a final time `t = 1.0` and a time step `Δt = 1.0`. Although this setup 
+# includes time parameters, they are primarily used for structural consistency 
+# with a generic solver function and do not affect the result, as the problem 
+# is linear and time-independent.
+#
+# Parameters:
+# - `final_time`: The end time for the simulation, set to 1.0.
+# - `dt`: The time step for the simulation, set to 1.0. In a static context, this
+#   only provides uniformity with dynamic cases but does not change the results.
+# - `path`: Optional path for saving results; set to `None` here to use the default.
+# - `quadrature_degree`: Defines the accuracy of numerical integration; set to 2 
+#   for this problem.
+#
+# Function Call:
+# The `solve` function is called with:
+# - `Data`: Simulation data and parameters.
+# - `msh`: Mesh of the domain.
+# - `V_phi`: Function space for `phi`.
+# - `bcs_list_phi`: List of boundary conditions.
+# - `update_boundary_conditions`, `update_loading`: Set to `None` as they are unused in this static problem.
+# - `ds_list`: Boundary measures for integration on specified boundaries.
+# - `dt` and `final_time` to define the static solution time window.
 
 final_time = 1.0
 dt = 1.0