diff --git a/examples/AllenCahn/plot_4000.py b/examples/AllenCahn/plot_4000.py
index 9d58eb4..30bf796 100644
--- a/examples/AllenCahn/plot_4000.py
+++ b/examples/AllenCahn/plot_4000.py
@@ -126,18 +126,21 @@ def right(x):
 
 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, 
@@ -166,11 +169,11 @@ def right(x):
 # 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$.
+# 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.
+# 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.
+# 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]
@@ -190,7 +193,7 @@ def right(x):
 # 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.
+# 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.
diff --git a/examples/Elasticity/plot_1100.py b/examples/Elasticity/plot_1100.py
index 23b71da..ae2be9a 100644
--- a/examples/Elasticity/plot_1100.py
+++ b/examples/Elasticity/plot_1100.py
@@ -75,9 +75,9 @@
 # - `E`: Young's modulus, set to 210 kN/mm².
 # - `nu`: Poisson's ratio, set to 0.3.
 # - `save_solution_xdmf` and `save_solution_vtu`: Set to `False` and `True`, respectively,
-#   specifying the file format to save displacement results (.vtu here).
+# specifying the file format to save displacement results (.vtu here).
 # - `results_folder_name`: Name of the folder for saving results. If it exists,
-#   it will be replaced with a new empty folder.
+# it will be replaced with a new empty folder.
 Data = Input(E=210.0,
              nu=0.3,
              save_solution_xdmf=False,
@@ -115,18 +115,21 @@ def top(x):
 
 fdim = msh.topology.dim - 1
 
+# %%
 # Using the `bottom` and `top` functions, we locate the facets on the bottom and top sides of the mesh,
 # where $y = 0$ and $y = ly$, respectively. The `locate_entities_boundary` function returns an array of facet
 # indices representing these identified boundary entities.
 bottom_facet_marker = dolfinx.mesh.locate_entities_boundary(msh, fdim, bottom)
 top_facet_marker = dolfinx.mesh.locate_entities_boundary(msh, fdim, top)
 
+# %%
 # The `get_ds_bound_from_marker` function generates a measure for applying boundary conditions 
 # specifically to the facets identified by `top_facet_marker` and `bottom_facet_marker`, respectively. 
 # This measure is then assigned to `ds_bottom` and `ds_top`.
 ds_bottom = get_ds_bound_from_marker(bottom_facet_marker, msh, fdim)
 ds_top = get_ds_bound_from_marker(top_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, 
@@ -155,6 +158,7 @@ def top(x):
 bc_bottom = bc_xy(bottom_facet_marker, V_u, fdim, value_x=0.0, value_y=0.0)
 bc_top = bc_xy(top_facet_marker, V_u, fdim, value_x=0.0, value_y=0.0)
 
+# %%
 # The bcs_list_u variable is a list that stores all boundary conditions for the displacement
 # field $\boldsymbol u$. This list facilitates easy management of multiple boundary
 # conditions and can be expanded if additional conditions are needed.
@@ -169,21 +173,21 @@ def top(x):
 #
 # Parameters:
 # - `bcs`: List of boundary conditions, where each entry corresponds to a boundary condition applied
-#   to a specific facet of the mesh.
+# to a specific facet of the mesh.
 # - `time`: Scalar representing the current time step in the analysis.
 #
 # Inside the function:
 # - `val` is calculated as a linear function of `time`, specifically `val = 0.0003 * time`,
-#   to simulate gradual displacement along the y-axis. This can be modified as needed for different
-#   quasi-static loading schemes.
+# to simulate gradual displacement along the y-axis. This can be modified as needed for different
+# quasi-static loading schemes.
 # - The value `val` is assigned to the y-component of the displacement field on the boundary,
-#   achieved by updating `bcs[0].g.value[1]`, where `bcs[0]` represents the top boundary condition.
+# achieved by updating `bcs[0].g.value[1]`, where `bcs[0]` represents the top boundary condition.
 #
 # Returns:
 # - A tuple `(0, val, 0)` where:
-#   - The first element is zero (indicating no update for the x component in this example).
-#   - The second element is `val`, the calculated y-displacement.
-#   - The third element is zero (indicating no z-component displacement in this 2D example).
+# - The first element is zero (indicating no update for the x component in this example).
+# - The second element is `val`, the calculated y-displacement.
+# - The third element is zero (indicating no z-component displacement in this 2D example).
 #
 # This function supports the quasi-static analysis by gradually updating the displacement boundary
 # condition over time, allowing for controlled loading in the simulation.
@@ -209,7 +213,7 @@ def update_boundary_conditions(bcs, time):
 # Parameters:
 # - `final_time`: The end time for the simulation, set to 10.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.
+# 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.
diff --git a/examples/Elasticity/plot_1101.py b/examples/Elasticity/plot_1101.py
index 5ee24d6..eb88bb3 100644
--- a/examples/Elasticity/plot_1101.py
+++ b/examples/Elasticity/plot_1101.py
@@ -77,9 +77,9 @@
 # - `E`: Young's modulus, set to 210 kN/mm².
 # - `nu`: Poisson's ratio, set to 0.3.
 # - `save_solution_xdmf` and `save_solution_vtu`: Set to `False` and `True`, respectively,
-#   specifying the file format to save displacement results (.vtu here).
+# specifying the file format to save displacement results (.vtu here).
 # - `results_folder_name`: Name of the folder for saving results. If it exists,
-#   it will be replaced with a new empty folder.
+# it will be replaced with a new empty folder.
 Data = Input(E=210.0,
              nu=0.3,
              save_solution_xdmf=False,
@@ -117,18 +117,21 @@ def top(x):
 
 fdim = msh.topology.dim - 1
 
+# %%
 # Using the `bottom` and `top` functions, we locate the facets on the bottom and top sides of the mesh,
 # where $y = 0$ and $y = ly$, respectively. The `locate_entities_boundary` function returns an array of facet
 # indices representing these identified boundary entities.
 bottom_facet_marker = dolfinx.mesh.locate_entities_boundary(msh, fdim, bottom)
 top_facet_marker = dolfinx.mesh.locate_entities_boundary(msh, fdim, top)
 
+# %%
 # The `get_ds_bound_from_marker` function generates a measure for applying boundary conditions 
 # specifically to the facets identified by `top_facet_marker` and `bottom_facet_marker`, respectively. 
 # This measure is then assigned to `ds_bottom` and `ds_top`.
 ds_bottom = get_ds_bound_from_marker(top_facet_marker, msh, fdim)
 ds_top = get_ds_bound_from_marker(top_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, 
@@ -155,6 +158,7 @@ def top(x):
 # - `bc_bottom`: Fixes x and y displacement to 0 on the bottom boundary.
 bc_bottom = bc_xy(bottom_facet_marker, V_u, fdim)
 
+# %%
 # The bcs_list_u variable is a list that stores all boundary conditions for the displacement
 # field $\boldsymbol u$. This list facilitates easy management of multiple boundary
 # conditions and can be expanded if additional conditions are needed.
@@ -171,6 +175,7 @@ def update_boundary_conditions(bcs, time):
 # `T_top` represents the external force applied in the y-direction.
 T_top = loading_Txy(V_u, msh, ds_top)
 
+# %%
 # The load is added to the list of external loads, `T_list_u`, which will be updated
 # incrementally in the `update_loading` function.
 T_list_u = [[T_top, ds_top]
@@ -187,22 +192,22 @@ def update_boundary_conditions(bcs, time):
 #
 # Parameters:
 # - `T_list_u`: List of tuples where each entry corresponds to a load applied to a specific 
-#   boundary or facet of the mesh.
+# boundary or facet of the mesh.
 # - `time`: Scalar representing the current time step in the analysis.
 #
 # Inside the function:
 # - `val` is calculated as `0.1 * time`, a linear function of `time`, which represents the 
-#   gradual application of force in the y-direction. This scaling factor (`0.1` in this case) can
-#   be adjusted to control the rate of force increase.
+# gradual application of force in the y-direction. This scaling factor (`0.1` in this case) can
+# be adjusted to control the rate of force increase.
 # - The value `val` is assigned to the y-component of the external force field on the top boundary
-#   by setting `T_list_u[0][0].value[1]`, where `T_list_u[0][0]` represents the load applied to
-#   the designated top boundary facet (`ds_top`).
+# by setting `T_list_u[0][0].value[1]`, where `T_list_u[0][0]` represents the load applied to
+# the designated top boundary facet (`ds_top`).
 #
 # Returns:
 # - A tuple `(0, val, 0)` where:
-#   - The first element is zero, indicating no load in the x-direction.
-#   - The second element is `val`, the calculated y-directional force.
-#   - The third element is zero, as this is a 2D example without z-component loading.
+# - The first element is zero, indicating no load in the x-direction.
+# - The second element is `val`, the calculated y-directional force.
+# - The third element is zero, as this is a 2D example without z-component loading.
 #
 # This function supports force-controlled quasi-static analysis by adjusting the applied load 
 # over time, ensuring a controlled force increase in the simulation.
@@ -225,10 +230,10 @@ def update_loading(T_list_u, time):
 # Parameters:
 # - `final_time`: The end time for the simulation, set to 10.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.
+# 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.
+# for this problem.
 #
 # Function Call:
 # The `solve` function is called with:
diff --git a/examples/Fatigue/plot_1800.py b/examples/Fatigue/plot_1800.py
index f683412..eee80c9 100644
--- a/examples/Fatigue/plot_1800.py
+++ b/examples/Fatigue/plot_1800.py
@@ -189,9 +189,9 @@
 # -------------------
 # Dirichlet boundary conditions are defined as follows:
 # - `bc_bottom`: Constrains both x and y displacements to 0 on the bottom boundary,
-#   ensuring that the bottom edge is fixed.
+# ensuring that the bottom edge is fixed.
 # - `bc_top`: The vertical displacement on the top boundary is updated dynamically
-#   to impose the cyclic load.
+# to impose the cyclic load.
 
 # %%
 # The `bcs_list_u` variable is a list containing all boundary conditions for the 
@@ -206,15 +206,11 @@
 # Definition of the Cyclic Load
 # -----------------------------
 # The cyclic load is applied by updating the boundary condition at the top of the specimen. 
-# The cyclic load follows the form:
-# 
-# $$ \frac{2}{\pi} A \arcsin[\sin(\omega t)] $$
-# 
+# The cyclic load follows the form $ \frac{2}{\pi} A \arcsin[\sin(\omega t)] $.
 # where:
-# - \( A = 0.002 \, \text{mm} \): Amplitude of the load.
-# - \( f = \frac{1}{8} \): Frequency of the load in Hz.
-# - \( \omega = 2 \pi f \): Angular frequency.
-#
+# - $A = 0.002 \, \text{mm}$: Amplitude of the load.
+# - $f = \frac{1}{8}$: Frequency of the load in Hz.
+# - $\omega = 2 \pi f $: Angular frequency.
 # This periodic loading is imposed on the top boundary to simulate the desired cyclic behavior.
 
 amplitude = 0.002
@@ -232,23 +228,23 @@
 #
 # Parameters:
 # - `bcs`: A list of boundary conditions, where each entry corresponds to a specific facet of the mesh.
-#          For this implementation, `bcs[0]` refers to the top boundary condition.
+# For this implementation, `bcs[0]` refers to the top boundary condition.
 # - `time`: A scalar representing the current time step in the simulation.
 #
 # Inside the function:
 # - `val` is computed based on the cyclic load formula:
-#   $$ \text{val} = \frac{2}{\pi} \cdot \text{amplitude} \cdot \arcsin(\sin(\omega \cdot \text{time})) $$
-#   where:
-#   - `amplitude` defines the maximum displacement of the cyclic load.
-#   - `w` (omega) is the angular frequency, given by \( \omega = 2 \pi f \), where \( f \) is the frequency.
+# $$ \text{val} = \frac{2}{\pi} \cdot \text{amplitude} \cdot \arcsin(\sin(\omega \cdot \text{time})) $$
+# where:
+# - `amplitude` defines the maximum displacement of the cyclic load.
+# - `w` (omega) is the angular frequency, given by \( \omega = 2 \pi f \), where \( f \) is the frequency.
 # - The computed `val` represents the y-displacement applied to the boundary at the current time step.
 # - This value is dynamically updated in `bcs[0].g.value[...]` to apply the displacement to the top boundary.
 #
 # Returns:
 # - A tuple `(0, val, 0)`, where:
-#   - The first element is `0` (indicating no x-displacement).
-#   - The second element is `val`, the calculated y-displacement.
-#   - The third element is `0` (indicating no z-displacement, as this is a 2D simulation).
+# - The first element is `0` (indicating no x-displacement).
+# - The second element is `val`, the calculated y-displacement.
+# - The third element is `0` (indicating no z-displacement, as this is a 2D simulation).
 #
 # The function enables controlled cyclic loading during the simulation, facilitating the study
 # of fatigue and quasi-static response under varying boundary conditions.
@@ -283,11 +279,11 @@ def update_boundary_conditions(bcs, time):
 # **Parameters:**
 # - `dt`: The time step for the simulation, set to 1.0.
 # - `final_time`: The total simulation time, computed as \( 8 \cdot 200 + 1 \), ensuring
-#   sufficient steps for the cyclic loading behavior.
+# sufficient steps for the cyclic loading behavior.
 # - `path`: Optional parameter for specifying the results folder; set to `None` here 
-#   to use the default location.
+# to use the default location.
 # - `quadrature_degree`: (Defined elsewhere) Specifies the accuracy of numerical integration 
-#   during the computation. For this problem, it is set to 2.
+# during the computation. For this problem, it is set to 2.
 #
 # **Function Call:**
 # The `solve` function is invoked with the following arguments:
diff --git a/examples/GmshGeoFiles/plot_9101.py b/examples/GmshGeoFiles/plot_9101.py
index 3ea5e59..c8f721d 100644
--- a/examples/GmshGeoFiles/plot_9101.py
+++ b/examples/GmshGeoFiles/plot_9101.py
@@ -41,16 +41,21 @@
 import pyvista as pv
 
 folder = "9101_SingleNotchedTensionTest"
+
+# %%
 # Initialize Gmsh
 gmsh.initialize()
 
+# %%
 # Open the .geo file
 geo_file = os.path.join(folder, "file.geo")
 gmsh.open(geo_file)
 
+# %%
 # Generate the mesh (2D example, for 3D use generate(3))
 gmsh.model.mesh.generate(3)
 
+# %%
 # Write the mesh to a .vtk file for visualization
 # Note that the input mesh file for the *phasefieldx* simulation should have the .msh extension.
 # Use "output_mesh_for_view.msh" to generate the mesh for the simulation input.
@@ -58,6 +63,7 @@
 vtu_file = os.path.join(folder, "output_mesh_for_view.vtk")
 gmsh.write(vtu_file)
 
+# %%
 # Finalize Gmsh
 gmsh.finalize()
 
diff --git a/examples/GmshGeoFiles/plot_9102.py b/examples/GmshGeoFiles/plot_9102.py
index 162ece3..4f91220 100644
--- a/examples/GmshGeoFiles/plot_9102.py
+++ b/examples/GmshGeoFiles/plot_9102.py
@@ -39,16 +39,21 @@
 import pyvista as pv
 
 folder = "9102_SingleNotchedShearTest"
+
+# %%
 # Initialize Gmsh
 gmsh.initialize()
 
+# %%
 # Open the .geo file
 geo_file = os.path.join(folder, "file.geo")
 gmsh.open(geo_file)
 
+# %%
 # Generate the mesh (2D example, for 3D use generate(3))
 gmsh.model.mesh.generate(3)
 
+# %%
 # Write the mesh to a .vtk file for visualization
 # Note that the input mesh file for the *phasefieldx* simulation should have the .msh extension.
 # Use "output_mesh_for_view.msh" to generate the mesh for the simulation input.
@@ -56,6 +61,7 @@
 vtu_file = os.path.join(folder, "output_mesh_for_view.vtk")
 gmsh.write(vtu_file)
 
+# %%
 # Finalize Gmsh
 gmsh.finalize()
 
diff --git a/examples/GmshGeoFiles/plot_9103.py b/examples/GmshGeoFiles/plot_9103.py
index 17e58fb..feaad44 100644
--- a/examples/GmshGeoFiles/plot_9103.py
+++ b/examples/GmshGeoFiles/plot_9103.py
@@ -42,20 +42,26 @@
 import pyvista as pv
 
 folder = "9103_ThreePointBendingTest"
+
+# %%
 # Initialize Gmsh
 gmsh.initialize()
 
+# %%
 # Open the .geo file
 geo_file = os.path.join(folder, "file.geo")
 gmsh.open(geo_file)
 
+# %%
 # Generate the mesh (2D example, for 3D use generate(3))
 gmsh.model.mesh.generate(3)
 
+# %%
 # Write the mesh to a .vtu file
 vtu_file = os.path.join(folder, "output_mesh_for_view.vtk")
 gmsh.write(vtu_file)
 
+# %%
 # Finalize Gmsh
 gmsh.finalize()
 
diff --git a/examples/GmshGeoFiles/plot_9104.py b/examples/GmshGeoFiles/plot_9104.py
index 4c6c502..48d8f53 100644
--- a/examples/GmshGeoFiles/plot_9104.py
+++ b/examples/GmshGeoFiles/plot_9104.py
@@ -42,16 +42,21 @@
 import pyvista as pv
 
 folder = "9104_Structured"
+
+# %%
 # Initialize Gmsh
 gmsh.initialize()
 
+# %%
 # Open the .geo file
 geo_file = os.path.join(folder, "file.geo")
 gmsh.open(geo_file)
 
+# %%
 # Generate the mesh (2D example, for 3D use generate(3))
 gmsh.model.mesh.generate(2)
 
+# %%
 # Write the mesh to a .vtk file for visualization
 # Note that the input mesh file for the *phasefieldx* simulation should have the .msh extension.
 # Use "output_mesh_for_view.msh" to generate the mesh for the simulation input.
@@ -59,6 +64,7 @@
 vtu_file = os.path.join(folder, "output_mesh_for_view.vtk")
 gmsh.write(vtu_file)
 
+# %%
 # Finalize Gmsh
 gmsh.finalize()
 
diff --git a/examples/GmshGeoFiles/plot_9106.py b/examples/GmshGeoFiles/plot_9106.py
new file mode 100644
index 0000000..32fb5d2
--- /dev/null
+++ b/examples/GmshGeoFiles/plot_9106.py
@@ -0,0 +1,54 @@
+"""
+.. _ref_9105:
+
+.geo File: Rectangle
+^^^^^^^^^^^^^^^^^^^^
+
+.geo file
+---------
+
+.. include::  ../../../../examples/GmshGeoFiles/9104_Structured/file.geo
+   :literal:
+
+"""
+
+###############################################################################
+# Mesh Visualization
+# ------------------
+# The purpose of this code is to visualize the mesh. The mesh is generated from
+# the .geo file and saved as output_mesh_for_view.vtu. It is then loaded and
+# visualized using PyVista.
+
+import os
+import gmsh
+import pyvista as pv
+
+folder = "9105_a"
+
+# %%
+# Initialize Gmsh
+gmsh.initialize()
+
+# %%
+# Open the .geo file
+geo_file = os.path.join(folder, "file.geo")
+gmsh.open(geo_file)
+
+# %%
+# Generate the mesh (2D example, for 3D use generate(3))
+gmsh.model.mesh.generate()
+
+# %%
+# Write the mesh to a .vtu file
+vtu_file = os.path.join(folder, "output_mesh_for_view.vtk")
+gmsh.write(vtu_file)
+
+# %%
+# Finalize Gmsh
+gmsh.finalize()
+
+print(f"Mesh successfully written to {vtu_file}")
+
+#pv.start_xvfb()
+file_vtu = pv.read(vtu_file)
+file_vtu.plot(cpos='xy', color='white', show_edges=True)
diff --git a/examples/PhaseField/plot_2000.py b/examples/PhaseField/plot_2000.py
index 19d2d6b..0cb5ec3 100644
--- a/examples/PhaseField/plot_2000.py
+++ b/examples/PhaseField/plot_2000.py
@@ -111,11 +111,11 @@
 # corresponding line, quadrilateral, or hexahedral elements.
 #
 # - `divx`, `divy`, and `divz` specify the number of divisions along the x, y, and z 
-#   axes, respectively. Here, `divx=100`, `divy=1`, and `divz=1` are set to divide 
-#   the x-axis primarily, as needed for a 2D or 1D mesh.
+# axes, respectively. Here, `divx=100`, `divy=1`, and `divz=1` are set to divide 
+# the x-axis primarily, as needed for a 2D or 1D mesh.
 #
 # - `lx`, `ly`, and `lz` define the physical dimensions of the domain in the x, y, and z 
-#   directions. In this example, we set `lx=5.0`, `ly=1.0`, and `lz=1.0`.
+# directions. In this example, we set `lx=5.0`, `ly=1.0`, and `lz=1.0`.
 #
 # Specify the simulation dimension with the `dimension` variable (`"1d"`, `"2d"`, or `"3d"`).
 # Here, we choose `"2d"`.
@@ -165,26 +165,26 @@
 def left(x):
     return np.equal(x[0], 0)
 
-#
+# %%
 # `fdim` represents the dimension of the boundary facets on the mesh, which is one 
 # less than the mesh's overall dimensionality (`msh.topology.dim`). For example, 
 # if the mesh is 2D, `fdim` will be 1, representing 1D boundary edges.
-#
 fdim = msh.topology.dim - 1
 
+# %%
 # Using the `left` function, we locate the facets on the left side of the mesh 
 # where `x=0`. The `locate_entities_boundary` function returns an array of facet 
 # indices that represent the identified boundary entities.
 left_facet_marker = dolfinx.mesh.locate_entities_boundary(msh, fdim, left)
 
-
+# %%
 # `get_ds_bound_from_marker` is a function that generates a measure for integrating 
 # boundary conditions specifically on the facets identified by `left_facet_marker`. 
 # This measure is assigned to `ds_left` and will be used for applying boundary 
 # conditions on the left side.
 ds_left = get_ds_bound_from_marker(left_facet_marker, msh, fdim)
 
-
+# %%
 # `ds_list` is an array that stores boundary condition measures and associated 
 # names for each boundary to facilitate result-saving processes. Each entry in 
 # `ds_list` is an array in the form `[ds_, "name"]`, where `ds_` is the boundary 
@@ -212,11 +212,11 @@ def left(x):
 # 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$.
+# 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.
+# 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.
+# unused in this static case with constant boundary conditions and loading.
 
 bc_left = bc_phi(left_facet_marker, V_phi, fdim, value=1.0)
 bcs_list_phi = [bc_left]
@@ -236,10 +236,10 @@ def left(x):
 # 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.
+# 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.
+# for this problem.
 #
 # Function Call:
 # The `solve` function is called with:
@@ -320,11 +320,11 @@ def left(x):
 # The energy components are calculated for the scalar field `phi` and its gradient. 
 # We compare the following energy terms:
 # - `W_phi`: The energy associated with the scalar field `phi`.
-#            $W_{\phi} = \frac{1}{2l} \int_{-a}^{a}   \left[ e^{-|x|/l} + \frac{1}{e^{\frac{2a}{l}}+1} 2 \sinh \left( \frac{|x|}{l} \right) \right]^2 dx$
+# $W_{\phi} = \frac{1}{2l} \int_{-a}^{a}   \left[ e^{-|x|/l} + \frac{1}{e^{\frac{2a}{l}}+1} 2 \sinh \left( \frac{|x|}{l} \right) \right]^2 dx$
 # - `W_gradphi`: The energy associated with the gradient of the scalar field.
-#                $W_{\nabla \phi}(\phi) = \frac{l}{2} \int_{-a}^{a}   \left[ \frac{-\text{sign}(x)}{l} e^{-|x|/l} + \frac{1}{e^{\frac{2a}{l}} +1} \frac{\text{sign}(x)}{l} 2 \cosh\left(\frac{|x|}{l}\right) \right]^2 dx$
+# $W_{\nabla \phi}(\phi) = \frac{l}{2} \int_{-a}^{a}   \left[ \frac{-\text{sign}(x)}{l} e^{-|x|/l} + \frac{1}{e^{\frac{2a}{l}} +1} \frac{\text{sign}(x)}{l} 2 \cosh\left(\frac{|x|}{l}\right) \right]^2 dx$
 # - `W`: The total energy, which includes contributions from both the field and its gradient.
-#        $W = \tanh \left( \frac{a}{l} \right)$
+# $W = \tanh \left( \frac{a}{l} \right)$
 #
 # Theoretical expressions are computed for these energies as functions of the 
 # length scale parameter `l`. Additionally, energy values from the simulation 
@@ -334,14 +334,17 @@ def left(x):
 a_div_l = lx / l_array  # Compute the ratio of `lx` to the length scale
 tanh_a_div_l = np.tanh(a_div_l)  # Compute the hyperbolic tangent of the ratio
 
+# %%
 # Compute the theoretical energy for the scalar field `phi`
 # The expression represents the energy of `phi` based on the length scale parameter `l`.
 W_phi = 0.5 * tanh_a_div_l + 0.5 * a_div_l * (1.0 - tanh_a_div_l**2)
 
+# %%
 # Compute the theoretical energy for the gradient of `phi`
 # This energy is based on the gradient of `phi` and the length scale parameter `l`.
 W_gradphi = 0.5 * tanh_a_div_l - 0.5 * a_div_l * (1.0 - tanh_a_div_l**2)
 
+# %%
 # Compute the total theoretical energy
 # This energy is related to the hyperbolic tangent of the ratio `a/l`.
 W = np.tanh(a_div_l)
@@ -357,20 +360,24 @@ def left(x):
 
 fig, energy = plt.subplots()  # Create a figure for plotting energy
 
+# %%
 # Enable LaTeX rendering for plot labels
 plt.rc('text', usetex=True)
 
+# %%
 # Plot the theoretical energy values
 energy.plot(l_array, W_phi, 'r-', label=r'$W_{\phi}$')  # Energy for `phi`
 energy.plot(l_array, W_gradphi, 'b-', label=r'$W_{\nabla \phi}$')  # Energy for gradient of `phi`
 energy.plot(l_array, W, 'k-', label=r'$W$')  # Total energy
 
+# %%
 # Plot the energy values obtained from the simulation
 # The values are scaled by 2 and correspond to the total energy components.
 energy.plot(Data.l, 2 * S.energy_files['total.energy']["gamma_phi"][0], 'k*', label=S.label)
 energy.plot(Data.l, 2 * S.energy_files['total.energy']["gamma_gradphi"][0], 'k*')
 energy.plot(Data.l, 2 * S.energy_files['total.energy']["gamma"][0], 'k*')
 
+# %%
 # Set plot labels and grid
 energy.set_xlabel('Length Scale Parameter $l$')  # X-axis label
 energy.set_ylabel('Energy')  # Y-axis label
diff --git a/examples/PhaseField/plot_2001.py b/examples/PhaseField/plot_2001.py
index b7a1688..c3ce47d 100644
--- a/examples/PhaseField/plot_2001.py
+++ b/examples/PhaseField/plot_2001.py
@@ -83,6 +83,7 @@
 divx, divy = 100, 50
 lx, ly = 1.0, 0.5
 
+# %%
 # Create a 2D mesh using the defined parameters.
 msh = dolfinx.mesh.create_rectangle(mpi4py.MPI.COMM_WORLD,
                                     [np.array([0, 0]),
@@ -100,22 +101,26 @@
 def bottom(x):
     return np.logical_and(np.isclose(x[1], 0), np.less(x[0], 0.5))
 
+# %%
 # `fdim` represents the dimension of the boundary facets on the mesh, which is one 
 # less than the mesh's overall dimensionality (`msh.topology.dim`). For example, 
 # if the mesh is 2D, `fdim` will be 1, representing 1D boundary edges.
 fdim = msh.topology.dim - 1
 
+# %%
 # Using the `bottom` function, we locate the facets on the bottom boundary side of the mesh.
 # The `locate_entities_boundary` function returns an array of facet 
 # indices that represent the identified boundary entities.
 bottom_facet_marker = dolfinx.mesh.locate_entities_boundary(msh, fdim, bottom)
 
+# %%
 # `get_ds_bound_from_marker` is a function that generates a measure for integrating 
 # boundary conditions specifically on the facets identified by `bottom_facet_marker`. 
 # This measure is assigned to `ds_bottom` and will be used for applying boundary 
 # conditions on the left side.
 ds_bottom = get_ds_bound_from_marker(bottom_facet_marker, msh, fdim)
 
+# %%
 # `ds_list` is an array that stores boundary condition measures and associated 
 # names for each boundary to facilitate result-saving processes. Each entry in 
 # `ds_list` is an array in the form `[ds_, "name"]`, where `ds_` is the boundary 
@@ -143,11 +148,11 @@ def bottom(x):
 # 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$.
+# 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.
+# 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.
+# unused in this static case with constant boundary conditions and loading.
 
 bc_bottom = bc_phi(bottom_facet_marker, V_phi, fdim, value=1.0)
 bcs_list_phi = [bc_bottom]
@@ -167,10 +172,10 @@ def bottom(x):
 # 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.
+# 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.
+# for this problem.
 #
 # Function Call:
 # The `solve` function is called with:
@@ -179,7 +184,7 @@ def bottom(x):
 # - `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.
+# 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.