diff --git a/doc/source/pre/whatsnew/1.3.rst b/doc/source/pre/whatsnew/1.3.rst
index 0ba00ce81f0e7286d6319863662b5acacb1faabc..8336e99d9c66dc1fb7628414fe28cf648c18319c 100644
--- a/doc/source/pre/whatsnew/1.3.rst
+++ b/doc/source/pre/whatsnew/1.3.rst
@@ -52,6 +52,9 @@ a quantum anomalous spin Hall phase:
     import kwant
     import kwant.continuum
 
+.. jupyter-execute:: ../../tutorial/boilerplate.py
+    :hide-code:
+
 .. jupyter-execute::
 
     def make_model(a):
diff --git a/doc/source/tutorial/boilerplate.py b/doc/source/tutorial/boilerplate.py
new file mode 100644
index 0000000000000000000000000000000000000000..db2fba37a6f8546b66d4ca8f06a8eda668438ca7
--- /dev/null
+++ b/doc/source/tutorial/boilerplate.py
@@ -0,0 +1,6 @@
+import matplotlib
+import matplotlib.pyplot
+from IPython.display import set_matplotlib_formats
+
+matplotlib.rcParams['figure.figsize'] = matplotlib.pyplot.figaspect(1) * 2
+set_matplotlib_formats('svg')
diff --git a/doc/source/tutorial/discretize.rst b/doc/source/tutorial/discretize.rst
index 27bab68f61345c84d0e89a6c95a40cd5b5b5b176..31c4b3bb3ca217888c3c0e1b53e444a270906705 100644
--- a/doc/source/tutorial/discretize.rst
+++ b/doc/source/tutorial/discretize.rst
@@ -42,6 +42,9 @@ continuum models and for discretizing them into tight-binding models.
 
     import kwant
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 .. _tutorial_discretizer_introduction:
 
 Discretizing by hand
diff --git a/doc/source/tutorial/faq.rst b/doc/source/tutorial/faq.rst
index 42afab7953224b281072ba8d8b06832cf52873a8..ec52448d7afaf0d3ee4a2440088829a8dfc1b921 100644
--- a/doc/source/tutorial/faq.rst
+++ b/doc/source/tutorial/faq.rst
@@ -14,6 +14,9 @@ into Kwant's structure.
     import matplotlib
     from matplotlib import pyplot as plt
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 What is a system, and what is a builder?
 ========================================
 A Kwant system represents a particular tight-binding model. It contains a graph
diff --git a/doc/source/tutorial/first_steps.rst b/doc/source/tutorial/first_steps.rst
index bfb60e0c0123e1c409328b3fa3c6a06d2b463701..f509aaa0c262be17a1eaa53b8551812e06e88f22 100644
--- a/doc/source/tutorial/first_steps.rst
+++ b/doc/source/tutorial/first_steps.rst
@@ -88,6 +88,9 @@ In order to use Kwant, we need to import it:
     # For plotting
     from matplotlib import pyplot
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 .. jupyter-execute::
 
     import kwant
@@ -452,6 +455,11 @@ file and defining the a square lattice and empty scattering region.
     # For plotting
     from matplotlib import pyplot
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
+.. jupyter-execute::
+
     a = 1
     t = 1.0
     W, L = 10, 30
@@ -661,6 +669,11 @@ return a Kwant ``Builder``:
     from matplotlib import pyplot
     import kwant
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
+.. jupyter-execute::
+
     def make_system(L, W, a=1, t=1.0):
         lat = kwant.lattice.square(a)
 
diff --git a/doc/source/tutorial/graphene.rst b/doc/source/tutorial/graphene.rst
index 6c271d887911ea1faeb4356f2b12c1018e67497a..6eaa93d6ed142f430bd078399b2850cdea61cf32 100644
--- a/doc/source/tutorial/graphene.rst
+++ b/doc/source/tutorial/graphene.rst
@@ -47,6 +47,9 @@ explicitly here to show how to define a new lattice:
 
     sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 .. jupyter-execute::
 
     graphene = kwant.lattice.general([(1, 0), (sin_30, cos_30)],
diff --git a/doc/source/tutorial/kpm.rst b/doc/source/tutorial/kpm.rst
index 767423f32d899915d846f6ab3594c904ff437a1b..7e1232581335aa1a77969899f59344e34c93eba7 100644
--- a/doc/source/tutorial/kpm.rst
+++ b/doc/source/tutorial/kpm.rst
@@ -39,6 +39,9 @@ KPM method `kwant.kpm`, that is based on the algorithms presented in Ref. [1]_.
     import scipy
     from matplotlib import pyplot
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 Introduction
 ************
 
@@ -151,7 +154,7 @@ We start by importing kwant and defining our system.
 
     # Plot several density of states curves on the same axes.
     def plot_dos(labels_to_data):
-        pyplot.figure(figsize=(5,4))
+        pyplot.figure()
         for label, (x, y) in labels_to_data:
             pyplot.plot(x, y.real, label=label, linewidth=2)
         pyplot.legend(loc=2, framealpha=0.5)
@@ -162,7 +165,7 @@ We start by importing kwant and defining our system.
 
     # Plot fill density of states plus curves on the same axes.
     def plot_dos_and_curves(dos, labels_to_data):
-        pyplot.figure(figsize=(5,4))
+        pyplot.figure()
         pyplot.fill_between(dos[0], dos[1], label="DoS [a.u.]",
                          alpha=0.5, color='gray')
         for label, (x, y) in labels_to_data:
diff --git a/doc/source/tutorial/operators.rst b/doc/source/tutorial/operators.rst
index c10b7780ea70c4e735baa1e13e7d77fe415c50ec..1646ae86e2cd5a47747843f6dd952144b0a92175 100644
--- a/doc/source/tutorial/operators.rst
+++ b/doc/source/tutorial/operators.rst
@@ -51,6 +51,9 @@ texture.
     # letters denote spinor indices
     sigma = np.rollaxis(np.array([sigma_x, sigma_y, sigma_z]), 1)
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 Introduction
 ------------
 Our starting point will be the following spinful tight-binding model on
diff --git a/doc/source/tutorial/plotting.rst b/doc/source/tutorial/plotting.rst
index 0b84f555e7e93436731c46b17ddbda99e26e7e00..ea04094a7381839540b5c5dfe82ea35258ff4767 100644
--- a/doc/source/tutorial/plotting.rst
+++ b/doc/source/tutorial/plotting.rst
@@ -35,6 +35,9 @@ these options can be used to achieve various very different objectives.
 
     import kwant
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 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:
@@ -216,6 +219,9 @@ visible. The hoppings are also plotted in order to show the underlying lattice.
 
     import kwant
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 Zincblende is a very common crystal structure of semiconductors. It is a
 face-centered cubic crystal with two inequivalent atoms in the unit cell
 (i.e. two different types of atoms, unlike diamond which has the same crystal
diff --git a/doc/source/tutorial/spectrum.rst b/doc/source/tutorial/spectrum.rst
index 1eeef13528f797a0e982f9bb1ae276c458b46c7b..b2cbc51d3bef2b378891baca292857360ca6d95d 100644
--- a/doc/source/tutorial/spectrum.rst
+++ b/doc/source/tutorial/spectrum.rst
@@ -30,6 +30,9 @@ Band structure calculations
     # For plotting
     from matplotlib import pyplot
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 When doing transport simulations, one also often needs to know the band
 structure of the leads, i.e. the energies of the propagating plane waves in the
 leads as a function of momentum. This band structure contains information about
@@ -137,6 +140,9 @@ Closed systems
     from matplotlib import pyplot
     import kwant
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 Although Kwant is (currently) mainly aimed towards transport problems, it
 can also easily be used to compute properties of closed systems -- after
 all, a closed system is nothing more than a scattering region without leads!
diff --git a/doc/source/tutorial/spin_potential_shape.rst b/doc/source/tutorial/spin_potential_shape.rst
index 9107f701a6ec55c878082da4a5c916106760b60c..770ba4a55f06f0c4a8bf0b89669cb66bb2fad407 100644
--- a/doc/source/tutorial/spin_potential_shape.rst
+++ b/doc/source/tutorial/spin_potential_shape.rst
@@ -68,6 +68,9 @@ for small arrays.)
     # For plotting
     from matplotlib import pyplot
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 .. jupyter-execute::
 
     # For matrix support
@@ -242,6 +245,9 @@ Spatially dependent values through functions
     # For plotting
     from matplotlib import pyplot
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 Up to now, all examples had position-independent matrix-elements
 (and thus translational invariance along the wire, which
 was the origin of the conductance steps). Now, we consider the
@@ -403,6 +409,9 @@ Nontrivial shapes
     # For plotting
     from matplotlib import pyplot
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 Up to now, we only dealt with simple wire geometries. Now we turn to the case
 of a more complex geometry, namely transport through a quantum ring
 that is pierced by a magnetic flux :math:`\Phi`:
@@ -594,6 +603,12 @@ period of one flux quantum.
         import kwant
         from matplotlib import pyplot
 
+    .. jupyter-execute:: boilerplate.py
+        :hide-code:
+
+    .. jupyter-execute::
+        :hide-code:
+
         a = 1
         t = 1.0
         W = 10
@@ -643,6 +658,12 @@ period of one flux quantum.
         import kwant
         from matplotlib import pyplot
 
+    .. jupyter-execute:: boilerplate.py
+        :hide-code:
+
+    .. jupyter-execute::
+        :hide-code:
+
         a = 1
         t = 1.0
         W = 10
diff --git a/doc/source/tutorial/superconductors.rst b/doc/source/tutorial/superconductors.rst
index 56f5a849d2b7679e174ab02e9f8ec029b58efff1..78f004d34820c3c4ef6219a855681ab542f0aa7a 100644
--- a/doc/source/tutorial/superconductors.rst
+++ b/doc/source/tutorial/superconductors.rst
@@ -35,6 +35,9 @@ Superconductors: orbital degrees of freedom, conservation laws and symmetries
     tau_y = tinyarray.array([[0, -1j], [1j, 0]])
     tau_z = tinyarray.array([[1, 0], [0, -1]])
 
+.. jupyter-execute:: boilerplate.py
+    :hide-code:
+
 This example deals with superconductivity on the level of the
 Bogoliubov-de Gennes (BdG) equation. In this framework, the Hamiltonian
 is given as