diff --git a/.gitignore b/.gitignore
index 4493ae39ec3787174e65a72cff84e750a739b940..6687006404604ac19439b0204412f88a12e61acc 100644
--- a/.gitignore
+++ b/.gitignore
@@ -5,6 +5,7 @@ MANIFEST
 /kwant/*.c
 /kwant/*/*.c
 /kwant/_static_version.py
+/tutorial
 /build
 /dist
 /doc/build
diff --git a/MANIFEST.in b/MANIFEST.in
index 64fcabc7caae394c5e478ace1d9d54fd489dfa7e..59f2d474632eeccc51d1cafb9b19da479af67868 100644
--- a/MANIFEST.in
+++ b/MANIFEST.in
@@ -12,9 +12,7 @@ recursive-include kwant *.pxd
 recursive-include kwant *.c
 recursive-include kwant *.h
 recursive-include kwant test_*.py
-
 recursive-include examples *.py
-recursive-include tutorial *.py
 
 include doc/other/*[a-zA-Z]
 include doc/Makefile
diff --git a/doc/Makefile b/doc/Makefile
index 3de1f635e7c90a57aa649306c4df2c4c0f2fa800..36a2cf09d1fa68cc17b2c4f2707b0acf0ce77238 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -119,8 +119,9 @@ doctest:
 
 
 # Make the image generation scripts by patching tutorial scipts.
+.SECONDARY:
 %.py: %.py.diff
-	@cp ../tutorial/$(notdir $@) $(dir $@)
+	@grep -v '^#HIDDEN' source/tutorial/$(notdir $@) >$@
 	@patch $@ $<
 
 # Generation of tutorial images.  This requires some make trickery, see
diff --git a/doc/source/images/generate-diffs.sh b/doc/source/images/generate-diffs.sh
index e68caceffb22fabb17a1e8da57bfd10f98acb776..aeeb479c44112754ccebf331cbef638f519b92db 100755
--- a/doc/source/images/generate-diffs.sh
+++ b/doc/source/images/generate-diffs.sh
@@ -6,5 +6,6 @@
 for f in [0-9]-*.py; do
     # We use custom labels to suppress the time stamps which are unnecessary
     # here and would only lead to noise in version control.
-    diff -u --label original --label modified ../../../tutorial/$f $f >$f.diff
+    grep -v '#HIDDEN' ../tutorial/$f |
+    diff -u --label original --label modified - $f >$f.diff
 done
diff --git a/tutorial/1-quantum_wire.py b/doc/source/tutorial/1-quantum_wire.py
similarity index 83%
rename from tutorial/1-quantum_wire.py
rename to doc/source/tutorial/1-quantum_wire.py
index 02935c56ec61b529b93b8e6fdc93b8dd9de719bd..9e8ad75a5d3ea08f54e337dc229327073c0b56e0 100644
--- a/tutorial/1-quantum_wire.py
+++ b/doc/source/tutorial/1-quantum_wire.py
@@ -8,16 +8,23 @@
 #  - Making scattering region and leads
 #  - Using the simple sparse solver for computing Landauer conductance
 
+#HIDDEN_BEGIN_dwhx
 import kwant
+#HIDDEN_END_dwhx
 
 # First, define the tight-binding system
 
+#HIDDEN_BEGIN_goiq
 sys = kwant.Builder()
+#HIDDEN_END_goiq
 
 # Here, we are only working with square lattices
+#HIDDEN_BEGIN_suwo
 a = 1
 lat = kwant.lattice.Square(a)
+#HIDDEN_END_suwo
 
+#HIDDEN_BEGIN_zfvr
 t = 1.0
 W = 10
 L = 30
@@ -35,6 +42,7 @@ for i in xrange(L):
         #hopping in x-direction
         if i > 0:
             sys[lat(i, j), lat(i - 1, j)] = -t
+#HIDDEN_END_zfvr
 
 # Then, define the leads:
 
@@ -42,9 +50,12 @@ for i in xrange(L):
 
 # (Note: in the current version, TranslationalSymmetry takes a
 # realspace vector)
+#HIDDEN_BEGIN_xcmc
 sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
 lead0 = kwant.Builder(sym_lead0)
+#HIDDEN_END_xcmc
 
+#HIDDEN_BEGIN_ndez
 for j in xrange(W):
     lead0[lat(0, j)] = 4 * t
 
@@ -52,8 +63,10 @@ for j in xrange(W):
         lead0[lat(0, j), lat(0, j - 1)] = -t
 
     lead0[lat(1, j), lat(0, j)] = -t
+#HIDDEN_END_ndez
 
 # Then the lead to the right
+#HIDDEN_BEGIN_xhqc
 
 sym_lead1 = kwant.TranslationalSymmetry([lat.vec((1, 0))])
 lead1 = kwant.Builder(sym_lead1)
@@ -65,22 +78,30 @@ for j in xrange(W):
         lead1[lat(0, j), lat(0, j - 1)] = -t
 
     lead1[lat(1, j), lat(0, j)] = -t
+#HIDDEN_END_xhqc
 
 # Then attach the leads to the system
 
+#HIDDEN_BEGIN_fskr
 sys.attach_lead(lead0)
 sys.attach_lead(lead1)
+#HIDDEN_END_fskr
 
 # Plot it, to make sure it's OK
 
+#HIDDEN_BEGIN_wsgh
 kwant.plot(sys)
+#HIDDEN_END_wsgh
 
 # Finalize the system
 
+#HIDDEN_BEGIN_dngj
 sys = sys.finalized()
+#HIDDEN_END_dngj
 
 # Now that we have the system, we can compute conductance
 
+#HIDDEN_BEGIN_buzn
 energies = []
 data = []
 for ie in xrange(100):
@@ -93,9 +114,11 @@ for ie in xrange(100):
     # lead 1
     energies.append(energy)
     data.append(smatrix.transmission(1, 0))
+#HIDDEN_END_buzn
 
 # Use matplotlib to write output
 # We should see conductance steps
+#HIDDEN_BEGIN_lliv
 from matplotlib import pyplot
 
 pyplot.figure()
@@ -103,3 +126,4 @@ pyplot.plot(energies, data)
 pyplot.xlabel("energy [in units of t]")
 pyplot.ylabel("conductance [in units of e^2/h]")
 pyplot.show()
+#HIDDEN_END_lliv
diff --git a/tutorial/1-quantum_wire_revisited.py b/doc/source/tutorial/1-quantum_wire_revisited.py
similarity index 87%
rename from tutorial/1-quantum_wire_revisited.py
rename to doc/source/tutorial/1-quantum_wire_revisited.py
index 8f3202a5b4c098d887dd9756ef6fd77935f6dec6..2e7da1aabef6fafdc47fb19f41c179c90623d363 100644
--- a/tutorial/1-quantum_wire_revisited.py
+++ b/doc/source/tutorial/1-quantum_wire_revisited.py
@@ -10,6 +10,7 @@
 # Note: Does the same as tutorial1a.py, but using other features of kwant
 #
 
+#HIDDEN_BEGIN_xkzy
 import kwant
 
 # For plotting
@@ -22,33 +23,45 @@ def make_system(a=1, t=1.0, W=10, L=30):
     lat = kwant.lattice.Square(a)
 
     sys = kwant.Builder()
+#HIDDEN_END_xkzy
 
     #### Define the scattering region. ####
+#HIDDEN_BEGIN_vvjt
     sys[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
+#HIDDEN_END_vvjt
+#HIDDEN_BEGIN_nooi
     for hopping in lat.nearest:
         sys[sys.possible_hoppings(*hopping)] = -t
+#HIDDEN_END_nooi
 
     #### Define the leads. ####
     # First the lead to the left, ...
     # (Note: in the current version, TranslationalSymmetry takes a
     # realspace vector)
+#HIDDEN_BEGIN_iepx
     sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead0 = kwant.Builder(sym_lead0)
 
     lead0[(lat(0, j) for j in xrange(W))] = 4 * t
     for hopping in lat.nearest:
         lead0[lead0.possible_hoppings(*hopping)] = -t
+#HIDDEN_END_iepx
 
     # ... then the lead to the right.  We use a method that returns a copy of
     # `lead0` with its direction reversed.
+#HIDDEN_BEGIN_xkdo
     lead1 = lead0.reversed()
+#HIDDEN_END_xkdo
 
     #### Attach the leads and return the system. ####
+#HIDDEN_BEGIN_yxot
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
 
     return sys
 
+#HIDDEN_END_yxot
+#HIDDEN_BEGIN_ayuk
 
 def plot_conductance(sys, energies):
     # Compute conductance
@@ -62,8 +75,10 @@ def plot_conductance(sys, energies):
     pyplot.xlabel("energy [in units of t]")
     pyplot.ylabel("conductance [in units of e^2/h]")
     pyplot.show()
+#HIDDEN_END_ayuk
 
 
+#HIDDEN_BEGIN_cjel
 def main():
     sys = make_system()
 
@@ -75,9 +90,12 @@ def main():
 
     # We should see conductance steps.
     plot_conductance(sys, energies=[0.01 * i for i in xrange(100)])
+#HIDDEN_END_cjel
 
 
 # Call the main function if the script gets executed (as opposed to imported).
 # See <http://docs.python.org/library/__main__.html>.
+#HIDDEN_BEGIN_ypbj
 if __name__ == '__main__':
     main()
+#HIDDEN_END_ypbj
diff --git a/tutorial/2-ab_ring.py b/doc/source/tutorial/2-ab_ring.py
similarity index 95%
rename from tutorial/2-ab_ring.py
rename to doc/source/tutorial/2-ab_ring.py
index 16559d6ab362b0ff93769cac74866ca86e6f2ddd..52c08259bd2f5f1239d2e760c65aabdd94185d99 100644
--- a/tutorial/2-ab_ring.py
+++ b/doc/source/tutorial/2-ab_ring.py
@@ -18,6 +18,7 @@ import kwant
 from matplotlib import pyplot
 
 
+#HIDDEN_BEGIN_eusz
 def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
     # Start with an empty tight-binding system and a single square lattice.
     # `a` is the lattice constant (by default set to 1 for simplicity).
@@ -32,11 +33,14 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
         (x, y) = pos
         rsq = x ** 2 + y ** 2
         return (r1 ** 2 < rsq < r2 ** 2)
+#HIDDEN_END_eusz
 
     # and add the corresponding lattice points using the `shape`-function
+#HIDDEN_BEGIN_lcak
     sys[lat.shape(ring, (0, r1 + 1))] = 4 * t
     for hopping in lat.nearest:
         sys[sys.possible_hoppings(*hopping)] = -t
+#HIDDEN_END_lcak
 
     # In order to introduce a flux through the ring, we introduce a phase
     # on the hoppings on the line cut through one of the arms
@@ -44,6 +48,7 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
     # since we want to change the flux without modifying Builder repeatedly,
     # we define the modified hoppings as a function that takes the flux
     # through the global variable phi.
+#HIDDEN_BEGIN_lvkt
     def fluxphase(site1, site2):
         return exp(1j * phi)
 
@@ -57,9 +62,11 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
     # Modify only those hopings in x-direction that cross the branch cut
     sys[(hop for hop in sys.possible_hoppings((1, 0), lat, lat)
          if crosses_branchcut(hop))] = fluxphase
+#HIDDEN_END_lvkt
 
     #### Define the leads. ####
     # left lead
+#HIDDEN_BEGIN_qwgr
     sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead0 = kwant.Builder(sym_lead0)
 
@@ -70,14 +77,17 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
     lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
     for hopping in lat.nearest:
         lead0[lead0.possible_hoppings(*hopping)] = -t
+#HIDDEN_END_qwgr
 
     # Then the lead to the right
     # [again, obtained using reversed()]
     lead1 = lead0.reversed()
 
     #### Attach the leads and return the system. ####
+#HIDDEN_BEGIN_skbz
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
+#HIDDEN_END_skbz
 
     return sys
 
diff --git a/tutorial/2-quantum_well.py b/doc/source/tutorial/2-quantum_well.py
similarity index 90%
rename from tutorial/2-quantum_well.py
rename to doc/source/tutorial/2-quantum_well.py
index d51ba0077dc47d0e5d3921caa895808ef88bf80b..5bbead7d95e082dff37593850c2b64c57703a239 100644
--- a/tutorial/2-quantum_well.py
+++ b/doc/source/tutorial/2-quantum_well.py
@@ -13,6 +13,10 @@ from matplotlib import pyplot
 
 # global variable governing the behavior of potential() in
 # make_system()
+#HIDDEN The following code line is included verbatim in the tutorial text
+#HIDDEN because nested code examples are not supported.  Remember to update
+#HIDDEN the tutorial text when you modify this line.
+#HIDDEN_BEGIN_ehso
 pot = 0
 
 def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
@@ -31,13 +35,16 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
             return pot
         else:
             return 0
+#HIDDEN_END_ehso
 
+#HIDDEN_BEGIN_coid
     def onsite(site):
         return 4 * t + potential(site)
 
     sys[(lat(x, y) for x in range(L) for y in range(W))] = onsite
     for hopping in lat.nearest:
         sys[sys.possible_hoppings(*hopping)] = -t
+#HIDDEN_END_coid
 
     #### Define the leads. ####
     # First the lead to the left, ...
@@ -64,6 +71,7 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
 def plot_conductance(sys, energy, welldepths):
     # We specify that we want to not only read, but also write to a
     # global variable.
+#HIDDEN_BEGIN_sqvr
     global pot
 
     # Compute conductance
@@ -80,6 +88,7 @@ def plot_conductance(sys, energy, welldepths):
     pyplot.xlabel("well depth [in units of t]")
     pyplot.ylabel("conductance [in units of e^2/h]")
     pyplot.show()
+#HIDDEN_END_sqvr
 
 
 def main():
diff --git a/tutorial/2-spin_orbit.py b/doc/source/tutorial/2-spin_orbit.py
similarity index 95%
rename from tutorial/2-spin_orbit.py
rename to doc/source/tutorial/2-spin_orbit.py
index 009132b21f00674856445203deaa7ad692a6397e..a37e5acc22eaeec84aa5f5cd782c697f452ca299 100644
--- a/tutorial/2-spin_orbit.py
+++ b/doc/source/tutorial/2-spin_orbit.py
@@ -16,13 +16,17 @@ import kwant
 from matplotlib import pyplot
 
 # For matrix support
+#HIDDEN_BEGIN_xumz
 import numpy
+#HIDDEN_END_xumz
 
 # define Pauli-matrices for convenience
+#HIDDEN_BEGIN_hwbt
 sigma_0 = numpy.eye(2)
 sigma_x = numpy.array([[0, 1], [1, 0]])
 sigma_y = numpy.array([[0, -1j], [1j, 0]])
 sigma_z = numpy.array([[1, 0], [0, -1]])
+#HIDDEN_END_hwbt
 
 
 def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
@@ -33,6 +37,7 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
     sys = kwant.Builder()
 
     #### Define the scattering region. ####
+#HIDDEN_BEGIN_uxrm
     sys[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t * sigma_0 + \
         e_z * sigma_z
     # hoppings in x-direction
@@ -41,12 +46,14 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
     # hoppings in y-directions
     sys[sys.possible_hoppings((0, 1), lat, lat)] = -t * sigma_0 + \
         1j * alpha * sigma_x
+#HIDDEN_END_uxrm
 
     #### Define the leads. ####
     # left lead
     sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead0 = kwant.Builder(sym_lead0)
 
+#HIDDEN_BEGIN_yliu
     lead0[(lat(0, j) for j in xrange(W))] = 4 * t * sigma_0 + e_z * sigma_z
     # hoppings in x-direction
     lead0[lead0.possible_hoppings((1, 0), lat, lat)] = -t * sigma_0 - \
@@ -54,6 +61,7 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
     # hoppings in y-directions
     lead0[lead0.possible_hoppings((0, 1), lat, lat)] = -t * sigma_0 + \
         1j * alpha * sigma_x
+#HIDDEN_END_yliu
 
     # Then the lead to the right
     # (again, obtained using reverse()
diff --git a/tutorial/3-band_structure.py b/doc/source/tutorial/3-band_structure.py
similarity index 95%
rename from tutorial/3-band_structure.py
rename to doc/source/tutorial/3-band_structure.py
index 8a17633c104dd21b1d94ef8f802111b066369675..3b8fe8e0182ef35c790eacf8d578e80a10271144 100644
--- a/tutorial/3-band_structure.py
+++ b/doc/source/tutorial/3-band_structure.py
@@ -14,6 +14,7 @@ from math import pi
 from matplotlib import pyplot
 
 
+#HIDDEN_BEGIN_zxip
 def make_lead(a=1, t=1.0, W=10):
     # Start with an empty lead with a single square lattice
     lat = kwant.lattice.Square(a)
@@ -32,8 +33,10 @@ def make_lead(a=1, t=1.0, W=10):
         lead[lat(1, j), lat(0, j)] = -t
 
     return lead
+#HIDDEN_END_zxip
 
 
+#HIDDEN_BEGIN_pejz
 def plot_bandstructure(lead, momenta):
     # Use the method ``energies`` of the finalized lead to compute
     # the bandstructure
@@ -53,6 +56,7 @@ def main():
     momenta = [-pi + 0.02 * pi * i for i in xrange(101)]
 
     plot_bandstructure(lead, momenta)
+#HIDDEN_END_pejz
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/tutorial/3-closed_system.py b/doc/source/tutorial/3-closed_system.py
similarity index 95%
rename from tutorial/3-closed_system.py
rename to doc/source/tutorial/3-closed_system.py
index 82361ebee728924816d91a11f7fa9391b3894092..97cd7405ef7981fb7d6b9599d70f979a65432ef9 100644
--- a/tutorial/3-closed_system.py
+++ b/doc/source/tutorial/3-closed_system.py
@@ -13,7 +13,9 @@ from cmath import exp
 import kwant
 
 # For eigenvalue computation
+#HIDDEN_BEGIN_tibv
 import scipy.linalg as la
+#HIDDEN_END_tibv
 
 # For plotting
 from matplotlib import pyplot
@@ -23,6 +25,7 @@ def make_system(a=1, t=1.0, r=10):
     # Start with an empty tight-binding system and a single square lattice.
     # `a` is the lattice constant (by default set to 1 for simplicity).
 
+#HIDDEN_BEGIN_qlyd
     lat = kwant.lattice.Square(a)
 
     sys = kwant.Builder()
@@ -46,8 +49,10 @@ def make_system(a=1, t=1.0, r=10):
 
     # It's a closed system for a change, so no leads
     return sys
+#HIDDEN_END_qlyd
 
 
+#HIDDEN_BEGIN_yvri
 def plot_spectrum(sys, Bfields):
     # global variable B controls the magnetic field
     global B
@@ -74,6 +79,7 @@ def plot_spectrum(sys, Bfields):
     pyplot.xlabel("magnetic field [some arbitrary units]")
     pyplot.ylabel("energy [in units of t]")
     pyplot.show()
+#HIDDEN_END_yvri
 
 
 def main():
diff --git a/tutorial/4-graphene.py b/doc/source/tutorial/4-graphene.py
similarity index 88%
rename from tutorial/4-graphene.py
rename to doc/source/tutorial/4-graphene.py
index 95782a819d213beec71f6f9ef89fc02cf28143b6..635dedbf9e238eb2ec78e7a439667522e4c3c833 100644
--- a/tutorial/4-graphene.py
+++ b/doc/source/tutorial/4-graphene.py
@@ -21,11 +21,14 @@ from matplotlib import pyplot
 
 # Define the graphene lattice
 sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
+#HIDDEN_BEGIN_hnla
 graphene = kwant.make_lattice([(1, 0), (sin_30, cos_30)],
                               [(0, 0), (0, 1 / sqrt(3))])
 a, b = graphene.sublattices
 
+#HIDDEN_END_hnla
 
+#HIDDEN_BEGIN_shzy
 def make_system(r=10, w=2.0, pot=0.1):
 
     #### Define the scattering region. ####
@@ -43,18 +46,26 @@ def make_system(r=10, w=2.0, pot=0.1):
         return pot * tanh(d / w)
 
     sys[graphene.shape(circle, (0, 0))] = potential
+#HIDDEN_END_shzy
 
     # specify the hoppings of the graphene lattice in the
     # format expected by possibe_hoppings()
+#HIDDEN_BEGIN_hsmc
     hoppings = (((0, 0), b, a), ((0, 1), b, a), ((-1, 1), b, a))
+#HIDDEN_END_hsmc
+#HIDDEN_BEGIN_bfwb
     for hopping in hoppings:
         sys[sys.possible_hoppings(*hopping)] = -1
+#HIDDEN_END_bfwb
 
     # Modify the scattering region
+#HIDDEN_BEGIN_efut
     del sys[a(0, 0)]
     sys[a(-2, 1), b(2, 2)] = -1
+#HIDDEN_END_efut
 
     #### Define the leads. ####
+#HIDDEN_BEGIN_aakh
     # left lead
     sym0 = kwant.TranslationalSymmetry([graphene.vec((-1, 0))])
 
@@ -80,10 +91,14 @@ def make_system(r=10, w=2.0, pot=0.1):
     lead1[graphene.shape(lead1_shape, (0, 0))] = pot
     for hopping in hoppings:
         lead1[lead1.possible_hoppings(*hopping)] = -1
+#HIDDEN_END_aakh
 
+#HIDDEN_BEGIN_kmmw
     return sys, [lead0, lead1]
+#HIDDEN_END_kmmw
 
 
+#HIDDEN_BEGIN_zydk
 def compute_evs(sys):
     # Compute some eigenvalues of the closed system
     sparse_mat = sys.hamiltonian_submatrix(sparse=True)
@@ -96,6 +111,7 @@ def compute_evs(sys):
         print evs
     except:
         pass
+#HIDDEN_END_zydk
 
 
 def plot_conductance(sys, energies):
@@ -124,6 +140,11 @@ def plot_bandstructure(flead, momenta):
     pyplot.show()
 
 
+#HIDDEN The part of the following code block which begins with plotter_symbols
+#HIDDEN is included verbatim in the tutorial text because nested code examples
+#HIDDEN are not supported.  Remember to update the tutorial text when you
+#HIDDEN modify this block.
+#HIDDEN_BEGIN_itkk
 def main():
     pot = 0.1
     sys, leads = make_system(pot=pot)
@@ -137,9 +158,12 @@ def main():
 
     # Plot the closed system without leads.
     kwant.plot(sys, symbols=plotter_symbols)
+#HIDDEN_END_itkk
 
     # Compute some eigenvalues.
+#HIDDEN_BEGIN_jmbi
     compute_evs(sys.finalized())
+#HIDDEN_END_jmbi
 
     # Attach the leads to the system.
     for lead in leads:
diff --git a/tutorial/5-superconductor_band_structure.py b/doc/source/tutorial/5-superconductor_band_structure.py
similarity index 98%
rename from tutorial/5-superconductor_band_structure.py
rename to doc/source/tutorial/5-superconductor_band_structure.py
index 27e39d3d298e9bf1cbed44b8cdbe940481160d42..a0dca9e708c30890d34cc0ca4916a98046f35687 100644
--- a/tutorial/5-superconductor_band_structure.py
+++ b/doc/source/tutorial/5-superconductor_band_structure.py
@@ -18,6 +18,7 @@ from math import pi
 # For plotting
 from matplotlib import pyplot
 
+#HIDDEN_BEGIN_nbvn
 tau_x = np.array([[0, 1], [1, 0]])
 tau_z = np.array([[1, 0], [0, -1]])
 
@@ -40,6 +41,7 @@ def make_lead(a=1, t=1.0, mu=0.7, Delta=0.1, W=10):
         lead[lat(1, j), lat(0, j)] = -t * tau_z
 
     return lead
+#HIDDEN_END_nbvn
 
 
 def plot_bandstructure(lead, momenta):
diff --git a/tutorial/5-superconductor_transport.py b/doc/source/tutorial/5-superconductor_transport.py
similarity index 94%
rename from tutorial/5-superconductor_transport.py
rename to doc/source/tutorial/5-superconductor_transport.py
index 977fb1e121960ced4562f25725b788b1103cb77b..1e0e3bb3f51ad418453311947b14c736808aea93 100644
--- a/tutorial/5-superconductor_transport.py
+++ b/doc/source/tutorial/5-superconductor_transport.py
@@ -13,13 +13,16 @@ import kwant
 from matplotlib import pyplot
 
 
+#HIDDEN_BEGIN_zuuw
 def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
                 mu=0.4, Delta=0.1, Deltapos=4, t=1.0):
     # Start with an empty tight-binding system and two square lattices,
     # corresponding to electron and hole degree of freedom
     lat_e = kwant.lattice.Square(a)
     lat_h = kwant.lattice.Square(a)
+#HIDDEN_END_zuuw
 
+#HIDDEN_BEGIN_pqmp
     sys = kwant.Builder()
 
     #### Define the scattering region. ####
@@ -42,8 +45,10 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
     # electrons and holes
     sys[((lat_e(x, y), lat_h(x, y)) for x in range(Deltapos, L)
          for y in range(W))] = Delta
+#HIDDEN_END_pqmp
 
     #### Define the leads. ####
+#HIDDEN_BEGIN_ttth
     # left electron lead
     sym_lead0 = kwant.TranslationalSymmetry([lat_e.vec((-1, 0))])
     lead0 = kwant.Builder(sym_lead0)
@@ -61,10 +66,12 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
     # hoppings in x and y-direction
     lead1[lead1.possible_hoppings((1, 0), lat_h, lat_h)] = t
     lead1[lead1.possible_hoppings((0, 1), lat_h, lat_h)] = t
+#HIDDEN_END_ttth
 
     # Then the lead to the right
     # this one is superconducting and thus is comprised of electrons
     # AND holes
+#HIDDEN_BEGIN_mhiw
     sym_lead2 = kwant.TranslationalSymmetry([lat_e.vec((1, 0))])
     lead2 = kwant.Builder(sym_lead2)
 
@@ -76,15 +83,19 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
     lead2[lead2.possible_hoppings((1, 0), lat_h, lat_h)] = t
     lead2[lead2.possible_hoppings((0, 1), lat_h, lat_h)] = t
     lead2[((lat_e(0, j), lat_h(0, j)) for j in xrange(W))] = Delta
+#HIDDEN_END_mhiw
 
     #### Attach the leads and return the system. ####
+#HIDDEN_BEGIN_ozsr
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
     sys.attach_lead(lead2)
 
     return sys
+#HIDDEN_END_ozsr
 
 
+#HIDDEN_BEGIN_jbjt
 def plot_conductance(sys, energies):
     # Compute conductance
     data = []
@@ -94,6 +105,7 @@ def plot_conductance(sys, energies):
         data.append(smatrix.submatrix(0, 0).shape[0] -
                     smatrix.transmission(0, 0) +
                     smatrix.transmission(1, 0))
+#HIDDEN_END_jbjt
 
     pyplot.figure()
     pyplot.plot(energies, data)
diff --git a/tutorial/README b/doc/source/tutorial/README
similarity index 91%
rename from tutorial/README
rename to doc/source/tutorial/README
index d77eef72386f20d7f22184cb427b4239f5cf6ddf..b7174b90ca4f4b03d20930042b58a4e9182daeb3 100644
--- a/tutorial/README
+++ b/doc/source/tutorial/README
@@ -1,6 +1,3 @@
-This directory contains the example scripts of the tutorial.
-
-
 Note for kwant developers
 -------------------------
 
diff --git a/doc/source/tutorial/tutorial1.rst b/doc/source/tutorial/tutorial1.rst
index 7e02b6a839d0811f962f52761992cfedf3bc394b..3cbc64e2d2d2f57cfd8b262d5d4b63d88be7ce94 100644
--- a/doc/source/tutorial/tutorial1.rst
+++ b/doc/source/tutorial/tutorial1.rst
@@ -16,8 +16,9 @@ with a hard wall confinement :math:`V(y)` in y-direction.
 
 In order to use kwant, we need to import it:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 11
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_dwhx
+    :end-before: #HIDDEN_END_dwhx
 
 Enabling kwant is as easy as this [#]_ !
 
@@ -26,15 +27,17 @@ and leads. For this we make use of the `~kwant.builder.Builder` class
 that allows for a convenient way to define the system. For this we need to
 create an instance of the `~kwant.builder.Builder` class:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 15
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_goiq
+    :end-before: #HIDDEN_END_goiq
 
 Apart from `~kwant.builder.Builder` we also need to specify
 what kind of sites we want to add to the system. Here we work with
 a square lattice. For simplicity, we set the lattice constant to unity:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 18-19
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_suwo
+    :end-before: #HIDDEN_END_suwo
 
 Since we work with a square lattice, we label the points with two
 integer coordinates `(i, j)`. `~kwant.builder.Builder` then
@@ -51,15 +54,17 @@ needed in Builder (more about that in the technical details below).
 We now build a rectangular scattering region that is `W`
 lattice points wide and `L` lattice points long:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 21-23, 26-37
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_zfvr
+    :end-before: #HIDDEN_END_zfvr
 
 Next, we define the leads. Leads are also constructed using
 `~kwant.builder.Builder`, but in this case, the
 system must have a translational symmetry:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 45-46
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_xcmc
+    :end-before: #HIDDEN_END_xcmc
 
 Here, the `~kwant.builder.Builder` takes the translational symmetry
 as an optional parameter. Note that the (real space)
@@ -73,8 +78,9 @@ as the hoppings inside one unit cell and to the next unit cell of the lead.
 For a square lattice, and a lead in y-direction the unit cell is
 simply a vertical line of points:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 48-54
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_ndez
+    :end-before: #HIDDEN_END_ndez
 
 Note that here it doesn't matter if you add the hoppings to the next or the
 previous unit cell -- the translational symmetry takes care of that.
@@ -83,8 +89,9 @@ We also want to add a lead on the right side. The only difference to
 the left lead is that the vector of the translational
 symmetry must point to the right, the remaining code is the same:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 57-67
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_xhqc
+    :end-before: #HIDDEN_END_xhqc
 
 Note that here we added points with x-coordinate 0, just as for the left lead.
 You might object that the right lead should be placed `L`
@@ -94,8 +101,9 @@ you do not need to worry about that. The `~kwant.builder.Builder` with
 infinitely extended. These isolated, infinite leads can then be simply
 attached at the right position using:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 71-72
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_fskr
+    :end-before: #HIDDEN_END_fskr
 
 More details about attaching leads can be found in the tutorial
 :ref:`tutorial-abring`.
@@ -103,8 +111,9 @@ More details about attaching leads can be found in the tutorial
 Now we have finished building our system! We plot it, to make sure we didn't
 make any mistakes:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 76
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_wsgh
+    :end-before: #HIDDEN_END_wsgh
 
 This should bring up this picture:
 
@@ -118,14 +127,16 @@ fading color.
 
 In order to use our system for a transport calculation, we need to finalize it
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 80
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_dngj
+    :end-before: #HIDDEN_END_dngj
 
 Having successfully created a system, we now can immediately start to compute
 its conductance as a function of energy:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 84-95
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_buzn
+    :end-before: #HIDDEN_END_buzn
 
 Currently **CHANGE**, there is only one algorithm implemented to compute the
 conductance: :func:`kwant.solve <kwant.solvers.common.SparseSolver.solve>`
@@ -138,8 +149,9 @@ Finally we can use `matplotlib` to make a plot of the computed data
 (although writing to file and using an external viewer such as
 gnuplot or xmgrace is just as viable)
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire.py
-    :lines: 99-105
+.. literalinclude:: 1-quantum_wire.py
+    :start-after: #HIDDEN_BEGIN_lliv
+    :end-before: #HIDDEN_END_lliv
 
 This should yield the result
 
@@ -217,9 +229,9 @@ subbands that increases with energy.
 
      ::
 
-	 fsys = sys.finalized()
-	 del sys
-	 sys = fsys
+         fsys = sys.finalized()
+         del sys
+         sys = fsys
 
    - Note that the vector passed to the `~kwant.lattice.TranslationalSymmetry`
      (in fact, what is passed is a list of vectors -- there could be more than
@@ -267,14 +279,16 @@ We begin the program collecting all imports in the beginning of the
 file and put the build-up of the system into a separate function
 `make_system`:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
-    :lines: 13-24
+.. literalinclude:: 1-quantum_wire_revisited.py
+    :start-after: #HIDDEN_BEGIN_xkzy
+    :end-before: #HIDDEN_END_xkzy
 
 Previously, the scattering region was build using two ``for``-loops.
 Instead, we now write:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
-    :lines: 27
+.. literalinclude:: 1-quantum_wire_revisited.py
+    :start-after: #HIDDEN_BEGIN_vvjt
+    :end-before: #HIDDEN_END_vvjt
 
 Here, all lattice points are added at once in the first line.  The
 construct ``((i, j) for i in xrange(L) for j in xrange(W))`` is a
@@ -290,8 +304,9 @@ hoppings. In this case, an iterable like for the lattice
 points becomes a bit cumbersome, and we use instead another
 feature of kwant:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
-    :lines: 28-29
+.. literalinclude:: 1-quantum_wire_revisited.py
+    :start-after: #HIDDEN_BEGIN_nooi
+    :end-before: #HIDDEN_END_nooi
 
 In regular lattices, one has only very few types of different hoppings
 (by one lattice point in x or y-direction in the case of a square
@@ -308,8 +323,9 @@ then sets all of those hopping matrix elements at once.
 
 The leads can be constructed in an analogous way:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
-    :lines: 35-40
+.. literalinclude:: 1-quantum_wire_revisited.py
+    :start-after: #HIDDEN_BEGIN_iepx
+    :end-before: #HIDDEN_END_iepx
 
 Note that in the previous example, we essentially used the same code
 for the right and the left lead, the only difference was the direction
@@ -320,34 +336,39 @@ lead, but with it's translational vector reversed.  This can thus be
 used to obtain a lead pointing in the opposite direction, i.e. makes a
 right lead from a left lead:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
-    :lines: 44
+.. literalinclude:: 1-quantum_wire_revisited.py
+    :start-after: #HIDDEN_BEGIN_xkdo
+    :end-before: #HIDDEN_END_xkdo
 
 The remainder of the code is identical to the previous example
 (except for a bit of reorganization into functions):
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
-    :lines: 47-51
+.. literalinclude:: 1-quantum_wire_revisited.py
+    :start-after: #HIDDEN_BEGIN_yxot
+    :end-before: #HIDDEN_END_yxot
 
 and
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
-    :lines: 52-64
+.. literalinclude:: 1-quantum_wire_revisited.py
+    :start-after: #HIDDEN_BEGIN_ayuk
+    :end-before: #HIDDEN_END_ayuk
 
 Finally, we use a python trick to make our example usable both
 as a script, as well as allowing it to be imported as a module.
 We collect all statements that should be executed in the script
 in a ``main()``-function:
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
-    :lines: 67-77
+.. literalinclude:: 1-quantum_wire_revisited.py
+    :start-after: #HIDDEN_BEGIN_cjel
+    :end-before: #HIDDEN_END_cjel
 
 Finally, we use the following python construct [#]_ that executes
 ``main()`` if the program is used as a script (i.e. executed as
 ``python tutorial1b.py``):
 
-.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
-    :lines: 82-83
+.. literalinclude:: 1-quantum_wire_revisited.py
+    :start-after: #HIDDEN_BEGIN_ypbj
+    :end-before: #HIDDEN_END_ypbj
 
 If the example however is imported using ``import tutorial1b``,
 ``main()`` is not executed automatically. Instead, you can execute it
@@ -365,8 +386,9 @@ The result of the example should be identical to the previous one.
 
    - In
 
-     .. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
-       :lines: 28-29
+     .. literalinclude:: 1-quantum_wire_revisited.py
+       :start-after: #HIDDEN_BEGIN_nooi
+       :end-before: #HIDDEN_END_nooi
 
      we write ``*hopping`` instead of ``hopping``. The reason is as follows:
      `~kwant.builder.Builder.possible_hoppings` expects the hopping to
@@ -395,8 +417,9 @@ The result of the example should be identical to the previous one.
      using a tuple of sites. Hence it is worth noting
      a subtle detail in
 
-     .. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
-         :lines: 27
+     .. literalinclude:: 1-quantum_wire_revisited.py
+         :start-after: #HIDDEN_BEGIN_vvjt
+         :end-before: #HIDDEN_END_vvjt
 
      Note that ``(lat(x, y) for x in range(L) for y in range(W))`` is not
      a tuple, but a generator.
diff --git a/doc/source/tutorial/tutorial2.rst b/doc/source/tutorial/tutorial2.rst
index e0046e16de36b6eb9ce9b5ff1ddb45cca2cc5e3c..f78e0d576efcb90f6b80357f7bbe3ca33517199c 100644
--- a/doc/source/tutorial/tutorial2.rst
+++ b/doc/source/tutorial/tutorial2.rst
@@ -37,21 +37,24 @@ In order to deal with matrices in python, kwant uses the `numpy package
 <numpy.scipy.org>`_. In order to use matrices in our program, we thus also
 have to import that package:
 
-.. literalinclude:: ../../../tutorial/2-spin_orbit.py
-    :lines: 19
+.. literalinclude:: 2-spin_orbit.py
+    :start-after: #HIDDEN_BEGIN_xumz
+    :end-before: #HIDDEN_END_xumz
 
 For convenience, we define the Pauli-matrices first (with `sigma_0` the
 unit matrix):
 
-.. literalinclude:: ../../../tutorial/2-spin_orbit.py
-    :lines: 22-25
+.. literalinclude:: 2-spin_orbit.py
+    :start-after: #HIDDEN_BEGIN_hwbt
+    :end-before: #HIDDEN_END_hwbt
 
 Previously, we used numbers as the values of our matrix elements.
 However, `~kwant.builder.Builder` also accepts matrices as values, and
 we can simply write:
 
-.. literalinclude:: ../../../tutorial/2-spin_orbit.py
-    :lines: 36-43
+.. literalinclude:: 2-spin_orbit.py
+    :start-after: #HIDDEN_BEGIN_uxrm
+    :end-before: #HIDDEN_END_uxrm
 
 Note that the Zeeman energy adds to the onsite term, whereas the Rashba
 spin-orbit term adds to the hoppings (due to the derivative operator).
@@ -78,8 +81,9 @@ when specifying `(1, 0)` it is not necessary to specify `(-1, 0)`),
 
 The leads also allow for a matrix structure,
 
-.. literalinclude:: ../../../tutorial/2-spin_orbit.py
-    :lines: 50-56
+.. literalinclude:: 2-spin_orbit.py
+    :start-after: #HIDDEN_BEGIN_yliu
+    :end-before: #HIDDEN_END_yliu
 
 The remainder of the code is unchanged, and as a result we should obtain
 the following, clearly non-monotonic conductance steps:
@@ -124,8 +128,9 @@ changing the potential then implies the need to build up the system again.
 Instead, we use a python *function* to define the onsite energies. We
 define the potential profile of a quantum well as:
 
-.. literalinclude:: ../../../tutorial/2-quantum_well.py
-    :lines: 16-18, 21-24, 27-33
+.. literalinclude:: 2-quantum_well.py
+    :start-after: #HIDDEN_BEGIN_ehso
+    :end-before: #HIDDEN_END_ehso
 
 This function takes one argument which is of type
 `~kwant.builder.Site`, from which you can get the realspace
@@ -140,8 +145,9 @@ the transmission as a function of well depth.
 kwant now allows us to pass a function as a value to
 `~kwant.builder.Builder`:
 
-.. literalinclude:: ../../../tutorial/2-quantum_well.py
-    :lines: 35-40
+.. literalinclude:: 2-quantum_well.py
+    :start-after: #HIDDEN_BEGIN_coid
+    :end-before: #HIDDEN_END_coid
 
 For each lattice point, the corresponding site is then passed to the
 function `onsite()`. Note that we had to define `onsite()`, as it is
@@ -155,8 +161,9 @@ of the lead -- this should be kept in mind.
 
 Finally, we compute the transmission probability:
 
-.. literalinclude:: ../../../tutorial/2-quantum_well.py
-    :lines: 67, 70-82
+.. literalinclude:: 2-quantum_well.py
+    :start-after: #HIDDEN_BEGIN_sqvr
+    :end-before: #HIDDEN_END_sqvr
 
 Since we change the value of the global variable `pot` to vary the
 well depth, python requires us to write ``global pot`` to `enable
@@ -185,10 +192,9 @@ oscillatory transmission behavior through resonances in the quantum well.
   - Functions can also be used for hoppings. In this case, they take
     two `~kwant.builder.Site`'s as arguments.
 
-  - In tutorial/2-quantum_well.py, line 16
+  - In tutorial/2-quantum_well.py, the line ::
 
-    .. literalinclude:: ../../../tutorial/2-quantum_well.py
-        :lines: 16
+        pot = 0
 
     is not really necessary. If this line was left out, the
     global variable `pot` would in fact be created by the
@@ -290,8 +296,9 @@ First, define a boolean function defining the desired shape, i.e. a function
 that returns ``True`` whenever a point is inside the shape, and
 ``False`` otherwise:
 
-.. literalinclude:: ../../../tutorial/2-ab_ring.py
-    :lines: 21, 24-27, 31-34
+.. literalinclude:: 2-ab_ring.py
+    :start-after: #HIDDEN_BEGIN_eusz
+    :end-before: #HIDDEN_END_eusz
 
 Note that this function takes a realspace position as argument (not a
 `~kwant.builder.Site`).
@@ -300,8 +307,9 @@ We can now simply add all of the lattice points inside this shape at
 once, using the function `~kwant.lattice.Square.shape`
 provided by the lattice:
 
-.. literalinclude:: ../../../tutorial/2-ab_ring.py
-    :lines: 37-39
+.. literalinclude:: 2-ab_ring.py
+    :start-after: #HIDDEN_BEGIN_lcak
+    :end-before: #HIDDEN_END_lcak
 
 Here, ``lat.shape()`` takes as a second parameter a (realspace) point
 that is inside the desired shape. The hoppings can still be added
@@ -314,8 +322,9 @@ along the branch cut in the lower arm of the ring. For this we select
 all hoppings in x-direction that are of the form `(lat(1, j), lat(0, j))`
 with ``j<0``:
 
-.. literalinclude:: ../../../tutorial/2-ab_ring.py
-    :lines: 47-59
+.. literalinclude:: 2-ab_ring.py
+    :start-after: #HIDDEN_BEGIN_lvkt
+    :end-before: #HIDDEN_END_lvkt
 
 Here, `crosses_branchcut` is a boolean function that returns ``True`` for
 the desired hoppings. We then use again a generator (this time with
@@ -328,16 +337,18 @@ by the global variable `phi`.
 
 For the leads, we can also use the ``lat.shape()``-functionality:
 
-.. literalinclude:: ../../../tutorial/2-ab_ring.py
-    :lines: 63-72
+.. literalinclude:: 2-ab_ring.py
+    :start-after: #HIDDEN_BEGIN_qwgr
+    :end-before: #HIDDEN_END_qwgr
 
 Here, the shape must cover *at least* one unit cell of the lead
 (it does not hurt if it covers more unit cells).
 
 Attaching the leads is done as before:
 
-.. literalinclude:: ../../../tutorial/2-ab_ring.py
-    :lines: 79-80
+.. literalinclude:: 2-ab_ring.py
+    :start-after: #HIDDEN_BEGIN_skbz
+    :end-before: #HIDDEN_END_skbz
 
 In fact, attaching leads seems not so simple any more for the current
 structure with a scattering region very much different from the lead
diff --git a/doc/source/tutorial/tutorial3.rst b/doc/source/tutorial/tutorial3.rst
index 5907ff947a7541830f8d2f2454bd8b45abdb5a75..ac334262cba648227a7343ee90fb499500c2737d 100644
--- a/doc/source/tutorial/tutorial3.rst
+++ b/doc/source/tutorial/tutorial3.rst
@@ -15,8 +15,9 @@ tight-binding wire.
 Computing band structures in kwant is easy. Just define a lead in the
 usual way:
 
-.. literalinclude:: ../../../tutorial/3-band_structure.py
-    :lines: 17-34
+.. literalinclude:: 3-band_structure.py
+    :start-after: #HIDDEN_BEGIN_zxip
+    :end-before: #HIDDEN_END_zxip
 
 "Usual way" means defining a translational symmetry vector, as well
 as one unit cell of the lead, and the hoppings to neighboring
@@ -30,8 +31,9 @@ that allows to compute the eigenenergies of the translational invariant system
 for a given momentum `k`. Computing these eigenenergies for different momenta
 `k` then yields the bandstructure:
 
-.. literalinclude:: ../../../tutorial/3-band_structure.py
-    :lines: 37 - 55
+.. literalinclude:: 3-band_structure.py
+    :start-after: #HIDDEN_BEGIN_pejz
+    :end-before: #HIDDEN_END_pejz
 
 This gives the result:
 
@@ -75,14 +77,16 @@ circular quantum dot as a function of magnetic field
 To compute the eigenenergies, we will make use of the linear algebra
 functionality of `scipy <www.scipy.org>`_:
 
-.. literalinclude:: ../../../tutorial/3-closed_system.py
-    :lines: 16
+.. literalinclude:: 3-closed_system.py
+    :start-after: #HIDDEN_BEGIN_tibv
+    :end-before: #HIDDEN_END_tibv
 
 We set up the system using the `shape`-function as in
 :ref:`tutorial-abring`, but do not add any leads:
 
-.. literalinclude:: ../../../tutorial/3-closed_system.py
-    :lines: 26-48
+.. literalinclude:: 3-closed_system.py
+    :start-after: #HIDDEN_BEGIN_qlyd
+    :end-before: #HIDDEN_END_qlyd
 
 We add the magnetic field using a function and a global variable as we
 did in the two previous tutorial. (Here, the gauge is chosen such that
@@ -92,8 +96,9 @@ The spectrum can be obtained by diagonalizing the Hamiltonian of the
 system, which in turn can be obtained from the finalized
 system using `~kwant.system.System.hamiltonian_submatrix`:
 
-.. literalinclude:: ../../../tutorial/3-closed_system.py
-    :lines: 51, 53, 60-76
+.. literalinclude:: 3-closed_system.py
+    :start-after: #HIDDEN_BEGIN_yvri
+    :end-before: #HIDDEN_END_yvri
 
 In this toy model we use dense matrices and dense matrix algebra since
 the system is very small. (In a real application one would probably
diff --git a/doc/source/tutorial/tutorial4.rst b/doc/source/tutorial/tutorial4.rst
index a93e4ae4bf7ed9bdaeded80b7a6d011f11bf0bea..4d9eb5f4402e377fba10f9616168b266e700bbb7 100644
--- a/doc/source/tutorial/tutorial4.rst
+++ b/doc/source/tutorial/tutorial4.rst
@@ -14,8 +14,9 @@ We begin by defining the honeycomb lattice of graphene. This is
 in principle already done in `kwant.lattice.Honeycomb`, but we do it
 explicitly here to show how to define a new lattice:
 
-.. literalinclude:: ../../../tutorial/4-graphene.py
-    :lines: 24-27
+.. literalinclude:: 4-graphene.py
+    :start-after: #HIDDEN_BEGIN_hnla
+    :end-before: #HIDDEN_END_hnla
 
 The first argument to the `~kwant.lattice.make_lattice` function is the list of
 primitive vectors of the lattice; the second one is the coordinates of basis
@@ -26,8 +27,9 @@ itself forms a regular lattice of the same type as well, and those
 In the next step we define the shape of the scattering region (circle again)
 and add all lattice points using the ``shape()``-functionality:
 
-.. literalinclude:: ../../../tutorial/4-graphene.py
-    :lines: 29-30, 33-38, 40-45
+.. literalinclude:: 4-graphene.py
+    :start-after: #HIDDEN_BEGIN_shzy
+    :end-before: #HIDDEN_END_shzy
 
 As you can see, this works exactly the same for any kind of lattice.
 We add the onsite energies using a function describing the p-n junction;
@@ -40,8 +42,9 @@ As a next step we add the hoppings, making use of
 lattice (instead of `kwant.lattice.Honeycomb`), we have to define
 the hoppings ourselves:
 
-.. literalinclude:: ../../../tutorial/4-graphene.py
-    :lines: 49
+.. literalinclude:: 4-graphene.py
+    :start-after: #HIDDEN_BEGIN_hsmc
+    :end-before: #HIDDEN_END_hsmc
 
 The nearest-neighbor model for graphene contains only
 hoppings between different basis atoms. For these type of
@@ -57,24 +60,27 @@ respect to the two primitive vectors ``[(1, 0), (sin_30, cos_30)]``.
 
 Adding the hoppings however still works the same way:
 
-.. literalinclude:: ../../../tutorial/4-graphene.py
-    :lines: 50-51
+.. literalinclude:: 4-graphene.py
+    :start-after: #HIDDEN_BEGIN_bfwb
+    :end-before: #HIDDEN_END_bfwb
 
 Modifying the scattering region is also possible as before. Let's
 do something crazy, and remove an atom in sublattice A
 (which removes also the hoppings from/to this site) as well
 as add an additional link:
 
-.. literalinclude:: ../../../tutorial/4-graphene.py
-    :lines: 54-55
+.. literalinclude:: 4-graphene.py
+    :start-after: #HIDDEN_BEGIN_efut
+    :end-before: #HIDDEN_END_efut
 
 Note again that the conversion from a tuple `(i,j)` to site
 is done by the sublattices `a` and `b`.
 
 The leads are defined as before:
 
-.. literalinclude:: ../../../tutorial/4-graphene.py
-    :lines: 58-82
+.. literalinclude:: 4-graphene.py
+    :start-after: #HIDDEN_BEGIN_aakh
+    :end-before: #HIDDEN_END_aakh
 
 Note that the translational vectors ``graphene.vec((-1, 0))`` and
 ``graphene.vec((0, 1))`` are *not* orthogonal any more as they would have been
@@ -85,14 +91,16 @@ Later, we will compute some eigenvalues of the closed scattering region without
 leads. This is why we postpone attaching the leads to the system. Instead,
 we return the scattering region and the leads separately.
 
-.. literalinclude:: ../../../tutorial/4-graphene.py
-    :lines: 84
+.. literalinclude:: 4-graphene.py
+    :start-after: #HIDDEN_BEGIN_kmmw
+    :end-before: #HIDDEN_END_kmmw
 
 The computation of some eigenvalues of the closed system is done
 in the following piece of code:
 
-.. literalinclude:: ../../../tutorial/4-graphene.py
-    :lines: 87-91, 95-98
+.. literalinclude:: 4-graphene.py
+    :start-after: #HIDDEN_BEGIN_zydk
+    :end-before: #HIDDEN_END_zydk
 
 Here we use in contrast to the previous example a sparse matrix and
 the sparse linear algebra functionality of scipy (this requires
@@ -106,16 +114,19 @@ to the previous examples, and needs not be further explained here.
 Finally, in the `main()` function we make and
 plot the system:
 
-.. literalinclude:: ../../../tutorial/4-graphene.py
-    :lines: 127-130, 139
+.. literalinclude:: 4-graphene.py
+    :start-after: #HIDDEN_BEGIN_itkk
+    :end-before: #HIDDEN_END_itkk
 
 We customize the plotting: `plotter_symbols` is a dictionary with the
 sublattice objects `a` and `b` as keys, and the `~kwant.plotter.Circle` objects
 specify that the sublattice `a` should be drawn using a filled black circle,
-and `b` using a white circle with a black outline.
+and `b` using a white circle with a black outline. ::
 
-.. literalinclude:: ../../../tutorial/4-graphene.py
-    :lines: 133-136
+    plotter_symbols = {a: kwant.plotter.Circle(r=0.3),
+                       b: kwant.plotter.Circle(r=0.3,
+                                               fcol=kwant.plotter.white,
+                                               lcol=kwant.plotter.black)}
 
 The radius of the circle is given in relative units: `~kwant.plotter.plot` uses
 a typical length scale as a reference length. By default, the typical length
@@ -135,8 +146,9 @@ Plotting the closed system gives this result:
 
 Computing the eigenvalues of largest magnitude,
 
-.. literalinclude:: ../../../tutorial/4-graphene.py
-    :lines: 142
+.. literalinclude:: 4-graphene.py
+    :start-after: #HIDDEN_BEGIN_jmbi
+    :end-before: #HIDDEN_END_jmbi
 
 should yield two eigenvalues similar to `[ 3.07869311 +1.02714523e-17j,
 -3.06233144 -6.68085759e-18j]` (round-off might change the imaginary part which
diff --git a/doc/source/tutorial/tutorial5.rst b/doc/source/tutorial/tutorial5.rst
index 3617286f8f6977cfb3b0a7377b674feac9ad2b84..a7dbbbfbd98bec7201e7b38e02820c4dd18be7a8 100644
--- a/doc/source/tutorial/tutorial5.rst
+++ b/doc/source/tutorial/tutorial5.rst
@@ -28,8 +28,9 @@ We begin by computing the band structure of a superconducting wire.
 The most natural way to implement the BdG Hamiltonian is by using a
 2x2 matrix structure for all Hamiltonian matrix elements:
 
-.. literalinclude:: ../../../tutorial/5-superconductor_band_structure.py
-    :lines: 21-42
+.. literalinclude:: 5-superconductor_band_structure.py
+    :start-after: #HIDDEN_BEGIN_nbvn
+    :end-before: #HIDDEN_END_nbvn
 
 As you see, the example is syntactically equivalent to our
 :ref:`spin example <tutorial_spinorbit>`, the only difference
@@ -80,8 +81,9 @@ a superconductor on the right, and a tunnel barrier inbetween:
 As already mentioned above, we begin by introducing two different
 square lattices representing electron and hole degrees of freedom:
 
-.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
-    :lines: 16-17,15,20-21
+.. literalinclude:: 5-superconductor_transport.py
+    :start-after: #HIDDEN_BEGIN_zuuw
+    :end-before: #HIDDEN_END_zuuw
 
 Any diagonal entry (kinetic energy, potentials, ...) in the BdG
 Hamiltonian then corresponds to on-site energies or hoppings within
@@ -89,8 +91,9 @@ the *same* lattice, whereas any off-diagonal entry (essentially, the
 superconducting order parameter :math:`\Delta`) corresponds
 to a hopping between *different* lattices:
 
-.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
-    :lines: 23-44
+.. literalinclude:: 5-superconductor_transport.py
+    :start-after: #HIDDEN_BEGIN_pqmp
+    :end-before: #HIDDEN_END_pqmp
 
 Note that the tunnel barrier is added by overwriting previously set
 on-site matrix elements.
@@ -103,8 +106,9 @@ part. We use this fact to attach purely electron and hole leads
 (comprised of only electron *or* hole lattices) to the
 system:
 
-.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
-    :lines: 47-63
+.. literalinclude:: 5-superconductor_transport.py
+    :start-after: #HIDDEN_BEGIN_ttth
+    :end-before: #HIDDEN_END_ttth
 
 This separation into two different leads allows us then later to compute the
 reflection probablities between electrons and holes explicitely.
@@ -112,15 +116,17 @@ reflection probablities between electrons and holes explicitely.
 On the superconducting side, we cannot do this separation, and can
 only define a single lead coupling electrons and holes:
 
-.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
-    :lines: 68-78
+.. literalinclude:: 5-superconductor_transport.py
+    :start-after: #HIDDEN_BEGIN_mhiw
+    :end-before: #HIDDEN_END_mhiw
 
 We now have on the left side two leads that are sitting in the same
 spatial position, but in different lattice spaces. This ensures that
 we can still attach all leads as before:
 
-.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
-    :lines: 81-85
+.. literalinclude:: 5-superconductor_transport.py
+    :start-after: #HIDDEN_BEGIN_ozsr
+    :end-before: #HIDDEN_END_ozsr
 
 When computing the conductance, we can now extract reflection from
 electrons to electrons as ``smatrix.transmission(0, 0)`` (Don't get
@@ -128,8 +134,9 @@ confused by the fact that it says ``transmission`` -- transmission
 into the same lead is reflection), and reflection from electrons to holes
 as ``smatrix.transmission(1, 0)``, by virtue of our electron and hole leads:
 
-.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
-    :lines: 88-96
+.. literalinclude:: 5-superconductor_transport.py
+    :start-after: #HIDDEN_BEGIN_jbjt
+    :end-before: #HIDDEN_END_jbjt
 
 Note that ``smatrix.submatrix(0,0)`` returns the block concerning reflection
 within (electron) lead 0, and from its size we can extract the number of modes
diff --git a/setup.py b/setup.py
index ef029c03c7a76ec97323682e63f3d878bfe7e7b6..e1a5f2dcb6ec4c2cd539a6a5c0d3e7766562a89a 100755
--- a/setup.py
+++ b/setup.py
@@ -5,14 +5,19 @@ STATIC_VERSION_FILE = 'kwant/_static_version.py'
 REQUIRED_CYTHON_VERSION = (0, 17, 1)
 MUMPS_DEBIAN_PACKAGE = 'libmumps-scotch-dev'
 NO_CYTHON_OPTION = '--no-cython'
+TUT_DIR = 'tutorial'
+TUT_GLOB = 'doc/source/tutorial/*.py'
+TUT_HIDDEN_PREFIX = '#HIDDEN'
 
 import sys
 import os
+import glob
 import subprocess
 import ConfigParser
-from distutils.core import setup
+from distutils.core import setup, Command
 from distutils.extension import Extension
 from distutils.errors import DistutilsError, CCompilerError
+from distutils.command.build import build as distutils_build
 import numpy
 
 try:
@@ -54,6 +59,36 @@ Build configuration was:
         print build_summary
 
 
+class build_tut(Command):
+    description = "build the tutorial scripts"
+    user_options = []
+
+    def initialize_options(self):
+        pass
+
+    def finalize_options(self):
+        pass
+
+    def run(self):
+        if not os.path.exists(TUT_DIR):
+            os.mkdir(TUT_DIR)
+        for in_fname in glob.glob(TUT_GLOB):
+            out_fname = os.path.join(TUT_DIR, os.path.basename(in_fname))
+            with open(in_fname) as in_file, open(out_fname, 'w') as out_file:
+                for line in in_file:
+                    if not line.startswith(TUT_HIDDEN_PREFIX):
+                        out_file.write(line)
+
+
+# Our version of the "build" command also makes sure the tutorial is made.
+# Even though the tutorial is not necessary for installation, and "build" is
+# supposed to make everything needed to install, this is a robust way to ensure
+# that the tutorial is present.
+class kwant_build(distutils_build):
+    sub_commands = [('build_tut', None)] + distutils_build.sub_commands
+    pass
+
+
 # This is an exact copy of the function from kwant/version.py.  We can't import
 # it here (because kwant is not yet built when this scipt is run), so we just
 # include a copy.
@@ -279,7 +314,9 @@ def main():
           license="not to be distributed",
           packages=["kwant", "kwant.graph", "kwant.linalg", "kwant.physics",
                     "kwant.solvers"],
-          cmdclass={'build_ext': kwant_build_ext},
+          cmdclass={'build': kwant_build,
+                    'build_ext': kwant_build_ext,
+                    'build_tut': build_tut},
           ext_modules=ext_modules(extensions()),
           include_dirs=[numpy.get_include()])