diff --git a/doc/source/tutorial/plotting.rst b/doc/source/tutorial/plotting.rst
index f8bb6586a1c2b7d6cd30cf084d983b5eef8b3ce5..ef1ec556e9c779a07895678ed617d7b851451774 100644
--- a/doc/source/tutorial/plotting.rst
+++ b/doc/source/tutorial/plotting.rst
@@ -12,15 +12,52 @@ these options can be used to achieve various very different objectives.
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`plot_graphene.py </code/download/plot_graphene.py>`
+    :jupyter-download:script:`plot_graphene`
+
+.. jupyter-kernel::
+    :id: plot_graphene
+
+.. jupyter-execute::
+    :hide-code:
+
+    # Tutorial 2.8.1. 2D example: graphene quantum dot
+    # ================================================
+    #
+    # Physics background
+    # ------------------
+    #  - graphene edge states
+    #
+    # Kwant features highlighted
+    # --------------------------
+    #  - demonstrate different ways of plotting
+
+    from matplotlib import pyplot
+
+    import kwant
 
 We begin by first considering a circular graphene quantum dot (similar to what
 has been used in parts of the tutorial :ref:`tutorial-graphene`.)  In contrast
 to previous examples, we will also use hoppings beyond next-nearest neighbors:
 
-.. literalinclude:: /code/include/plot_graphene.py
-    :start-after: #HIDDEN_BEGIN_makesyst
-    :end-before: #HIDDEN_END_makesyst
+.. jupyter-execute::
+
+    lat = kwant.lattice.honeycomb()
+    a, b = lat.sublattices
+
+    def make_system(r=8, t=-1, tp=-0.1):
+
+        def circle(pos):
+            x, y = pos
+            return x**2 + y**2 < r**2
+
+        syst = kwant.Builder()
+        syst[lat.shape(circle, (0, 0))] = 0
+        syst[lat.neighbors()] = t
+        syst.eradicate_dangling()
+        if tp:
+            syst[lat.neighbors(2)] = tp
+
+        return syst
 
 Note that adding hoppings hoppings to the `n`-th nearest neighbors can be
 simply done by passing `n` as an argument to
@@ -29,16 +66,14 @@ simply done by passing `n` as an argument to
 out of the shape. It is necessary to do so *before* adding the
 next-nearest-neighbor hopping [#]_.
 
-Of course, the system can be plotted simply with default settings:
-
-.. literalinclude:: /code/include/plot_graphene.py
-    :start-after: #HIDDEN_BEGIN_plotsyst1
-    :end-before: #HIDDEN_END_plotsyst1
-
-However, due to the richer structure of the lattice, this results in a rather
+Of course, the system can be plotted simply with default settings,
+however, due to the richer structure of the lattice, this results in a rather
 busy plot:
 
-.. image:: /code/figure/plot_graphene_syst1.*
+.. jupyter-execute::
+
+    syst = make_system()
+    kwant.plot(syst);
 
 A much clearer plot can be obtained by using different colors for both
 sublattices, and by having different line widths for different hoppings.  This
@@ -47,15 +82,19 @@ can be achieved by passing a function to the arguments of
 must be a function taking one site as argument, for hoppings a function taking
 the start end end site of hopping as arguments:
 
-.. literalinclude:: /code/include/plot_graphene.py
-    :start-after: #HIDDEN_BEGIN_plotsyst2
-    :end-before: #HIDDEN_END_plotsyst2
+.. jupyter-execute::
+
+    def family_color(site):
+        return 'black' if site.family == a else 'white'
+
+    def hopping_lw(site1, site2):
+        return 0.04 if site1.family == site2.family else 0.1
+
+    kwant.plot(syst, site_lw=0.1, site_color=family_color, hop_lw=hopping_lw);
 
 Note that since we are using an unfinalized Builder, a ``site`` is really an
 instance of `~kwant.builder.Site`. With these adjustments we arrive at a plot
-that carries the same information, but is much easier to interpret:
-
-.. image:: /code/figure/plot_graphene_syst2.*
+that carries the same information, but is much easier to interpret.
 
 Apart from plotting the *system* itself, `~kwant.plotter.plot` can also be used
 to plot *data* living on the system.
@@ -67,23 +106,27 @@ nearest-neighbor hopping.  Computing the wave functions is done in the usual
 way (note that for a large-scale system, one would probably want to use sparse
 linear algebra):
 
-.. literalinclude:: /code/include/plot_graphene.py
-    :start-after: #HIDDEN_BEGIN_plotdata1
-    :end-before: #HIDDEN_END_plotdata1
 
-In most cases, to plot the wave function probability, one wouldn't use
-`~kwant.plotter.plot`, but rather `~kwant.plotter.map`. Here, we plot the
-`n`-th wave function using it:
+.. jupyter-execute::
+
+    import scipy.linalg as la
+
+    syst = make_system(tp=0).finalized()
+    ham = syst.hamiltonian_submatrix()
+    evecs = la.eigh(ham)[1]
 
-.. literalinclude:: /code/include/plot_graphene.py
-    :start-after: #HIDDEN_BEGIN_plotdata2
-    :end-before: #HIDDEN_END_plotdata2
+    wf = abs(evecs[:, 225])**2
 
+In most cases, to plot the wave function probability, one wouldn't use
+`~kwant.plotter.plot`, but rather `~kwant.plotter.map`. Here, we plot the
+`n`-th wave function using it.
 This results in a standard pseudocolor plot, showing in this case (``n=225``) a
 graphene edge state, i.e. a wave function mostly localized at the zigzag edges
 of the quantum dot.
 
-.. image:: /code/figure/plot_graphene_data1.*
+.. jupyter-execute::
+
+    kwant.plotter.map(syst, wf, oversampling=10, cmap='gist_heat_r');
 
 However although in general preferable, `~kwant.plotter.map` has a few
 deficiencies for this small system: For example, there are a few distortions at
@@ -99,9 +142,17 @@ the symbol shape depending on the sublattice. With a triangle pointing up and
 down on the respective sublattice, the symbols used by plot fill the space
 completely:
 
-.. literalinclude:: /code/include/plot_graphene.py
-    :start-after: #HIDDEN_BEGIN_plotdata3
-    :end-before: #HIDDEN_END_plotdata3
+.. jupyter-execute::
+
+    def family_shape(i):
+        site = syst.sites[i]
+        return ('p', 3, 180) if site.family == a else ('p', 3, 0)
+
+    def family_color(i):
+        return 'black' if syst.sites[i].family == a else 'white'
+
+    kwant.plot(syst, site_color=wf, site_symbol=family_shape,
+               site_size=0.5, hop_lw=0, cmap='gist_heat_r');
 
 Note that with ``hop_lw=0`` we deactivate plotting the hoppings (that would not
 serve any purpose here). Moreover, ``site_size=0.5`` guarantees that the two
@@ -114,31 +165,23 @@ Finally, note that since we are dealing with a finalized system now, a site `i`
 is represented by an integer. In order to obtain the original
 `~kwant.builder.Site`, ``syst.sites[i]`` can be used.
 
-With this we arrive at
-
-.. image:: /code/figure/plot_graphene_data2.*
-
-with the same information as `~kwant.plotter.map`, but with a cleaner look.
-
 The way how data is presented of course influences what features of the data
 are best visible in a given plot. With `~kwant.plotter.plot` one can easily go
 beyond pseudocolor-like plots. For example, we can represent the wave function
 probability using the symbols itself:
 
-.. literalinclude:: /code/include/plot_graphene.py
-    :start-after: #HIDDEN_BEGIN_plotdata4
-    :end-before: #HIDDEN_END_plotdata4
+.. jupyter-execute::
+
+    def site_size(i):
+        return 3 * wf[i] / wf.max()
+
+    kwant.plot(syst, site_size=site_size, site_color=(0, 0, 1, 0.3),
+               hop_lw=0.1);
 
 Here, we choose the symbol size proportional to the wave function probability,
 while the site color is transparent to also allow for overlapping symbols to be
 visible. The hoppings are also plotted in order to show the underlying lattice.
 
-With this, we arrive at
-
-.. image:: /code/figure/plot_graphene_data3.*
-
-which shows the edge state nature of the wave function most clearly.
-
 .. rubric:: Footnotes
 
 .. [#] A dangling site is defined as having only one hopping connecting it to