diff --git a/chapter1/poisson.ipynb b/chapter1/poisson.ipynb
index c8ec097..281da08 100644
--- a/chapter1/poisson.ipynb
+++ b/chapter1/poisson.ipynb
@@ -55,12 +55,16 @@
     "from sympde.topology import ScalarFunctionSpace, Square, element_of\n",
     "from sympde.calculus import grad, dot\n",
     "\n",
-    "from sympy import pi, sin\n",
+    "import matplotlib.pyplot as plt\n",
+    "import numpy as np\n",
+    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
     "\n",
-    "domain = Square()\n",
+    "from sympy import pi, sin, lambdify\n",
     "\n",
-    "V = ScalarFunctionSpace('V', domain)\n",
+    "domain = Square()\n",
     "\n",
+    "V = ScalarFunctionSpace('V', domain) \n",
+    "                                     \n",
     "x,y = domain.coordinates\n",
     "\n",
     "u,v = [element_of(V, name=i) for i in ['u', 'v']]\n",
@@ -202,11 +206,11 @@
    "outputs": [],
    "source": [
     "ue = sin(pi*x)*sin(pi*y)\n",
-    "\n",
     "u = element_of(V, name='u')\n",
     "\n",
     "# create the formal Norm object\n",
     "l2norm = Norm(u - ue, domain, kind='l2')\n",
+    "print(l2norm)\n",
     "\n",
     "# discretize the norm\n",
     "l2norm_h = discretize(l2norm, domain_h, Vh)\n",
@@ -257,7 +261,7 @@
   {
    "cell_type": "code",
    "execution_count": null,
-   "id": "d829e410",
+   "id": "827c3e69-77ac-4312-a4dd-a1c26a40b27c",
    "metadata": {},
    "outputs": [],
    "source": [
@@ -273,9 +277,120 @@
     "# 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": "968a2165-ab75-4bf0-b943-daab3bd9fa2d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pyplot as plt\n",
+    "from mpl_toolkits.axes_grid1 import make_axes_locatable\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, ue, 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",
+    "        ue: The exact solution function.\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",
+    "ue = lambdify((x, y), u, 'numpy') # convert sympy function to a callable function \n",
+    "plot_solutions(Vh, uh, ue, N=10)"
+   ]
   }
  ],
- "metadata": {},
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.12"
+  }
+ },
  "nbformat": 4,
  "nbformat_minor": 5
 }