diff --git a/notebooks/geoclaw/topotools_examples.ipynb b/notebooks/geoclaw/topotools_examples.ipynb
index c2aa674..61af6dd 100644
--- a/notebooks/geoclaw/topotools_examples.ipynb
+++ b/notebooks/geoclaw/topotools_examples.ipynb
@@ -238,6 +238,7 @@
    "metadata": {},
    "source": [
     "<div id=crop></div>\n",
+    "\n",
     "## Cropping topography\n",
     "\n",
     "Sometimes it is useful to crop a large topo file to create a smaller one that contains a sub-region, e.g."
@@ -302,6 +303,7 @@
    "metadata": {},
    "source": [
     "<div id=coarsen></div>\n",
+    "\n",
     "## Coarsening topography\n",
     "\n",
     "The `crop` function also has an optional argument `coarsen` with default value 1.  If the value is  larger integer than it coarsens the grid by that factor in each direction.\n",
@@ -332,6 +334,7 @@
    "metadata": {},
    "source": [
     "<div id=unstructured></div>\n",
+    "\n",
     "## Unstructured Topography"
    ]
   },
@@ -350,18 +353,30 @@
    "source": [
     "import scipy.interpolate as interpolate\n",
     "\n",
+    "# Here we need to transpose the Z array so that the coordinate evaluation is correct\n",
+    "bathy = interpolate.RegularGridInterpolator((topo.x, topo.y), topo.Z.transpose())"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now choose `N` random points and evaluate the interpolant at these points:"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
     "N = 10000\n",
     "\n",
     "x = numpy.random.rand(N) * (topo.x[-1] - topo.x[0]) + topo.x[0]\n",
     "y = numpy.random.rand(N) * (topo.y[-1] - topo.y[0]) + topo.y[0]\n",
     "\n",
-    "# Here we need to transpose the Z array so that the coordinate evaluation is correct\n",
-    "bathy = interpolate.RectBivariateSpline(topo.x, topo.y, topo.Z.transpose())\n",
-    "\n",
-    "# strange hack to evaluate the function `bathy` produced by RectBivariateSpline since bathy(x, y) does not work\n",
-    "z = numpy.empty(x.shape)\n",
-    "for i in range(x.shape[0]):\n",
-    "    z[i] = bathy(x[i], y[i])"
+    "xypts = array(list(zip(x,y)))\n",
+    "z = bathy(xypts)"
    ]
   },
   {
@@ -475,9 +490,9 @@
  ],
  "metadata": {
   "kernelspec": {
-   "display_name": "Python 2",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
-   "name": "python2"
+   "name": "python3"
   },
   "language_info": {
    "codemirror_mode": {
@@ -489,9 +504,9 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.6.2"
+   "version": "3.10.13"
   }
  },
  "nbformat": 4,
- "nbformat_minor": 1
+ "nbformat_minor": 4
 }