diff --git a/chapter1/poisson.ipynb b/chapter1/poisson.ipynb
index c8ec097..4e39175 100644
--- a/chapter1/poisson.ipynb
+++ b/chapter1/poisson.ipynb
@@ -257,7 +257,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "d829e410",
+   "id": "827c3e69-77ac-4312-a4dd-a1c26a40b27c",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -273,6 +273,104 @@
     "# print the result\n",
     "print(h1_error)"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "id": "14d92a71-8c0f-4a91-8f56-c405f3fc76d1",
+   "metadata": {},
+   "source": [
+    "### Visualization\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2a33eab4-f4b8-428d-ae35-2d086d5645c3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
+    "from sympy import lambdify\n",
+    "\n",
+    "def refine_array_1d(breaks, N):\n",
+    "    \"\"\"Refine a 1D array by inserting N points between each pair of original values.\"\"\"\n",
+    "    refined = np.concatenate([np.linspace(breaks[i], breaks[i+1], N+1)[:-1] for i in range(len(breaks)-1)])\n",
+    "    return np.append(refined, breaks[-1])  # Add the last point\n",
+    "\n",
+    "def plot_solutions(Vh, uh, phi_exact, N=10):\n",
+    "    \"\"\"\n",
+    "    Refine the grid, evaluate solutions, compute errors, and plot results.\n",
+    "\n",
+    "    Parameters:\n",
+    "        Vh: The finite element space (must have `spaces` attribute with `breaks`).\n",
+    "        uh: The numerical solution function.\n",
+    "        phi_exact: The exact solution function (callable).\n",
+    "        N: Number of points to insert between breaks (default: 10).\n",
+    "    \"\"\"\n",
+    "    # Generate refined grid for visualization\n",
+    "    eta1 = refine_array_1d(Vh.spaces[0].breaks, N)\n",
+    "    eta2 = refine_array_1d(Vh.spaces[1].breaks, N)\n",
+    "\n",
+    "    # Evaluate numerical solution on the refined grid\n",
+    "    num = np.array([[uh(e1, e2) for e2 in eta2] for e1 in eta1])\n",
+    "\n",
+    "    # Compute exact solution\n",
+    "    ex = np.array([[phi_exact(e1, e2) for e2 in eta2] for e1 in eta1])\n",
+    "    err = num - ex\n",
+    "\n",
+    "    # Generate physical grid coordinates\n",
+    "    xx, yy = np.meshgrid(eta1, eta2, indexing='ij')\n",
+    "\n",
+    "    # Create figure with 3 subplots\n",
+    "    fig, axes = plt.subplots(1, 3, figsize=(12.8, 4.8))\n",
+    "\n",
+    "    def add_colorbar(im, ax):\n",
+    "        \"\"\"Add a colorbar to the given axis.\"\"\"\n",
+    "        divider = make_axes_locatable(ax)\n",
+    "        cax = divider.append_axes(\"right\", size=0.2, pad=0.2)\n",
+    "        cbar = ax.get_figure().colorbar(im, cax=cax)\n",
+    "        return cbar\n",
+    "\n",
+    "    # Plot exact solution\n",
+    "    ax = axes[0]\n",
+    "    im = ax.contourf(xx, yy, ex, 40, cmap='jet')\n",
+    "    add_colorbar(im, ax)\n",
+    "    ax.set_xlabel(r'$x$', rotation='horizontal')\n",
+    "    ax.set_ylabel(r'$y$', rotation='horizontal')\n",
+    "    ax.set_title(r'$\\phi_{exact}(x,y)$')\n",
+    "    ax.set_aspect('equal')\n",
+    "\n",
+    "    # Plot numerical solution\n",
+    "    ax = axes[1]\n",
+    "    im = ax.contourf(xx, yy, num, 40, cmap='jet')\n",
+    "    add_colorbar(im, ax)\n",
+    "    ax.set_xlabel(r'$x$', rotation='horizontal')\n",
+    "    ax.set_ylabel(r'$y$', rotation='horizontal')\n",
+    "    ax.set_title(r'$\\phi(x,y)$')\n",
+    "    ax.set_aspect('equal')\n",
+    "\n",
+    "    # Plot numerical error\n",
+    "    ax = axes[2]\n",
+    "    im = ax.contourf(xx, yy, err, 40, cmap='jet')\n",
+    "    add_colorbar(im, ax)\n",
+    "    ax.set_xlabel(r'$x$', rotation='horizontal')\n",
+    "    ax.set_ylabel(r'$y$', rotation='horizontal')\n",
+    "    ax.set_title(r'$\\phi(x,y) - \\phi_{exact}(x,y)$')\n",
+    "    ax.set_aspect('equal')\n",
+    "\n",
+    "    # Show figure\n",
+    "    plt.tight_layout()\n",
+    "    fig.show()\n",
+    "\n",
+    "# Define the exact solution as a symbolic expression (example: phi_exact = x^2 + y^2)\n",
+    "ue = sin(np.pi*x)*sin(np.pi*y)\n",
+    "\n",
+    "# Convert the symbolic expression to a callable function using lambdify\n",
+    "phi_exact = lambdify((x, y), ue, 'numpy')\n",
+    "plot_solutions(Vh, uh, phi_exact, N=10)"
+   ]
   }
  ],
  "metadata": {},