diff --git a/figures.ipynb b/figures.ipynb
index 799e8c61bfd61240cc16285cb9431f8c2b71a4ee..49bdf4d2be74c375ee3c4ce09db70a3c60e77533 100644
--- a/figures.ipynb
+++ b/figures.ipynb
@@ -111,15 +111,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import adaptive"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
+    "import adaptive\n",
+    "\n",
     "def f(x, offset=0.123):\n",
     "    a = 0.02\n",
     "    return x + a**2 / (a**2 + (x - offset)**2)\n",
@@ -130,102 +123,33 @@
     "def h(x):\n",
     "    return np.sin(100*x) * np.exp(-x**2 / 0.1**2)\n",
     "\n",
-    "funcs = [dict(function=f, bounds=(-1, 1)), dict(function=g, bounds=(-1, 1)), dict(function=h, bounds=(-0.3, 0.3))]\n",
+    "funcs = [dict(function=f, bounds=(-1, 1), title=\"peak\"), dict(function=g, bounds=(-1, 1), title=\"tanh\"), dict(function=h, bounds=(-0.3, 0.3), title=\"wave packet\")]\n",
     "fig, axs = plt.subplots(2, len(funcs), figsize=(fig_width, 1.5*fig_height))\n",
     "n_points = 50\n",
     "for i, ax in enumerate(axs.T.flatten()):\n",
+    "    ax.xaxis.set_ticks([])\n",
+    "    ax.yaxis.set_ticks([])\n",
     "    if i % 2 == 0:\n",
     "        d = funcs[i // 2]\n",
     "        # homogeneous\n",
     "        xs = np.linspace(*d['bounds'], n_points)\n",
     "        ys = d['function'](xs)\n",
+    "        ax.set_title(rf\"\\textrm{{{d['title']}}}\")\n",
     "    else:\n",
     "        d = funcs[(i - 1) // 2]\n",
     "        loss = adaptive.learner.learner1D.curvature_loss_function()\n",
-    "        learner = adaptive.Learner1D(**d, loss_per_interval=loss)\n",
+    "        learner = adaptive.Learner1D(d['function'], bounds=d['bounds'], loss_per_interval=loss)\n",
     "        adaptive.runner.simple(learner, goal=lambda l: l.npoints >= n_points)\n",
     "        # adaptive\n",
     "        xs, ys = zip(*sorted(learner.data.items()))\n",
     "    xs_dense = np.linspace(*d['bounds'], 1000)\n",
     "    ax.plot(xs_dense, d['function'](xs_dense), c='red', alpha=0.3, lw=0.5)\n",
-    "#     ax.plot(xs, ys, c='k', alpha=0.3, lw=0.3)\n",
-    "    ax.scatter(xs, ys, s=0.5, c='k')"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "import adaptive\n",
-    "adaptive.notebook_extension()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "def f(x, offset=0.12312):\n",
-    "    a = 0.01\n",
-    "    return x + a**2 / (a**2 + (x - offset)**2)\n",
-    "\n",
-    "learner = adaptive.Learner1D(f, bounds=(-1, 1))\n",
-    "adaptive.runner.simple(learner, goal=lambda l: l.npoints > 100)\n"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "xs, ys = zip(*learner.data.items())"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "fig, ax = plt.subplots()\n",
-    "for i in range(1, len(xs)):\n",
-    "    if i % 10 != 0:\n",
-    "        continue\n",
-    "    alpha = np.linspace(0.2, 1, 101)[i]\n",
-    "    offset = i / len(xs)\n",
-    "    xs_part, ys_part = xs[:i], ys[:i]\n",
-    "    xs_part, ys_part = zip(*sorted(zip(xs_part, ys_part)))\n",
-    "    ax.plot(xs_part, offset + np.array(ys_part), alpha=alpha, c='grey', lw=0.5)\n",
-    "\n",
-    "plt.show()"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": [
-    "xs_part, ys_part"
+    "    ax.scatter(xs, ys, s=0.5, c='k')\n",
+    "    \n",
+    "axs[0][0].set_ylabel(r'$\\textrm{homogeneous}$')\n",
+    "axs[1][0].set_ylabel(r'$\\textrm{adaptive}$')\n",
+    "plt.savefig(\"figures/adaptive_vs_grid.pdf\", bbox_inches=\"tight\", transparent=True)\n"
    ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": null,
-   "metadata": {},
-   "outputs": [],
-   "source": []
   }
  ],
  "metadata": {
diff --git a/paper.md b/paper.md
index d066a98e0b69d73299d12de2581017e10cf593db..684493df2e4238f1f6e4b6b50c2f093e75862fd3 100755
--- a/paper.md
+++ b/paper.md
@@ -53,6 +53,10 @@ Each candidate point has a loss $L$ indicated by the size of the red dots.
 The candidate point with the largest loss will be chosen, which in this case is the one with $L_{1,2}$.
 ](figures/loss_1D.pdf){#fig:loss_1D}
 
+![Comparison of homogeneous sampling (top) with adaptive sampling (bottom) for different one-dimensional functions (red).
+We see that when the function has a distince feature---such as with the peak and tanh---adaptive sampling performs much better.
+When the features are homogeneously spaced, such as with the wave packet, adaptive sampling is not as effective as in the other cases.](figures/adaptive_vs_grid.pdf){#fig:adaptive_vs_grid}
+
 #### We provide a reference implementation, the Adaptive package, and demonstrate its performance.
 We provide a reference implementation, the open-source Python package called Adaptive[@Nijholt2019a], which has previously been used in several scientific publications[@vuik2018reproducing; @laeven2019enhanced; @bommer2019spin; @melo2019supercurrent].
 It has algorithms for $f \colon \R^N \to \R^M$, where $N, M \in \mathbb{Z}^+$ but which work best when $N$ is small; integration in $\R$; and the averaging of stochastic functions.