diff --git a/doc/Makefile b/doc/Makefile
index da1b700b126f26a5036bb6655e62950aff3cda13..7be69e65d7b196e470d1c4bad1d9aadf8afc83b9 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -19,15 +19,15 @@ IMAGESOURCES    = $(shell find source -name "*.svg")
 GENERATEDPDF    = $(patsubst %.svg, %.pdf, $(IMAGESOURCES))
 
 # Tutorial images.
-TUTORIAL1A_IMAGES = source/images/tutorial1a_result.png source/images/tutorial1a_result.pdf source/images/tutorial1a_sys.png source/images/tutorial1a_sys.pdf
-TUTORIAL2A_IMAGES = source/images/tutorial2a_result.png source/images/tutorial2a_result.pdf
-TUTORIAL2B_IMAGES = source/images/tutorial2b_result.png source/images/tutorial2b_result.pdf
-TUTORIAL2C_IMAGES = source/images/tutorial2c_result.png source/images/tutorial2c_result.pdf source/images/tutorial2c_sys.png source/images/tutorial2c_sys.pdf source/images/tutorial2c_note1.png source/images/tutorial2c_note1.pdf source/images/tutorial2c_note2.png source/images/tutorial2c_note2.pdf
-TUTORIAL3A_IMAGES = source/images/tutorial3a_result.png source/images/tutorial3a_result.pdf
-TUTORIAL3B_IMAGES = source/images/tutorial3b_result.png source/images/tutorial3b_result.pdf source/images/tutorial3b_sys.png source/images/tutorial3b_sys.pdf
-TUTORIAL4_IMAGES = source/images/tutorial4_result.png source/images/tutorial4_result.pdf source/images/tutorial4_sys1.png source/images/tutorial4_sys1.pdf source/images/tutorial4_sys2.png source/images/tutorial4_sys2.pdf source/images/tutorial4_bs.png source/images/tutorial4_bs.pdf
-TUTORIAL5A_IMAGES = source/images/tutorial5a_result.png source/images/tutorial5a_result.pdf
-TUTORIAL5B_IMAGES = source/images/tutorial5b_result.png source/images/tutorial5b_result.pdf
+TUTORIAL1A_IMAGES = source/images/1-quantum_wire_result.png source/images/1-quantum_wire_result.pdf source/images/1-quantum_wire_sys.png source/images/1-quantum_wire_sys.pdf
+TUTORIAL2A_IMAGES = source/images/2-spin_orbit_result.png source/images/2-spin_orbit_result.pdf
+TUTORIAL2B_IMAGES = source/images/2-quantum_well_result.png source/images/2-quantum_well_result.pdf
+TUTORIAL2C_IMAGES = source/images/2-ab_ring_result.png source/images/2-ab_ring_result.pdf source/images/2-ab_ring_sys.png source/images/2-ab_ring_sys.pdf source/images/2-ab_ring_note1.png source/images/2-ab_ring_note1.pdf source/images/2-ab_ring_note2.png source/images/2-ab_ring_note2.pdf
+TUTORIAL3A_IMAGES = source/images/3-band_structure_result.png source/images/3-band_structure_result.pdf
+TUTORIAL3B_IMAGES = source/images/3-closed_system_result.png source/images/3-closed_system_result.pdf source/images/3-closed_system_sys.png source/images/3-closed_system_sys.pdf
+TUTORIAL4_IMAGES = source/images/4-graphene_result.png source/images/4-graphene_result.pdf source/images/4-graphene_sys1.png source/images/4-graphene_sys1.pdf source/images/4-graphene_sys2.png source/images/4-graphene_sys2.pdf source/images/4-graphene_bs.png source/images/4-graphene_bs.pdf
+TUTORIAL5A_IMAGES = source/images/5-superconductor_band_structure_result.png source/images/5-superconductor_band_structure_result.pdf
+TUTORIAL5B_IMAGES = source/images/5-superconductor_transport_result.png source/images/5-superconductor_transport_result.pdf
 ALL_TUTORIAL_IMAGES = $(TUTORIAL1A_IMAGES)  $(TUTORIAL2A_IMAGES) $(TUTORIAL2B_IMAGES) $(TUTORIAL2C_IMAGES) $(TUTORIAL3A_IMAGES) $(TUTORIAL3B_IMAGES) $(TUTORIAL4_IMAGES) $(TUTORIAL5A_IMAGES) $(TUTORIAL5B_IMAGES)
 
 .PHONY: help clean realclean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest
@@ -115,56 +115,56 @@ doctest:
 
 # Generation of tutorial images.  This requires some serious make trickery, see
 # http://article.gmane.org/gmane.comp.gnu.make.general/5806
-$(TUTORIAL1A_IMAGES): source/images/.tutorial1a_flag
+$(TUTORIAL1A_IMAGES): source/images/.1-quantum_wire_flag
 	@:
-source/images/.tutorial1a_flag: source/images/tutorial1a.py
-	cd source/images/ && python tutorial1a.py
-	touch source/images/.tutorial1a_flag
+source/images/.1-quantum_wire_flag: source/images/1-quantum_wire.py
+	cd source/images/ && python 1-quantum_wire.py
+	touch source/images/.1-quantum_wire_flag
 
-$(TUTORIAL2A_IMAGES): source/images/.tutorial2a_flag
+$(TUTORIAL2A_IMAGES): source/images/.2-spin_orbit_flag
 	@:
-source/images/.tutorial2a_flag: source/images/tutorial2a.py
-	cd source/images/ && python tutorial2a.py
-	touch source/images/.tutorial2a_flag
+source/images/.2-spin_orbit_flag: source/images/2-spin_orbit.py
+	cd source/images/ && python 2-spin_orbit.py
+	touch source/images/.2-spin_orbit_flag
 
-$(TUTORIAL2B_IMAGES): source/images/.tutorial2b_flag
+$(TUTORIAL2B_IMAGES): source/images/.2-quantum_well_flag
 	@:
-source/images/.tutorial2b_flag: source/images/tutorial2b.py
-	cd source/images/ && python tutorial2b.py
-	touch source/images/.tutorial2b_flag
+source/images/.2-quantum_well_flag: source/images/2-quantum_well.py
+	cd source/images/ && python 2-quantum_well.py
+	touch source/images/.2-quantum_well_flag
 
-$(TUTORIAL2C_IMAGES): source/images/.tutorial2c_flag
+$(TUTORIAL2C_IMAGES): source/images/.2-ab_ring_flag
 	@:
-source/images/.tutorial2c_flag: source/images/tutorial2c.py
-	cd source/images/ && python tutorial2c.py
-	touch source/images/.tutorial2c_flag
+source/images/.2-ab_ring_flag: source/images/2-ab_ring.py
+	cd source/images/ && python 2-ab_ring.py
+	touch source/images/.2-ab_ring_flag
 
-$(TUTORIAL3A_IMAGES): source/images/.tutorial3a_flag
+$(TUTORIAL3A_IMAGES): source/images/.3-band_structure_flag
 	@:
-source/images/.tutorial3a_flag: source/images/tutorial3a.py
-	cd source/images/ && python tutorial3a.py
-	touch source/images/.tutorial3a_flag
+source/images/.3-band_structure_flag: source/images/3-band_structure.py
+	cd source/images/ && python 3-band_structure.py
+	touch source/images/.3-band_structure_flag
 
-$(TUTORIAL3B_IMAGES): source/images/.tutorial3b_flag
+$(TUTORIAL3B_IMAGES): source/images/.3-closed_system_flag
 	@:
-source/images/.tutorial3b_flag: source/images/tutorial3b.py
-	cd source/images/ && python tutorial3b.py
-	touch source/images/.tutorial3b_flag
+source/images/.3-closed_system_flag: source/images/3-closed_system.py
+	cd source/images/ && python 3-closed_system.py
+	touch source/images/.3-closed_system_flag
 
-$(TUTORIAL4_IMAGES): source/images/.tutorial4_flag
+$(TUTORIAL4_IMAGES): source/images/.4-graphene_flag
 	@:
-source/images/.tutorial4_flag: source/images/tutorial4.py
-	cd source/images/ && python tutorial4.py
-	touch source/images/.tutorial4_flag
+source/images/.4-graphene_flag: source/images/4-graphene.py
+	cd source/images/ && python 4-graphene.py
+	touch source/images/.4-graphene_flag
 
-$(TUTORIAL5A_IMAGES): source/images/.tutorial5a_flag
+$(TUTORIAL5A_IMAGES): source/images/.5-superconductor_band_structure_flag
 	@:
-source/images/.tutorial5a_flag: source/images/tutorial5a.py
-	cd source/images/ && python tutorial5a.py
-	touch source/images/.tutorial5a_flag
+source/images/.5-superconductor_band_structure_flag: source/images/5-superconductor_band_structure.py
+	cd source/images/ && python 5-superconductor_band_structure.py
+	touch source/images/.5-superconductor_band_structure_flag
 
-$(TUTORIAL5B_IMAGES): source/images/.tutorial5b_flag
+$(TUTORIAL5B_IMAGES): source/images/.5-superconductor_transport_flag
 	@:
-source/images/.tutorial5b_flag: source/images/tutorial5b.py
-	cd source/images/ && python tutorial5b.py
-	touch source/images/.tutorial5b_flag
+source/images/.5-superconductor_transport_flag: source/images/5-superconductor_transport.py
+	cd source/images/ && python 5-superconductor_transport.py
+	touch source/images/.5-superconductor_transport_flag
diff --git a/doc/source/images/tutorial1a.py b/doc/source/images/1-quantum_wire.py
similarity index 65%
rename from doc/source/images/tutorial1a.py
rename to doc/source/images/1-quantum_wire.py
index 41c5fefbc51ec126a975fad2ed1f0e3eab4922b8..835b0a9653110f7f7b53cd1fb315a8097d4dfcb5 100644
--- a/doc/source/images/tutorial1a.py
+++ b/doc/source/images/1-quantum_wire.py
@@ -9,7 +9,6 @@
 #  - Using the simple sparse solver for computing Landauer conductance
 
 import kwant
-
 import latex, html
 
 # First, define the tight-binding system
@@ -17,9 +16,8 @@ import latex, html
 sys = kwant.Builder()
 
 # Here, we are only working with square lattices
-
-lat = kwant.lattice.Square()
-sys.default_site_group = lat
+a = 1
+lat = kwant.lattice.Square(a)
 
 t = 1.0
 W = 10
@@ -29,15 +27,15 @@ L = 30
 
 for i in xrange(L):
     for j in xrange(W):
-        sys[(i, j)] = 4 * t
+        sys[lat(i, j)] = 4 * t
 
         # hoppig in y-direction
-        if j > 0 :
-            sys[(i, j), (i, j-1)] = - t
+        if j > 0:
+            sys[lat(i, j), lat(i, j - 1)] = - t
 
         #hopping in x-direction
         if i > 0:
-            sys[(i, j), (i-1, j)] = -t
+            sys[lat(i, j), lat(i - 1, j)] = -t
 
 # Then, define the leads:
 
@@ -47,29 +45,27 @@ for i in xrange(L):
 # realspace vector)
 sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
 lead0 = kwant.Builder(sym_lead0)
-lead0.default_site_group = lat
 
 for j in xrange(W):
-    lead0[(0, j)] = 4 * t
+    lead0[lat(0, j)] = 4 * t
 
     if j > 0:
-        lead0[(0, j), (0, j-1)] = - t
+        lead0[lat(0, j), lat(0, j - 1)] = - t
 
-    lead0[(1, j), (0, j)] = - t
+    lead0[lat(1, j), lat(0, j)] = - t
 
 # Then the lead to the right
 
 sym_lead1 = kwant.TranslationalSymmetry([lat.vec((1, 0))])
 lead1 = kwant.Builder(sym_lead1)
-lead1.default_site_group = lat
 
 for j in xrange(W):
-    lead1[(0, j)] = 4 * t
+    lead1[lat(0, j)] = 4 * t
 
     if j > 0:
-        lead1[(0, j), (0, j-1)] = - t
+        lead1[lat(0, j), lat(0, j - 1)] = - t
 
-    lead1[(1, j), (0, j)] = - t
+    lead1[lat(1, j), lat(0, j)] = - t
 
 # Then attach the leads to the system
 
@@ -78,12 +74,12 @@ sys.attach_lead(lead1)
 
 # Plot it, to make sure it's OK
 
-kwant.plot(sys, "tutorial1a_sys.pdf", width=latex.figwidth_pt)
-kwant.plot(sys, "tutorial1a_sys.png", width=html.figwidth_px)
+kwant.plot(sys, "1-quantum_wire_sys.pdf", width=latex.figwidth_pt)
+kwant.plot(sys, "1-quantum_wire_sys.png", width=html.figwidth_px)
 
-# finalize the system
+# Finalize the system
 
-fsys = sys.finalized()
+sys = sys.finalized()
 
 # Now that we have the system, we can compute conductance
 
@@ -93,7 +89,7 @@ for ie in xrange(100):
     energy = ie * 0.01
 
     # compute the scattering matrix at energy energy
-    smatrix = kwant.solve(fsys, energy)
+    smatrix = kwant.solve(sys, energy)
 
     # compute the transmission probability from lead 0 to
     # lead 1
@@ -102,19 +98,19 @@ for ie in xrange(100):
 
 # Use matplotlib to write output
 # We should see conductance steps
-import pylab
+from matplotlib import pyplot
 
-pylab.plot(energies, data)
-pylab.xlabel("energy [in units of t]",
+fig = pyplot.figure()
+pyplot.plot(energies, data)
+pyplot.xlabel("energy [in units of t]",
                  fontsize=latex.mpl_label_size)
-pylab.ylabel("conductance [in units of e^2/h]",
+pyplot.ylabel("conductance [in units of e^2/h]",
                  fontsize=latex.mpl_label_size)
-fig = pylab.gcf()
-pylab.setp(fig.get_axes()[0].get_xticklabels(),
+pyplot.setp(fig.get_axes()[0].get_xticklabels(),
            fontsize=latex.mpl_tick_size)
-pylab.setp(fig.get_axes()[0].get_yticklabels(),
+pyplot.setp(fig.get_axes()[0].get_yticklabels(),
            fontsize=latex.mpl_tick_size)
 fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
 fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-fig.savefig("tutorial1a_result.pdf")
-fig.savefig("tutorial1a_result.png", dpi=(html.figwidth_px/latex.mpl_width_in))
+fig.savefig("1-quantum_wire_result.pdf")
+fig.savefig("1-quantum_wire_result.png", dpi=(html.figwidth_px/latex.mpl_width_in))
diff --git a/doc/source/images/tutorial2c.py b/doc/source/images/2-ab_ring.py
similarity index 77%
rename from doc/source/images/tutorial2c.py
rename to doc/source/images/2-ab_ring.py
index e876d9f343c818d86dc404bad76481a79497ef1b..ad6ba6e14ef7feebfa220eedf98cb5e7d4beff13 100644
--- a/doc/source/images/tutorial2c.py
+++ b/doc/source/images/2-ab_ring.py
@@ -15,7 +15,7 @@ from math import pi
 import kwant
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
 import latex, html
 
@@ -31,11 +31,11 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
     # Now, we aim for a more complex shape, namely a ring (or annulus)
     def ring(pos):
         (x, y) = pos
-        rsq = x**2 + y**2
-        return ( r1**2 < rsq < r2**2)
+        rsq = x ** 2 + y ** 2
+        return (r1 ** 2 < rsq < r2 ** 2)
 
     # and add the corresponding lattice points using the `shape`-function
-    sys[lat.shape(ring, (0, 11))] = 4 * t
+    sys[lat.shape(ring, (0, r1 + 1))] = 4 * t
     for hopping in lat.nearest:
         sys[sys.possible_hoppings(*hopping)] = - t
 
@@ -53,10 +53,10 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
 
         # possible_hoppings with the argument (1, 0) below
         # returns hoppings ordered as ((i+1, j), (i, j))
-        return iy0 < 0 and ix0 == 1 # ix1 == 0 then implied
+        return iy0 < 0 and ix0 == 1  # ix1 == 0 then implied
 
     # Modify only those hopings in x-direction that cross the branch cut
-    sys[(hop for hop in sys.possible_hoppings((1,0), lat, lat)
+    sys[(hop for hop in sys.possible_hoppings((1, 0), lat, lat)
          if crosses_branchcut(hop))] = fluxphase
 
     #### Define the leads. ####
@@ -66,17 +66,17 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
 
     def lead_shape(pos):
         (x, y) = pos
-        return (-1 < x < 1) and ( -W/2 < y < W/2  )
+        return (-1 < x < 1) and (-W / 2 < y < W / 2)
 
     lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
     for hopping in lat.nearest:
         lead0[lead0.possible_hoppings(*hopping)] = - t
 
     # Then the lead to the right
-    # (again, obtained using reverse()
+    # [again, obtained using reversed()]
     lead1 = lead0.reversed()
 
-    #### Attach the leads and return the finalized system. ####
+    #### Attach the leads and return the system. ####
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
 
@@ -131,33 +131,33 @@ def make_system_note2(a=1, t=1.0, W=10, r1=10, r2=20):
     return sys
 
 
-def plot_conductance(fsys, energy, fluxes):
+def plot_conductance(sys, energy, fluxes):
     # compute conductance
     # global variable phi controls the flux
     global phi
 
-    normalized_fluxes = [flux/(2 * pi) for flux in fluxes]
+    normalized_fluxes = [flux / (2 * pi) for flux in fluxes]
     data = []
     for flux in fluxes:
         phi = flux
 
-        smatrix = kwant.solve(fsys, energy)
+        smatrix = kwant.solve(sys, energy)
         data.append(smatrix.transmission(1, 0))
 
-    pylab.plot(normalized_fluxes, data)
-    pylab.xlabel("flux [in units of the flux quantum]",
+    fig = pyplot.figure()
+    pyplot.plot(normalized_fluxes, data)
+    pyplot.xlabel("flux [in units of the flux quantum]",
                  fontsize=latex.mpl_label_size)
-    pylab.ylabel("conductance [in units of e^2/h]",
+    pyplot.ylabel("conductance [in units of e^2/h]",
                  fontsize=latex.mpl_label_size)
-    fig = pylab.gcf()
-    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
                fontsize=latex.mpl_tick_size)
-    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
                fontsize=latex.mpl_tick_size)
     fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
     fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-    fig.savefig("tutorial2c_result.pdf")
-    fig.savefig("tutorial2c_result.png",
+    fig.savefig("2-ab_ring_result.pdf")
+    fig.savefig("2-ab_ring_result.png",
                 dpi=(html.figwidth_px/latex.mpl_width_in))
 
 
@@ -165,23 +165,23 @@ def main():
     sys = make_system()
 
     # Check that the system looks as intended.
-    kwant.plot(sys, "tutorial2c_sys.pdf", width=latex.figwidth_pt)
-    kwant.plot(sys, "tutorial2c_sys.png", width=html.figwidth_px)
+    kwant.plot(sys, "2-ab_ring_sys.pdf", width=latex.figwidth_pt)
+    kwant.plot(sys, "2-ab_ring_sys.png", width=html.figwidth_px)
 
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
     # We should see a conductance that is periodic with the flux quantum
-    plot_conductance(fsys, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
+    plot_conductance(sys, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
                                                 for i in xrange(100)])
 
     # Finally, some plots needed for the notes
     sys = make_system_note1()
-    kwant.plot(sys, "tutorial2c_note1.pdf", width=latex.figwidth_small_pt)
-    kwant.plot(sys, "tutorial2c_note1.png", width=html.figwidth_small_px)
+    kwant.plot(sys, "2-ab_ring_note1.pdf", width=latex.figwidth_small_pt)
+    kwant.plot(sys, "2-ab_ring_note1.png", width=html.figwidth_small_px)
     sys = make_system_note2()
-    kwant.plot(sys, "tutorial2c_note2.pdf", width=latex.figwidth_small_pt)
-    kwant.plot(sys, "tutorial2c_note2.png", width=html.figwidth_small_px)
+    kwant.plot(sys, "2-ab_ring_note2.pdf", width=latex.figwidth_small_pt)
+    kwant.plot(sys, "2-ab_ring_note2.png", width=html.figwidth_small_px)
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/doc/source/images/tutorial2c_sketch.svg b/doc/source/images/2-ab_ring_sketch.svg
similarity index 100%
rename from doc/source/images/tutorial2c_sketch.svg
rename to doc/source/images/2-ab_ring_sketch.svg
diff --git a/doc/source/images/tutorial2c_sketch2.svg b/doc/source/images/2-ab_ring_sketch2.svg
similarity index 100%
rename from doc/source/images/tutorial2c_sketch2.svg
rename to doc/source/images/2-ab_ring_sketch2.svg
diff --git a/doc/source/images/tutorial2b.py b/doc/source/images/2-quantum_well.py
similarity index 79%
rename from doc/source/images/tutorial2b.py
rename to doc/source/images/2-quantum_well.py
index 99861e3b7e2647ab22101774fe4d1cd8edb45f19..1b90a7e1af8d48de75db8425e56d1f4acc980429 100644
--- a/doc/source/images/tutorial2b.py
+++ b/doc/source/images/2-quantum_well.py
@@ -9,7 +9,7 @@
 import kwant
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
 import latex, html
 
@@ -17,13 +17,13 @@ import latex, html
 # make_system()
 pot = 0
 
+
 def make_system(a=1, t=1.0, W=10, L=30, L_well=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.
     lat = kwant.lattice.Square(a)
 
     sys = kwant.Builder()
-    sys.default_site_group = lat
 
     #### Define the scattering region. ####
     # Potential profile
@@ -38,7 +38,7 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
     def onsite(site):
         return 4 * t + potential(site)
 
-    sys[((x, y) for x in range(L) for y in range(W))] = onsite
+    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
 
@@ -48,9 +48,8 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
     # realspace vector)
     sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead0 = kwant.Builder(sym_lead0)
-    lead0.default_site_group = lat
 
-    lead0[((0, j) for j in xrange(W))] = 4 * t
+    lead0[(lat(0, j) for j in xrange(W))] = 4 * t
     for hopping in lat.nearest:
         lead0[lead0.possible_hoppings(*hopping)] = - t
 
@@ -64,7 +63,8 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
 
     return sys
 
-def plot_conductance(fsys, energy, welldepths):
+
+def plot_conductance(sys, energy, welldepths):
     # We specify that we want to not only read, but also write to a
     # global variable.
     global pot
@@ -75,23 +75,23 @@ def plot_conductance(fsys, energy, welldepths):
         # Set the global variable that defines the potential well depth
         pot = -welldepth
 
-        smatrix = kwant.solve(fsys, energy)
+        smatrix = kwant.solve(sys, energy)
         data.append(smatrix.transmission(1, 0))
 
-    pylab.plot(welldepths, data)
-    pylab.xlabel("well depth [in units of t]",
+    fig = pyplot.figure()
+    pyplot.plot(welldepths, data)
+    pyplot.xlabel("well depth [in units of t]",
                  fontsize=latex.mpl_label_size)
-    pylab.ylabel("conductance [in units of e^2/h]",
+    pyplot.ylabel("conductance [in units of e^2/h]",
                  fontsize=latex.mpl_label_size)
-    fig = pylab.gcf()
-    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
                fontsize=latex.mpl_tick_size)
-    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
                fontsize=latex.mpl_tick_size)
     fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
     fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-    fig.savefig("tutorial2b_result.pdf")
-    fig.savefig("tutorial2b_result.png",
+    fig.savefig("2-quantum_well_result.pdf")
+    fig.savefig("2-quantum_well_result.png",
                 dpi=(html.figwidth_px/latex.mpl_width_in))
 
 
@@ -99,10 +99,10 @@ def main():
     sys = make_system()
 
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
     # We should see conductance steps.
-    plot_conductance(fsys, energy=0.2,
+    plot_conductance(sys, energy=0.2,
                      welldepths=[0.01 * i for i in xrange(100)])
 
 
diff --git a/doc/source/images/tutorial2a.py b/doc/source/images/2-spin_orbit.py
similarity index 78%
rename from doc/source/images/tutorial2a.py
rename to doc/source/images/2-spin_orbit.py
index d0a1238fbae42caa3b7cb05bb1dd5a7a01e220ca..cab6df712432e338165568f2ccfc662af1462ec7 100644
--- a/doc/source/images/tutorial2a.py
+++ b/doc/source/images/2-spin_orbit.py
@@ -13,7 +13,7 @@
 import kwant
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
 # For matrix support
 import numpy
@@ -33,10 +33,9 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
     lat = kwant.lattice.Square(a)
 
     sys = kwant.Builder()
-    sys.default_site_group = lat
 
     #### Define the scattering region. ####
-    sys[((x, y) for x in range(L) for y in range(W))] = 4 * t * sigma_0 + \
+    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
     sys[sys.possible_hoppings((1, 0), lat, lat)] = - t * sigma_0 - \
@@ -49,9 +48,8 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
     # left lead
     sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead0 = kwant.Builder(sym_lead0)
-    lead0.default_site_group = lat
 
-    lead0[((0, j) for j in xrange(W))] = 4 * t * sigma_0 + e_z * sigma_z
+    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 - \
         1j * alpha * sigma_y
@@ -69,27 +67,27 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
 
     return sys
 
-def plot_conductance(fsys, energies):
+def plot_conductance(sys, energies):
     # Compute conductance
     data = []
     for energy in energies:
-        smatrix = kwant.solve(fsys, energy)
+        smatrix = kwant.solve(sys, energy)
         data.append(smatrix.transmission(1, 0))
 
-    pylab.plot(energies, data)
-    pylab.xlabel("energy [in units of t]",
+    fig = pyplot.figure()
+    pyplot.plot(energies, data)
+    pyplot.xlabel("energy [in units of t]",
                  fontsize=latex.mpl_label_size)
-    pylab.ylabel("conductance [in units of e^2/h]",
+    pyplot.ylabel("conductance [in units of e^2/h]",
                  fontsize=latex.mpl_label_size)
-    fig = pylab.gcf()
-    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
                fontsize=latex.mpl_tick_size)
-    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
                fontsize=latex.mpl_tick_size)
     fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
     fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-    fig.savefig("tutorial2a_result.pdf")
-    fig.savefig("tutorial2a_result.png",
+    fig.savefig("2-spin_orbit_result.pdf")
+    fig.savefig("2-spin_orbit_result.png",
                 dpi=(html.figwidth_px/latex.mpl_width_in))
 
 
@@ -97,10 +95,10 @@ def main():
     sys = make_system()
 
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
     # We should see non-monotonic conductance steps.
-    plot_conductance(fsys, energies=[0.01 * i - 0.3 for i in xrange(100)])
+    plot_conductance(sys, energies=[0.01 * i - 0.3 for i in xrange(100)])
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/doc/source/images/tutorial3a.py b/doc/source/images/3-band_structure.py
similarity index 64%
rename from doc/source/images/tutorial3a.py
rename to doc/source/images/3-band_structure.py
index a4068960e397e945db8e4a56e026330426cf92e2..1174523a301c08f98b2e8b235c8b7004ab638095 100644
--- a/doc/source/images/tutorial3a.py
+++ b/doc/source/images/3-band_structure.py
@@ -8,11 +8,10 @@
 
 import kwant
 
-import numpy as np
 from math import pi
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
 import latex, html
 
@@ -23,51 +22,49 @@ def make_lead(a=1, t=1.0, W=10):
 
     sym_lead = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead = kwant.Builder(sym_lead)
-    lead.default_site_group = lat
 
     # build up one unit cell of the lead, and add the hoppings
     # to the next unit cell
     for j in xrange(W):
-        lead[(0, j)] = 4 * t
+        lead[lat(0, j)] = 4 * t
 
         if j > 0:
-            lead[(0, j), (0, j-1)] = - t
+            lead[lat(0, j), lat(0, j - 1)] = - t
 
-        lead[(1, j), (0, j)] = - t
+        lead[lat(1, j), lat(0, j)] = - t
 
-    # return a finalized lead
     return lead
 
 
-def plot_bandstructure(flead, momenta):
+def plot_bandstructure(lead, momenta):
     # Use the method ``energies`` of the finalized lead to compute
     # the bandstructure
-    energy_list = [flead.energies(k) for k in momenta]
+    energy_list = [lead.energies(k) for k in momenta]
 
-    pylab.plot(momenta, energy_list)
-    pylab.xlabel("momentum [in units of (lattice constant)^-1]",
+    fig = pyplot.figure()
+    pyplot.plot(momenta, energy_list)
+    pyplot.xlabel("momentum [in units of (lattice constant)^-1]",
                  fontsize=latex.mpl_label_size)
-    pylab.ylabel("energy [in units of t]",
+    pyplot.ylabel("energy [in units of t]",
                  fontsize=latex.mpl_label_size)
-    fig = pylab.gcf()
-    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
                fontsize=latex.mpl_tick_size)
-    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
                fontsize=latex.mpl_tick_size)
     fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
     fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-    fig.savefig("tutorial3a_result.pdf")
-    fig.savefig("tutorial3a_result.png",
+    fig.savefig("3-band_structure_result.pdf")
+    fig.savefig("3-band_structure_result.png",
                 dpi=(html.figwidth_px/latex.mpl_width_in))
 
 
 def main():
-    flead = make_lead().finalized()
+    lead = make_lead().finalized()
 
     # list of momenta at which the bands should be computed
-    momenta = np.arange(-pi, pi + .01, 0.02 * pi)
+    momenta = [-pi + 0.02 * pi * i for i in xrange(101)]
 
-    plot_bandstructure(flead, momenta)
+    plot_bandstructure(lead, momenta)
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/doc/source/images/tutorial3b.py b/doc/source/images/3-closed_system.py
similarity index 74%
rename from doc/source/images/tutorial3b.py
rename to doc/source/images/3-closed_system.py
index bac129f2decfae627388395ee0f74fc0a3908805..c4dd0dc33b5b55e48591a2a6c362a62f3c873336 100644
--- a/doc/source/images/tutorial3b.py
+++ b/doc/source/images/3-closed_system.py
@@ -16,7 +16,7 @@ import kwant
 import scipy.linalg as la
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
 import latex, html
 
@@ -32,8 +32,8 @@ def make_system(a=1, t=1.0, r=10):
     # Define the quantum dot
     def circle(pos):
         (x, y) = pos
-        rsq = x**2 + y**2
-        return rsq < r**2
+        rsq = x ** 2 + y ** 2
+        return rsq < r ** 2
 
     def hopx(site1, site2):
         # The magnetic field is controlled by the global variable B
@@ -50,7 +50,7 @@ def make_system(a=1, t=1.0, r=10):
     return sys
 
 
-def plot_spectrum(fsys, Bfields):
+def plot_spectrum(sys, Bfields):
     # global variable B controls the magnetic field
     global B
 
@@ -64,39 +64,45 @@ def plot_spectrum(fsys, Bfields):
         B = Bfield
 
         # Obtain the Hamiltonian as a dense matrix
-        ham_mat = fsys.hamiltonian_submatrix()[0]
+        ham_mat = sys.hamiltonian_submatrix()[0]
 
         ev = la.eigh(ham_mat, eigvals_only=True)
 
         # we only plot the 15 lowest eigenvalues
         energies.append(ev[:15])
 
-    pylab.plot(Bfields, energies)
-    pylab.xlabel("magnetic field [some arbitrary units]",
+    fig = pyplot.figure()
+    pyplot.plot(Bfields, energies)
+    pyplot.xlabel("magnetic field [some arbitrary units]",
                  fontsize=latex.mpl_label_size)
-    pylab.ylabel("energy [in units of t]",
+    pyplot.ylabel("energy [in units of t]",
                  fontsize=latex.mpl_label_size)
-    fig = pylab.gcf()
-    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
                fontsize=latex.mpl_tick_size)
-    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
                fontsize=latex.mpl_tick_size)
     fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
     fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-    fig.savefig("tutorial3b_result.pdf")
-    fig.savefig("tutorial3b_result.png",
+    fig.savefig("3-closed_system_result.pdf")
+    fig.savefig("3-closed_system_result.png",
                 dpi=(html.figwidth_px/latex.mpl_width_in))
 
 
 def main():
     sys = make_system()
 
+    # Check that the system looks as intended.
+    kwant.plot(sys, filename="3-closed_system_sys.pdf",
+               width=latex.figwidth_pt)
+    kwant.plot(sys, filename="3-closed_system_sys.png",
+               width=html.figwidth_px)
+
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
     # We should observe energy levels that flow towards Landau
     # level energies with increasing magnetic field
-    plot_spectrum(fsys, [iB * 0.002 for iB in xrange(100)])
+    plot_spectrum(sys, [iB * 0.002 for iB in xrange(100)])
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/doc/source/images/tutorial4.py b/doc/source/images/4-graphene.py
similarity index 70%
rename from doc/source/images/tutorial4.py
rename to doc/source/images/4-graphene.py
index e93014920ddbcac207e97e07e38c5a6dbadaf3ad..7eabb367818c0e07b05ab75f659af659559a06d5 100644
--- a/doc/source/images/tutorial4.py
+++ b/doc/source/images/4-graphene.py
@@ -7,24 +7,23 @@
 #  - Application of all the aspects of tutorials 1-3 to a more complicated
 #    lattice, namely graphene
 
-from __future__ import division # so that 1/2 == 0.5, and not 0
-from math import pi, sqrt
-import numpy as np
+from __future__ import division  # so that 1/2 == 0.5, and not 0
+from math import pi, sqrt, tanh
 
 import kwant
-import latex, html
 
 # For computing eigenvalues
 import scipy.sparse.linalg as sla
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
+import latex, html
 
 # Define the graphene lattice
-sin_30, cos_30 = (1/2, np.sqrt(3)/2)
+sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
 graphene = kwant.make_lattice([(1, 0), (sin_30, cos_30)],
-                              [(0, 0), (0, 1/np.sqrt(3))])
+                              [(0, 0), (0, 1 / sqrt(3))])
 a, b = graphene.sublattices
 
 
@@ -34,17 +33,17 @@ def make_system(r=10, w=2.0, pot=0.1):
     # circular scattering region
     def circle(pos):
         x, y = pos
-        return x**2 + y**2 < r**2
+        return x ** 2 + y ** 2 < r ** 2
 
-    sys= kwant.Builder()
+    sys = kwant.Builder()
 
     # w: width and pot: potential maximum of the p-n junction
     def potential(site):
         (x, y) = site.pos
         d = y * cos_30 + x * sin_30
-        return pot * np.tanh(d / w)
+        return pot * tanh(d / w)
 
-    sys[graphene.shape(circle, (0,0))] = potential
+    sys[graphene.shape(circle, (0, 0))] = potential
 
     # specify the hoppings of the graphene lattice in the
     # format expected by possibe_hoppings()
@@ -53,8 +52,8 @@ def make_system(r=10, w=2.0, pot=0.1):
         sys[sys.possible_hoppings(*hopping)] = - 1
 
     # Modify the scattering region
-    del sys[a(0,0)]
-    sys[a(-2,1), b(2, 2)] = -1
+    del sys[a(0, 0)]
+    sys[a(-2, 1), b(2, 2)] = -1
 
     #### Define the leads. ####
     # left lead
@@ -65,7 +64,7 @@ def make_system(r=10, w=2.0, pot=0.1):
         return (-1 < x < 1) and (-0.4 * r < y < 0.4 * r)
 
     lead0 = kwant.Builder(sym0)
-    lead0[graphene.shape(lead0_shape, (0,0))] = - pot
+    lead0[graphene.shape(lead0_shape, (0, 0))] = - pot
     for hopping in hoppings:
         lead0[lead0.possible_hoppings(*hopping)] = - 1
 
@@ -79,7 +78,7 @@ def make_system(r=10, w=2.0, pot=0.1):
         return (-1 < u < 1) and (-0.4 * r < v < 0.4 * r)
 
     lead1 = kwant.Builder(sym1)
-    lead1[graphene.shape(lead1_shape, (0,0))] = pot
+    lead1[graphene.shape(lead1_shape, (0, 0))] = pot
     for hopping in hoppings:
         lead1[lead1.possible_hoppings(*hopping)] = - 1
 
@@ -94,34 +93,33 @@ def compute_evs(sys):
         # This requires scipy version >= 0.9.0
         # Failure (i.e. insufficient scipy version) is not critical
         # for the remainder of the tutorial, hence the try-block
-        evs = scipy.sparse.linalg.eigs(sparse_mat, 2)[0]
+        evs = sla.eigs(sparse_mat, 2)[0]
         print evs
     except:
         pass
 
 
-def plot_conductance(fsys, energies):
+def plot_conductance(sys, energies):
     # Compute transmission as a function of energy
     data = []
     for energy in energies:
-        smatrix = kwant.solve(fsys, energy)
+        smatrix = kwant.solve(sys, energy)
         data.append(smatrix.transmission(0, 1))
 
-    pylab.clf()
-    pylab.plot(energies, data)
-    pylab.xlabel("energy [in units of t]",
+    fig = pyplot.figure()
+    pyplot.plot(energies, data)
+    pyplot.xlabel("energy [in units of t]",
                  fontsize=latex.mpl_label_size)
-    pylab.ylabel("conductance [in units of e^2/h]",
+    pyplot.ylabel("conductance [in units of e^2/h]",
                  fontsize=latex.mpl_label_size)
-    fig = pylab.gcf()
-    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
                fontsize=latex.mpl_tick_size)
-    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
                fontsize=latex.mpl_tick_size)
     fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
     fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-    fig.savefig("tutorial4_result.pdf")
-    fig.savefig("tutorial4_result.png",
+    fig.savefig("4-graphene_result.pdf")
+    fig.savefig("4-graphene_result.png",
                 dpi=(html.figwidth_px/latex.mpl_width_in))
 
 
@@ -130,21 +128,20 @@ def plot_bandstructure(flead, momenta):
     # the bandstructure
     energy_list = [flead.energies(k) for k in momenta]
 
-    pylab.clf()
-    pylab.plot(momenta, energy_list)
-    pylab.xlabel("momentum [in untis of (lattice constant)^-1]",
+    fig = pyplot.figure()
+    pyplot.plot(momenta, energy_list)
+    pyplot.xlabel("momentum [in untis of (lattice constant)^-1]",
                  fontsize=latex.mpl_label_size)
-    pylab.ylabel("energy [in units of t]",
+    pyplot.ylabel("energy [in units of t]",
                  fontsize=latex.mpl_label_size)
-    fig = pylab.gcf()
-    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
                fontsize=latex.mpl_tick_size)
-    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
                fontsize=latex.mpl_tick_size)
     fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
     fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-    fig.savefig("tutorial4_bs.pdf")
-    fig.savefig("tutorial4_bs.png",
+    fig.savefig("4-graphene_bs.pdf")
+    fig.savefig("4-graphene_bs.png",
                 dpi=(html.figwidth_px/latex.mpl_width_in))
 
 
@@ -161,9 +158,9 @@ def main():
 
     # Plot the closed system without leads.
     kwant.plot(sys, symbols=plotter_symbols,
-               filename="tutorial4_sys1.pdf", width=latex.figwidth_pt)
+               filename="4-graphene_sys1.pdf", width=latex.figwidth_pt)
     kwant.plot(sys, symbols=plotter_symbols,
-               filename="tutorial4_sys1.png", width=html.figwidth_px)
+               filename="4-graphene_sys1.png", width=html.figwidth_px)
 
     # Compute some eigenvalues.
     compute_evs(sys.finalized())
@@ -174,20 +171,20 @@ def main():
 
     # Then, plot the system with leads.
     kwant.plot(sys, symbols=plotter_symbols,
-               filename="tutorial4_sys2.pdf", width=latex.figwidth_pt)
+               filename="4-graphene_sys2.pdf", width=latex.figwidth_pt)
     kwant.plot(sys, symbols=plotter_symbols,
-               filename="tutorial4_sys2.png", width=html.figwidth_px)
+               filename="4-graphene_sys2.png", width=html.figwidth_px)
 
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
     # Compute the band structure of lead 0.
-    momenta = np.arange(-pi, pi + .01, 0.1 * pi)
-    plot_bandstructure(fsys.leads[0], momenta)
+    momenta = [-pi + 0.02 * pi * i for i in xrange(101)]
+    plot_bandstructure(sys.leads[0], momenta)
 
     # Plot conductance.
-    energies = np.arange(-2 * pot, 2 * pot, pot / 10.5)
-    plot_conductance(fsys, energies)
+    energies = [-2 * pot + 4. / 50. * pot * i for i in xrange(51)]
+    plot_conductance(sys, energies)
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/doc/source/images/tutorial5a.py b/doc/source/images/5-superconductor_band_structure.py
similarity index 63%
rename from doc/source/images/tutorial5a.py
rename to doc/source/images/5-superconductor_band_structure.py
index 359c3a3530cdc91c50207d76662a6a9fd9dcbe38..d2bc32767e3879f48b9df81f364c602aa3c12825 100644
--- a/doc/source/images/tutorial5a.py
+++ b/doc/source/images/5-superconductor_band_structure.py
@@ -11,16 +11,18 @@
 #    in tutorial5b.py
 
 import kwant
-import latex, html
 
 import numpy as np
 from math import pi
 
 # For plotting
-import pylab
+from matplotlib import pyplot
+
+import latex, html
+
+tau_x = np.array([[0, 1], [1, 0]])
+tau_z = np.array([[1, 0], [0, -1]])
 
-tau_x = np.array([[0,1],[1,0]])
-tau_z = np.array([[1,0],[0,-1]])
 
 def make_lead(a=1, t=1.0, mu=0.7, Delta=0.1, W=10):
     # Start with an empty lead with a single square lattice
@@ -28,50 +30,49 @@ def make_lead(a=1, t=1.0, mu=0.7, Delta=0.1, W=10):
 
     sym_lead = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead = kwant.Builder(sym_lead)
-    lead.default_site_group = lat
 
     # build up one unit cell of the lead, and add the hoppings
     # to the next unit cell
     for j in xrange(W):
-        lead[(0, j)] = (4 * t - mu) * tau_z + Delta * tau_x
+        lead[lat(0, j)] = (4 * t - mu) * tau_z + Delta * tau_x
 
         if j > 0:
-            lead[(0, j), (0, j-1)] = - t * tau_z
+            lead[lat(0, j), lat(0, j - 1)] = - t * tau_z
 
-        lead[(1, j), (0, j)] = - t * tau_z
+        lead[lat(1, j), lat(0, j)] = - t * tau_z
 
     return lead
 
 
-def plot_bandstructure(flead, momenta):
+def plot_bandstructure(lead, momenta):
     # Use the method ``energies`` of the finalized lead to compute
     # the bandstructure
-    energy_list = [flead.energies(k) for k in momenta]
-
-    pylab.plot(momenta, energy_list)
-    pylab.xlabel("momentum [in untis of (lattice constant)^-1]")
-    pylab.ylabel("energy [in units of t]")
-    pylab.ylim([-0.8, 0.8])
-    fig = pylab.gcf()
-    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+    energy_list = [lead.energies(k) for k in momenta]
+
+    fig = pyplot.figure()
+    pyplot.plot(momenta, energy_list)
+    pyplot.xlabel("momentum [in untis of (lattice constant)^-1]")
+    pyplot.ylabel("energy [in units of t]")
+    pyplot.ylim([-0.8, 0.8])
+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
                fontsize=latex.mpl_tick_size)
-    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
                fontsize=latex.mpl_tick_size)
     fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
     fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-    fig.savefig("tutorial5a_result.pdf")
-    fig.savefig("tutorial5a_result.png",
+    fig.savefig("5-superconductor_band_structure_result.pdf")
+    fig.savefig("5-superconductor_band_structure_result.png",
                 dpi=(html.figwidth_px/latex.mpl_width_in))
 
 
 def main():
     # Make system and finalize it right away.
-    flead = make_lead().finalized()
+    lead = make_lead().finalized()
 
     # list of momenta at which the bands should be computed
-    momenta = np.arange(-1.5, 1.5 + .0001, 0.002 * pi)
+    momenta = np.linspace(-1.5, 1.5, 201)
 
-    plot_bandstructure(flead, momenta)
+    plot_bandstructure(lead, momenta)
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/doc/source/images/tutorial5b.py b/doc/source/images/5-superconductor_transport.py
similarity index 84%
rename from doc/source/images/tutorial5b.py
rename to doc/source/images/5-superconductor_transport.py
index 1ed8213f8797b1d911938bb064a9948f585b436d..b9a90a61462377473bc6bfc301de5298cabc4306 100644
--- a/doc/source/images/tutorial5b.py
+++ b/doc/source/images/5-superconductor_transport.py
@@ -8,13 +8,11 @@
 #   using different lattices
 
 import kwant
-import latex, html
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
-# For matrix support
-import numpy
+import latex, html
 
 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):
@@ -43,7 +41,7 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
 
     # Superconducting order parameter enters as hopping between
     # electrons and holes
-    sys[((lat_e(x,y), lat_h(x, y)) for x in range(Deltapos, L)
+    sys[((lat_e(x, y), lat_h(x, y)) for x in range(Deltapos, L)
          for y in range(W))] = Delta
 
     #### Define the leads. ####
@@ -87,36 +85,39 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
 
     return sys
 
-def plot_conductance(fsys, energies):
+
+def plot_conductance(sys, energies):
     # Compute conductance
     data = []
     for energy in energies:
-        smatrix = kwant.solve(fsys, energy)
+        smatrix = kwant.solve(sys, energy)
         # Conductance is N - R_ee + R_he
         data.append(smatrix.submatrix(0, 0).shape[0] -
                     smatrix.transmission(0, 0) +
                     smatrix.transmission(1, 0))
 
-    pylab.plot(energies, data)
-    pylab.xlabel("energy [in units of t]")
-    pylab.ylabel("conductance [in units of e^2/h]")
-    fig = pylab.gcf()
-    pylab.setp(fig.get_axes()[0].get_xticklabels(),
+    fig = pyplot.figure()
+    pyplot.plot(energies, data)
+    pyplot.xlabel("energy [in units of t]")
+    pyplot.ylabel("conductance [in units of e^2/h]")
+    pyplot.setp(fig.get_axes()[0].get_xticklabels(),
                fontsize=latex.mpl_tick_size)
-    pylab.setp(fig.get_axes()[0].get_yticklabels(),
+    pyplot.setp(fig.get_axes()[0].get_yticklabels(),
                fontsize=latex.mpl_tick_size)
     fig.set_size_inches(latex.mpl_width_in, latex.mpl_width_in*3./4.)
     fig.subplots_adjust(left=0.15, right=0.95, top=0.95, bottom=0.15)
-    fig.savefig("tutorial5b_result.pdf")
-    fig.savefig("tutorial5b_result.png",
+    fig.savefig("5-superconductor_transport_result.pdf")
+    fig.savefig("5-superconductor_transport_result.png",
                 dpi=(html.figwidth_px/latex.mpl_width_in))
 
 
-
 def main():
-    fsys = make_system().finalized()
+    sys = make_system()
+
+    # Finalize the system.
+    sys = sys.finalized()
 
-    plot_conductance(fsys, energies=[0.002 * i for i in xrange(100)])
+    plot_conductance(sys, energies=[0.002 * i for i in xrange(100)])
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/doc/source/images/tutorial5b_sketch.svg b/doc/source/images/5-superconductor_transport_sketch.svg
similarity index 100%
rename from doc/source/images/tutorial5b_sketch.svg
rename to doc/source/images/5-superconductor_transport_sketch.svg
diff --git a/doc/source/tutorial/index.rst b/doc/source/tutorial/index.rst
index 255772a847dc62251969c2acebaa7ff14ffa7b8e..883dec65b7eb62b776e372293d1ae6af3048d153 100644
--- a/doc/source/tutorial/index.rst
+++ b/doc/source/tutorial/index.rst
@@ -5,7 +5,7 @@ In the following, the most important features of kwant are explained using
 simple, but still physically meaningful examples. Each of the examples
 is commented extensively. In addition, you will find notes about more subtle,
 technical details at the end of each example. At first reading,
-these notes maybe safely skipped.
+these notes may be safely skipped.
 
 .. toctree::
     introduction
diff --git a/doc/source/tutorial/tutorial1.rst b/doc/source/tutorial/tutorial1.rst
index 697b58780bc96517fee6f4ad3ba9e3382c1bc7bf..7e02b6a839d0811f962f52761992cfedf3bc394b 100644
--- a/doc/source/tutorial/tutorial1.rst
+++ b/doc/source/tutorial/tutorial1.rst
@@ -16,7 +16,7 @@ with a hard wall confinement :math:`V(y)` in y-direction.
 
 In order to use kwant, we need to import it:
 
-.. literalinclude:: ../../../examples/tutorial1a.py
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
     :lines: 11
 
 Enabling kwant is as easy as this [#]_ !
@@ -26,42 +26,40 @@ 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:: ../../../examples/tutorial1a.py
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
     :lines: 15
 
-Next, we tell `~kwant.builder.Builder` that we want to work
-with a square lattice (more about the details of this code snippet in
-the notes below).  For simplicity, we set the lattice constant to
-unity:
+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:: ../../../examples/tutorial1a.py
-    :lines: 18-20
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
+    :lines: 18-19
 
 Since we work with a square lattice, we label the points with two
 integer coordinates `(i, j)`. `~kwant.builder.Builder` then
 allows us to add matrix elements corresponding to lattice points:
-``sys[(i, j)] = ...`` sets the on-site energy for the point `(i, j)`,
-and ``sys[(i1, j1), (i2, j2)] = ...`` the hopping matrix element
+``sys[lat(i, j)] = ...`` sets the on-site energy for the point `(i, j)`,
+and ``sys[lat(i1, j1), lat(i2, j2)] = ...`` the hopping matrix element
 **from** point `(i2, j2)` **to** point `(i1, j1)`.
 
+Note that we need to specify sites for `~kwant.builder.Builder`
+in the form ``lat(i, j)``. The lattice object `lat` does the
+translation from integer coordinates to proper site format
+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:: ../../../examples/tutorial1a.py
-    :lines: 22-24, 27-38
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
+    :lines: 21-23, 26-37
 
 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:: ../../../examples/tutorial1a.py
-    :lines: 46-48
-
-.. note::
-
-    Here it is essential that we write ``lead0.default_site_group = lat``
-    instead of ``lead0.default_site_group = kwant.lattice.Square(a)``.
-    For details see the notes below.
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
+    :lines: 45-46
 
 Here, the `~kwant.builder.Builder` takes the translational symmetry
 as an optional parameter. Note that the (real space)
@@ -75,8 +73,8 @@ 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:: ../../../examples/tutorial1a.py
-    :lines: 50-56
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
+    :lines: 48-54
 
 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.
@@ -85,8 +83,8 @@ 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:: ../../../examples/tutorial1a.py
-    :lines: 60-70
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
+    :lines: 57-67
 
 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`
@@ -96,8 +94,8 @@ 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:: ../../../examples/tutorial1a.py
-    :lines: 74-75
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
+    :lines: 71-72
 
 More details about attaching leads can be found in the tutorial
 :ref:`tutorial-abring`.
@@ -105,12 +103,12 @@ 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:: ../../../examples/tutorial1a.py
-    :lines: 79
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
+    :lines: 76
 
 This should bring up this picture:
 
-.. image:: /images/tutorial1a_sys.*
+.. image:: /images/1-quantum_wire_sys.*
 
 The system is represented in the usual way for tight-binding systems:
 dots represent the lattice points `(i, j)`, and for every
@@ -120,18 +118,18 @@ fading color.
 
 In order to use our system for a transport calculation, we need to finalize it
 
-.. literalinclude:: ../../../examples/tutorial1a.py
-    :lines: 83
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
+    :lines: 80
 
 Having successfully created a system, we now can immediately start to compute
 its conductance as a function of energy:
 
- .. literalinclude:: ../../../examples/tutorial1a.py
-    :lines: 87-98
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
+    :lines: 84-95
 
-Currently, there is only one algorithm implemented to compute the
-conductance: :func:`kwant.solve <kwant.solvers.sparse.solve>` which computes
-the scattering matrix `smatrix` solving a sparse linear system.
+Currently **CHANGE**, there is only one algorithm implemented to compute the
+conductance: :func:`kwant.solve <kwant.solvers.common.SparseSolver.solve>`
+which computes the scattering matrix `smatrix` solving a sparse linear system.
 `smatrix` itself allows you to directly compute the total
 transmission probability from lead 0 to lead 1 as
 ``smatrix.transmission(1, 0)``.
@@ -140,12 +138,12 @@ 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:: ../../../examples/tutorial1a.py
-    :lines: 102-108
+.. literalinclude:: ../../../tutorial/1-quantum_wire.py
+    :lines: 99-105
 
 This should yield the result
 
-.. image:: /images/tutorial1a_result.*
+.. image:: /images/1-quantum_wire_result.*
 
 We see a conductance quantized in units of :math:`e^2/h`,
 increasing in steps as the energy is increased. The
@@ -155,43 +153,35 @@ subbands that increases with energy.
 
 .. seealso::
      The full source code can be found in
-     :download:`example/tutorial1a.py <../../../examples/tutorial1a.py>`
+     :download:`tutorial/1-quantum_wire.py <../../../tutorial/1-quantum_wire.py>`
 
 .. specialnote:: Technical details
 
    - In the example above, when building the system, only one direction
-     of hopping is given, i.e. ``sys[(i, j), (i, j-1)] = ...`` and
-     not also ``sys[(i, j-1), (i, j)] = ...``. The reason is that
+     of hopping is given, i.e. ``sys[lat(i, j), lat(i, j-1)] = ...`` and
+     not also ``sys[lat(i, j-1), lat(i, j)] = ...``. The reason is that
      `~kwant.builder.Builder` automatically adds the other
      direction of the hopping such that the resulting system is Hermitian.
 
      However, it does not hurt to define the opposite direction of hopping as
-     well.
+     well::
 
-         sys[(1, 0), (0, 0)] = - t
-         sys[(0, 0), (1, 0)] = - t.conj()
+         sys[lat(1, 0), lat(0, 0)] = - t
+         sys[lat(0, 0), lat(1, 0)] = - t.conj()
 
      (assuming that `t` is complex) is perfectly fine. However,
      be aware that also
 
      ::
 
-         sys[(1, 0), (0, 0)] = - 1
-         sys[(0, 0), (1, 0)] = - 2
-
-     is valid code. In the latter case, the hopping ``sys[(1, 0), (0, 0)]``
-     is overwritten by the last line and also equals to -2.
-
-   - Some more details about
-
-     ::
+         sys[lat(1, 0), lat(0, 0)] = - 1
+         sys[lat(0, 0), lat(1, 0)] = - 2
 
-         lat = kwant.lattices.Square(a)
-         sys.default_site_group = lat
+     is valid code. In the latter case, the hopping ``sys[lat(1, 0),
+     lat(0, 0)]`` is overwritten by the last line and also equals to -2.
 
-     By setting ``sys.default_site_group = lat`` you specify to
-     `~kwant.builder.Builder` that it should interpret tuples like
-     `(i, j)` as indices in a square lattice.
+   - Some more details the relation between `~kwant.builder.Builder`
+     and the square lattice `lat` in the example:
 
      Technically, `~kwant.builder.Builder` expects
      **sites** as indices. Sites themselves have a certain type, and
@@ -200,48 +190,36 @@ subbands that increases with energy.
      proper `~kwant.builder.Site` object that can be used with
      `~kwant.builder.Builder`.
 
-     In the above example, `lat` is the site group. By specifying it
-     as the `default_site_group`, `~kwant.builder.Builder`
-     knows that it should use `lat` to interpret any input that is not of
-     type `~kwant.builder.Site`. Instead of using
-     `default_site_group`, one could have manually converted the
-     tuples `(i, j)` into sites ``lat(i, j)``::
-
-         for i in xrange(L):
-             for j in xrange(W):
-                 sys[lat(i, j)] = 4 * t
-
-                 # hoppig in y-direction
-                 if j > 0 :
-                     sys[lat(i, j), lat(i, j-1)] = - t
+     In the above example, `lat` is the site group. ``lat(i, j)``
+     then translates the description of a lattice site in terms of two
+     integer indices (which is the natural way to do here) into
+     a proper `~kwant.builder.Site` object.
 
-                 #hopping in x-direction
-                 if i > 0:
-                     sys[lat(i, j), lat(i-1, j)] = -t
+     The concept of site groups and sites allows `~kwant.builder.Builder`
+     to mix arbitrary lattices and site groups
 
-     (The concept of site groups and sites allows `~kwant.builder.Builder`
-     to mix arbitrary lattices and site groups)
+   - In the example, we wrote
 
-   - Note that we wrote::
-
-         lat = kwant.lattices.Square(a)
+     ::
 
-         sys.default_site_group = lat
-         lead0.default_site_group = lat
+         sys = sys.finalized()
 
-     instead of::
+     In doing so, we transform the `~kwant.builder.Builder` object (with which
+     we built up the system step by step) into a `~kwant.system.System`
+     that has a fixed structure (which we cannot change any more).
 
-         sys.default_site_group = kwant.lattices.Square(a)
-         lead0.default_site_group = kwant.lattices.Square(a)
+     Note that this means that we cannot access the `~kwant.builder.Builder`
+     object any more. This is not necesarry any more, as the computational
+     routines all expect finalized systems. It even has the advantage
+     that python is now free to release the memory occupied by the
+     `~kwant.builder.Builder` which, for large systems, can be considerable.
+     Roughly speaking, the above code corresponds to
 
-     The reason is that in the latter case, `sys` and `lead0` have two
-     different site groups (although both representing a
-     square lattice), since a site group is represented by a particular
-     instance of the class, not the class itself.
+     ::
 
-     Hence, the latter example is interpreted as two different
-     square lattices, which will fail when the lead is attached to the
-     system.
+	 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
@@ -259,7 +237,8 @@ subbands that increases with energy.
      with the lattice symmetry.
 
    - Instead of plotting to the screen (which is standard, if the
-     Python Image Library PIL is installed), :func:`plot <kwant.plotter.plot>`
+     Python Image Library PIL **CHANGE [if plotter is changed]** is installed),
+     :func:`plot <kwant.plotter.plot>`
      can also write to the file specified by the argument `filename`.
      (for details, see the documentation of :func:`plot <kwant.plotter.plot>`.)
 
@@ -288,13 +267,13 @@ 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:: ../../../examples/tutorial1b.py
+.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
     :lines: 13-24
 
 Previously, the scattering region was build using two ``for``-loops.
 Instead, we now write:
 
-.. literalinclude:: ../../../examples/tutorial1b.py
+.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
     :lines: 27
 
 Here, all lattice points are added at once in the first line.  The
@@ -311,7 +290,7 @@ hoppings. In this case, an iterable like for the lattice
 points becomes a bit cumbersome, and we use instead another
 feature of kwant:
 
-.. literalinclude:: ../../../examples/tutorial1b.py
+.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
     :lines: 28-29
 
 In regular lattices, one has only very few types of different hoppings
@@ -329,8 +308,8 @@ then sets all of those hopping matrix elements at once.
 
 The leads can be constructed in an analogous way:
 
-.. literalinclude:: ../../../examples/tutorial1b.py
-    :lines: 35-41
+.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
+    :lines: 35-40
 
 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
@@ -341,34 +320,34 @@ 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:: ../../../examples/tutorial1b.py
-    :lines: 45
+.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
+    :lines: 44
 
 The remainder of the code is identical to the previous example
 (except for a bit of reorganization into functions):
 
-.. literalinclude:: ../../../examples/tutorial1b.py
-    :lines: 48-52
+.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
+    :lines: 47-51
 
 and
 
-.. literalinclude:: ../../../examples/tutorial1b.py
-    :lines: 53-63
+.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
+    :lines: 52-64
 
 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:: ../../../examples/tutorial1b.py
-    :lines: 66-76
+.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
+    :lines: 67-77
 
 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:: ../../../examples/tutorial1b.py
-    :lines: 81-82
+.. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
+    :lines: 82-83
 
 If the example however is imported using ``import tutorial1b``,
 ``main()`` is not executed automatically. Instead, you can execute it
@@ -380,13 +359,13 @@ The result of the example should be identical to the previous one.
 
 .. seealso::
     The full source code can be found in
-    :download:`examples/tutorial1b.py <../../../examples/tutorial1b.py>`
+    :download:`tutorial/1-quantum_wire_revisited.py <../../../tutorial/1-quantum_wire_revisited.py>`
 
 .. specialnote:: Technical details
 
    - In
 
-     .. literalinclude:: ../../../examples/tutorial1b.py
+     .. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
        :lines: 28-29
 
      we write ``*hopping`` instead of ``hopping``. The reason is as follows:
@@ -407,9 +386,7 @@ The result of the example should be identical to the previous one.
    - We have seen different ways to add lattice points to a
      `~kwant.builder.Builder`. It allows to
 
-     * add single points, specified as sites (or tuples, if
-       a `default_site_group` is specified as in the previous
-       example).
+     * add single points, specified as sites
      * add several points at once using a generator (as in this example)
      * add several points at once using a list (typically less
        effective compared to a generator)
@@ -418,10 +395,10 @@ 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:: ../../../examples/tutorial1b.py
+     .. literalinclude:: ../../../tutorial/1-quantum_wire_revisited.py
          :lines: 27
 
-     Note that ``((x, y) for x in range(L) for y in range(W))`` is not
+     Note that ``(lat(x, y) for x in range(L) for y in range(W))`` is not
      a tuple, but a generator.
 
      Let us elaborate a bit more on this using a simpler example:
@@ -453,7 +430,7 @@ The result of the example should be identical to the previous one.
      added at once, these two sites should be encapsulated in a tuple.
      In particular, one must write::
 
-         sys[(((0,j+1), (0, j)) for j in xrange(W-1)] = ...
+         sys[((lat(0,j+1), lat(0, j)) for j in xrange(W-1)] = ...
 
      or::
 
diff --git a/doc/source/tutorial/tutorial2.rst b/doc/source/tutorial/tutorial2.rst
index eba8cc01b130166b9b2e3ffbb9f0682be0d2ed18..e0046e16de36b6eb9ce9b5ff1ddb45cca2cc5e3c 100644
--- a/doc/source/tutorial/tutorial2.rst
+++ b/doc/source/tutorial/tutorial2.rst
@@ -37,21 +37,21 @@ 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:: ../../../examples/tutorial2a.py
+.. literalinclude:: ../../../tutorial/2-spin_orbit.py
     :lines: 19
 
 For convenience, we define the Pauli-matrices first (with `sigma_0` the
 unit matrix):
 
-.. literalinclude:: ../../../examples/tutorial2a.py
+.. literalinclude:: ../../../tutorial/2-spin_orbit.py
     :lines: 22-25
 
 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:: ../../../examples/tutorial2a.py
-    :lines: 37-44
+.. literalinclude:: ../../../tutorial/2-spin_orbit.py
+    :lines: 36-43
 
 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,17 +78,17 @@ when specifying `(1, 0)` it is not necessary to specify `(-1, 0)`),
 
 The leads also allow for a matrix structure,
 
-.. literalinclude:: ../../../examples/tutorial2a.py
-    :lines: 52-58
+.. literalinclude:: ../../../tutorial/2-spin_orbit.py
+    :lines: 50-56
 
 The remainder of the code is unchanged, and as a result we should obtain
 the following, clearly non-monotonic conductance steps:
 
-.. image:: ../images/tutorial2a_result.*
+.. image:: ../images/2-spin_orbit_result.*
 
 .. seealso::
      The full source code can be found in
-     :download:`example/tutorial2a.py <../../../examples/tutorial2a.py>`
+     :download:`tutorial/2-spin_orbit.py <../../../tutorial/2-spin_orbit.py>`
 
 .. specialnote:: Technical details
 
@@ -124,8 +124,8 @@ 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:: ../../../examples/tutorial2b.py
-    :lines: 16-18, 22, 28-34
+.. literalinclude:: ../../../tutorial/2-quantum_well.py
+    :lines: 16-18, 21-24, 27-33
 
 This function takes one argument which is of type
 `~kwant.builder.Site`, from which you can get the realspace
@@ -140,8 +140,8 @@ the transmission as a function of well depth.
 kwant now allows us to pass a function as a value to
 `~kwant.builder.Builder`:
 
-.. literalinclude:: ../../../examples/tutorial2b.py
-    :lines: 36-41
+.. literalinclude:: ../../../tutorial/2-quantum_well.py
+    :lines: 35-40
 
 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 +155,8 @@ of the lead -- this should be kept in mind.
 
 Finally, we compute the transmission probability:
 
-.. literalinclude:: ../../../examples/tutorial2b.py
-    :lines: 65, 68-77
+.. literalinclude:: ../../../tutorial/2-quantum_well.py
+    :lines: 67, 70-82
 
 Since we change the value of the global variable `pot` to vary the
 well depth, python requires us to write ``global pot`` to `enable
@@ -165,14 +165,14 @@ access to it
 Subsequent calls to :func:`kwant.solve <kwant.solvers.sparse.solve>`
 then will use the updated value of pot, and we get the result:
 
-.. image:: ../images/tutorial2b_result.*
+.. image:: ../images/2-quantum_well_result.*
 
 Starting from no potential (well depth = 0), we observe the typical
 oscillatory transmission behavior through resonances in the quantum well.
 
 .. seealso::
      The full source code can be found in
-     :download:`example/tutorial2b.py <../../../examples/tutorial2b.py>`
+     :download:`tutorial/2-quantum_well.py <../../../tutorial/2-quantum_well.py>`
 
 .. warning::
 
@@ -185,9 +185,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 example/tutorial2b.py, line 16
+  - In tutorial/2-quantum_well.py, line 16
 
-    .. literalinclude:: ../../../examples/tutorial2b.py
+    .. literalinclude:: ../../../tutorial/2-quantum_well.py
         :lines: 16
 
     is not really necessary. If this line was left out, the
@@ -276,7 +276,7 @@ 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`:
 
-.. image:: ../images/tutorial2c_sketch.*
+.. image:: ../images/2-ab_ring_sketch.*
 
 For a flux line, it is possible to choose a gauge such that a
 charged particle acquires a phase :math:`e\Phi/h` whenever it
@@ -290,8 +290,8 @@ 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:: ../../../examples/tutorial2c.py
-    :lines: 20, 24-27, 30-33
+.. literalinclude:: ../../../tutorial/2-ab_ring.py
+    :lines: 21, 24-27, 31-34
 
 Note that this function takes a realspace position as argument (not a
 `~kwant.builder.Site`).
@@ -300,8 +300,8 @@ 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:: ../../../examples/tutorial2c.py
-    :lines: 36-38
+.. literalinclude:: ../../../tutorial/2-ab_ring.py
+    :lines: 37-39
 
 Here, ``lat.shape()`` takes as a second parameter a (realspace) point
 that is inside the desired shape. The hoppings can still be added
@@ -311,11 +311,11 @@ Up to now, the system contains constant hoppings and onsite energies,
 and we still need to include the phase shift due to the magnetic flux.
 This is done by **overwriting** the values of hoppings in x-direction
 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 `((1, j), (0, j))`
+all hoppings in x-direction that are of the form `(lat(1, j), lat(0, j))`
 with ``j<0``:
 
-.. literalinclude:: ../../../examples/tutorial2c.py
-    :lines: 46-58
+.. literalinclude:: ../../../tutorial/2-ab_ring.py
+    :lines: 47-59
 
 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 +328,16 @@ by the global variable `phi`.
 
 For the leads, we can also use the ``lat.shape()``-functionality:
 
-.. literalinclude:: ../../../examples/tutorial2c.py
-    :lines: 62-71
+.. literalinclude:: ../../../tutorial/2-ab_ring.py
+    :lines: 63-72
 
 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:: ../../../examples/tutorial2c.py
-    :lines: 78-79
+.. literalinclude:: ../../../tutorial/2-ab_ring.py
+    :lines: 79-80
 
 In fact, attaching leads seems not so simple any more for the current
 structure with a scattering region very much different from the lead
@@ -350,23 +350,23 @@ back (going opposite to the direction of the translational vector)
 until it intersects the scattering region. At this intersection,
 the lead is attached:
 
-.. image:: ../images/tutorial2c_sketch2.*
+.. image:: ../images/2-ab_ring_sketch2.*
 
 After the lead has been attached, the system should look like this:
 
-.. image:: ../images/tutorial2c_sys.*
+.. image:: ../images/2-ab_ring_sys.*
 
 The computation of the conductance goes in the same fashion as before.
 Finally you should get the following result:
 
-.. image:: ../images/tutorial2c_result.*
+.. image:: ../images/2-ab_ring_result.*
 
 where one can observe the conductance oscillations with the
 period of one flux quantum.
 
 .. seealso::
      The full source code can be found in
-     :download:`example/tutorial2c.py <../../../examples/tutorial2c.py>`
+     :download:`tutorial/2-ab_ring.py <../../../tutorial/2-ab_ring.py>`
 
 .. specialnote:: Technical details
 
@@ -384,7 +384,7 @@ period of one flux quantum.
     becomes more apparent if we attach the leads a bit further away
     from the central axis o the ring, as was done in this example:
 
-    .. image:: ../images/tutorial2c_note1.*
+    .. image:: ../images/2-ab_ring_note1.*
 
   - Per default, `~kwant.builder.Builder.attach_lead` attaches
     the lead to the "outside" of the structure, by tracing the
@@ -399,7 +399,7 @@ period of one flux quantum.
     starts the trace-back in the middle of the ring, resulting
     in the lead being attached to the inner circle:
 
-    .. image:: ../images/tutorial2c_note2.*
+    .. image:: ../images/2-ab_ring_note2.*
 
     Note that here the lead is treated as if it would pass over
     the other arm of the ring, without intersecting it.
diff --git a/doc/source/tutorial/tutorial3.rst b/doc/source/tutorial/tutorial3.rst
index 6e5b8e0b953105d2214d3e8a912b7e17ed7e5926..5907ff947a7541830f8d2f2454bd8b45abdb5a75 100644
--- a/doc/source/tutorial/tutorial3.rst
+++ b/doc/source/tutorial/tutorial3.rst
@@ -15,8 +15,8 @@ tight-binding wire.
 Computing band structures in kwant is easy. Just define a lead in the
 usual way:
 
-.. literalinclude:: ../../../examples/tutorial3a.py
-    :lines: 18-37
+.. literalinclude:: ../../../tutorial/3-band_structure.py
+    :lines: 17-34
 
 "Usual way" means defining a translational symmetry vector, as well
 as one unit cell of the lead, and the hoppings to neighboring
@@ -25,16 +25,17 @@ invariant system needed for band structure calculations.
 
 In contrast to previous usage however, you have to *finalize* the lead.  A
 finalized lead has a method/attribute `~kwant.system.InfiniteSystem.energies`
+***CHANGE [once energies is in its own band structure module] ***
 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:: ../../../examples/tutorial3a.py
-    :lines: 40 - 57
+.. literalinclude:: ../../../tutorial/3-band_structure.py
+    :lines: 37 - 55
 
 This gives the result:
 
-.. image:: ../images/tutorial3a_result.*
+.. image:: ../images/3-band_structure_result.*
 
 where we observe the cosine-like dispersion of the square lattice. Close
 to ``k=0`` this agrees well with the quadratic dispersion this tight-binding
@@ -42,7 +43,7 @@ Hamiltonian is approximating.
 
 .. seealso::
      The full source code can be found in
-     :download:`example/tutorial3a.py <../../../examples/tutorial3a.py>`
+     :download:`tutorial/3-band_structure.py <../../../tutorial/3-band_structure.py>`
 
 .. specialnote:: Technical details
 
@@ -74,31 +75,31 @@ 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:: ../../../examples/tutorial3b.py
+.. literalinclude:: ../../../tutorial/3-closed_system.py
     :lines: 16
 
 We set up the system using the `shape`-function as in
 :ref:`tutorial-abring`, but do not add any leads:
 
-.. literalinclude:: ../../../examples/tutorial3b.py
-    :lines: 30-47
+.. literalinclude:: ../../../tutorial/3-closed_system.py
+    :lines: 26-48
 
 We add the magnetic field using a function and a global variable as we
-did in the two previous examples. (Here, the gauge is chosen such that
+did in the two previous tutorial. (Here, the gauge is chosen such that
 :math:`A_x(y) = - B y` and :math:`A_y=0`.)
 
 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:: ../../../examples/tutorial3b.py
-    :lines: 50, 52, 58-69
+.. literalinclude:: ../../../tutorial/3-closed_system.py
+    :lines: 51, 53, 60-76
 
 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
 want to use sparse matrix methods.) Finally, we obtain the result:
 
-.. image:: ../images/tutorial3b_result.*
+.. image:: ../images/3-closed_system_result.*
 
 At zero magnetic field several energy levels are degenerate (since our
 quantum dot is rather symmetric). These degeneracies are split
@@ -107,7 +108,7 @@ Landau level energies at higher magnetic fields [#]
 
 .. seealso::
     The full source code can be found in
-    :download:`examples/tutorial3b.py <../../../examples/tutorial3b.py>`
+    :download:`tutorial/3-closed_system.py <../../../tutorial/3-closed_system.py>`
 
 .. specialnote:: Technical details
 
diff --git a/doc/source/tutorial/tutorial4.rst b/doc/source/tutorial/tutorial4.rst
index 85ae898e7ee5b2c19048197ca5ec09882cffe7f3..a93e4ae4bf7ed9bdaeded80b7a6d011f11bf0bea 100644
--- a/doc/source/tutorial/tutorial4.rst
+++ b/doc/source/tutorial/tutorial4.rst
@@ -14,7 +14,7 @@ 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:: ../../../examples/tutorial4.py
+.. literalinclude:: ../../../tutorial/4-graphene.py
     :lines: 24-27
 
 The first argument to the `~kwant.lattice.make_lattice` function is the list of
@@ -26,12 +26,12 @@ 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:: ../../../examples/tutorial4.py
-    :lines: 30-31, 34-39, 41-46
+.. literalinclude:: ../../../tutorial/4-graphene.py
+    :lines: 29-30, 33-38, 40-45
 
 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;
-in contrast to the previous examples, the potential value is this time taken
+in contrast to the previous tutorial, the potential value is this time taken
 from the scope of `make_system()`, since we keep the potential fixed
 in this example.
 
@@ -40,8 +40,8 @@ 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:: ../../../examples/tutorial4.py
-    :lines: 50
+.. literalinclude:: ../../../tutorial/4-graphene.py
+    :lines: 49
 
 The nearest-neighbor model for graphene contains only
 hoppings between different basis atoms. For these type of
@@ -57,24 +57,24 @@ respect to the two primitive vectors ``[(1, 0), (sin_30, cos_30)]``.
 
 Adding the hoppings however still works the same way:
 
-.. literalinclude:: ../../../examples/tutorial4.py
-    :lines: 51-52
+.. literalinclude:: ../../../tutorial/4-graphene.py
+    :lines: 50-51
 
 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:: ../../../examples/tutorial4.py
-    :lines: 55-56
+.. literalinclude:: ../../../tutorial/4-graphene.py
+    :lines: 54-55
 
-Note that the conversion from a tuple `(i,j)` to site
+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:: ../../../examples/tutorial4.py
-    :lines: 60-83
+.. literalinclude:: ../../../tutorial/4-graphene.py
+    :lines: 58-82
 
 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 +85,14 @@ 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:: ../../../examples/tutorial4.py
-    :lines: 85
+.. literalinclude:: ../../../tutorial/4-graphene.py
+    :lines: 84
 
 The computation of some eigenvalues of the closed system is done
 in the following piece of code:
 
-.. literalinclude:: ../../../examples/tutorial4.py
-    :lines: 88-93, 96-99
+.. literalinclude:: ../../../tutorial/4-graphene.py
+    :lines: 87-91, 95-98
 
 Here we use in contrast to the previous example a sparse matrix and
 the sparse linear algebra functionality of scipy (this requires
@@ -106,16 +106,16 @@ to the previous examples, and needs not be further explained here.
 Finally, in the `main()` function we make and
 plot the system:
 
-.. literalinclude:: ../../../examples/tutorial4.py
-    :lines: 126-129, 138
+.. literalinclude:: ../../../tutorial/4-graphene.py
+    :lines: 127-130, 139
 
 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.
 
-.. literalinclude:: ../../../examples/tutorial4.py
-    :lines: 132-135
+.. literalinclude:: ../../../tutorial/4-graphene.py
+    :lines: 133-136
 
 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
@@ -131,12 +131,12 @@ well-defined part of the plot.
 
 Plotting the closed system gives this result:
 
-.. image:: ../images/tutorial4_sys1.*
+.. image:: ../images/4-graphene_sys1.*
 
 Computing the eigenvalues of largest magnitude,
 
-.. literalinclude:: ../../../examples/tutorial4.py
-    :lines: 141
+.. literalinclude:: ../../../tutorial/4-graphene.py
+    :lines: 142
 
 should yield two eigenvalues similar to `[ 3.07869311 +1.02714523e-17j,
 -3.06233144 -6.68085759e-18j]` (round-off might change the imaginary part which
@@ -145,11 +145,11 @@ would be equal to zero for exact arithmetics).
 The remaining code of `main()` attaches the leads to the system and plots it
 again:
 
-.. image:: ../images/tutorial4_sys2.*
+.. image:: ../images/4-graphene_sys2.*
 
 It computes the band structure of one of lead 0:
 
-.. image:: ../images/tutorial4_bs.*
+.. image:: ../images/4-graphene_bs.*
 
 showing all the features of a zigzag lead, including the flat
 edge state bands (note that the band structure is not symmetric around
@@ -157,14 +157,14 @@ zero energy, as we have a potential in the leads).
 
 Finally the transmission through the system is computed,
 
-.. image:: ../images/tutorial4_result.*
+.. image:: ../images/4-graphene_result.*
 
 showing the typical resonance-like transmission probability through
 an open quantum dot
 
 .. seealso::
     The full source code can be found in
-    :download:`examples/tutorial4.py <../../../examples/tutorial4.py>`
+    :download:`tutorial/4-graphene.py <../../../tutorial/4-graphene.py>`
 
 .. specialnote:: Technical details
 
diff --git a/doc/source/tutorial/tutorial5.rst b/doc/source/tutorial/tutorial5.rst
index 6b77648fafdddd3e1547db7609c8ba14572dbb6e..3617286f8f6977cfb3b0a7377b674feac9ad2b84 100644
--- a/doc/source/tutorial/tutorial5.rst
+++ b/doc/source/tutorial/tutorial5.rst
@@ -1,5 +1,3 @@
-.. _tutorial-superconductor:
-
 Superconductors: orbital vs lattice degrees of freedom
 ------------------------------------------------------
 
@@ -30,7 +28,7 @@ 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:: ../../../examples/tutorial5a.py
+.. literalinclude:: ../../../tutorial/5-superconductor_band_structure.py
     :lines: 21-42
 
 As you see, the example is syntactically equivalent to our
@@ -39,14 +37,14 @@ is now that the Pauli matrices act in electron-hole space.
 
 Computing the band structure then yields the result
 
-.. image:: ../images/tutorial5a_result.*
+.. image:: ../images/5-superconductor_band_structure_result.*
 
 We clearly observe the superconducting gap in the spectrum. That was easy,
 he?
 
 .. seealso::
     The full source code can be found in
-    :download:`examples/tutorial5a.py <../../../examples/tutorial5a.py>`
+    :download:`tutorial/5-superconductor_band_structure.py <../../../tutorial/5-superconductor_band_structure.py>`
 
 
 "Lattice description": Using different lattices
@@ -70,21 +68,20 @@ to electrons in the normal lead, and :math:`R_\text{eh}` the total
 probability of reflection from electrons to holes in the normal
 lead. However, the current version of kwant does not allow for an easy
 and elegant partitioning of the scattering matrix in these two degrees
-of freedom (well, there is one since v0.1.3, see the technical notes
-below).
+of freedom [#]_.
 
 In the following, we will circumvent this problem by introducing
 separate "leads" for electrons and holes, making use of different
 lattices. The system we consider consists of a normal lead on the left,
 a superconductor on the right, and a tunnel barrier inbetween:
 
-.. image:: ../images/tutorial5b_sketch.*
+.. image:: ../images/5-superconductor_transport_sketch.*
 
 As already mentioned above, we begin by introducing two different
 square lattices representing electron and hole degrees of freedom:
 
-.. literalinclude:: ../../../examples/tutorial5b.py
-    :lines: 18-19,17,23-24
+.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
+    :lines: 16-17,15,20-21
 
 Any diagonal entry (kinetic energy, potentials, ...) in the BdG
 Hamiltonian then corresponds to on-site energies or hoppings within
@@ -92,8 +89,8 @@ the *same* lattice, whereas any off-diagonal entry (essentially, the
 superconducting order parameter :math:`\Delta`) corresponds
 to a hopping between *different* lattices:
 
-.. literalinclude:: ../../../examples/tutorial5b.py
-    :lines: 25-46
+.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
+    :lines: 23-44
 
 Note that the tunnel barrier is added by overwriting previously set
 on-site matrix elements.
@@ -106,8 +103,8 @@ part. We use this fact to attach purely electron and hole leads
 (comprised of only electron *or* hole lattices) to the
 system:
 
-.. literalinclude:: ../../../examples/tutorial5b.py
-    :lines: 49-65
+.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
+    :lines: 47-63
 
 This separation into two different leads allows us then later to compute the
 reflection probablities between electrons and holes explicitely.
@@ -115,15 +112,15 @@ 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:: ../../../examples/tutorial5b.py
-    :lines: 70-80
+.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
+    :lines: 68-78
 
 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:: ../../../examples/tutorial5b.py
-    :lines: 83-87
+.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
+    :lines: 81-85
 
 When computing the conductance, we can now extract reflection from
 electrons to electrons as ``smatrix.transmission(0, 0)`` (Don't get
@@ -131,8 +128,8 @@ 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:: ../../../examples/tutorial5b.py
-    :lines: 89-97
+.. literalinclude:: ../../../tutorial/5-superconductor_transport.py
+    :lines: 88-96
 
 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
@@ -140,7 +137,7 @@ within (electron) lead 0, and from its size we can extract the number of modes
 
 Finally, for the default parameters, we obtain the following result:
 
-.. image:: ../images/tutorial5b_result.*
+.. image:: ../images/5-superconductor_transport_result.*
 
 We a see a conductance that is proportional to the square of the tunneling
 probability within the gap, and proportional to the tunneling probability
@@ -148,7 +145,7 @@ above the gap. At the gap edge, we observe a resonant Andreev reflection.
 
 .. seealso::
     The full source code can be found in
-    :download:`examples/tutorial5b.py <../../../examples/tutorial5b.py>`
+    :download:`tutorial/5-superconductor_transport.py <../../../tutorial/5-superconductor_transport.py>`
 
 .. specialnote:: Technical details
 
@@ -161,11 +158,18 @@ above the gap. At the gap edge, we observe a resonant Andreev reflection.
     - It is in fact possible to separate electron and hole degrees of
       freedom in the scattering matrix, even if one uses matrices for
       these degrees of freedom. In the solve step,
-      `~kwant.solvers.sparse.solve` returns an array containing the
-      transverse wave functions of the lead modes. By inspecting the wave
-      functions, electron and hole wave functions can be distinguished (they
-      only have entries in either the electron part *or* the hole part. If you
-      encounter modes with entries in both parts, you hit a very unlikely
-      situation in which the standard procedure to compute the modes gave you
-      a superposition of electron and hole modes. That is still OK for
+      `~kwant.solvers.common.SparseSolver.solve` returns an array containing
+      the transverse wave functions of the lead modes (in
+      `BlockResult.lead_info <kwant.solvers.common.BlockResult.lead_info>`.
+      By inspecting the wave functions, electron and hole wave
+      functions can be distinguished (they only have entries in either
+      the electron part *or* the hole part. If you encounter modes
+      with entries in both parts, you hit a very unlikely situation in
+      which the standard procedure to compute the modes gave you a
+      superposition of electron and hole modes. That is still OK for
       computing particle current, but not for electrical current).
+
+.. rubric:: Footnotes
+
+.. [#] Well, there is a not so elegant way to do it still. See the technical
+       details
diff --git a/examples/tutorial1a.py b/tutorial/1-quantum_wire.py
similarity index 72%
rename from examples/tutorial1a.py
rename to tutorial/1-quantum_wire.py
index 0919347329eddc85eb66281169e65450b7660c6a..02935c56ec61b529b93b8e6fdc93b8dd9de719bd 100644
--- a/examples/tutorial1a.py
+++ b/tutorial/1-quantum_wire.py
@@ -17,7 +17,6 @@ sys = kwant.Builder()
 # Here, we are only working with square lattices
 a = 1
 lat = kwant.lattice.Square(a)
-sys.default_site_group = lat
 
 t = 1.0
 W = 10
@@ -27,15 +26,15 @@ L = 30
 
 for i in xrange(L):
     for j in xrange(W):
-        sys[(i, j)] = 4 * t
+        sys[lat(i, j)] = 4 * t
 
         # hoppig in y-direction
-        if j > 0 :
-            sys[(i, j), (i, j-1)] = - t
+        if j > 0:
+            sys[lat(i, j), lat(i, j - 1)] = -t
 
         #hopping in x-direction
         if i > 0:
-            sys[(i, j), (i-1, j)] = -t
+            sys[lat(i, j), lat(i - 1, j)] = -t
 
 # Then, define the leads:
 
@@ -45,29 +44,27 @@ for i in xrange(L):
 # realspace vector)
 sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
 lead0 = kwant.Builder(sym_lead0)
-lead0.default_site_group = lat
 
 for j in xrange(W):
-    lead0[(0, j)] = 4 * t
+    lead0[lat(0, j)] = 4 * t
 
     if j > 0:
-        lead0[(0, j), (0, j-1)] = - t
+        lead0[lat(0, j), lat(0, j - 1)] = -t
 
-    lead0[(1, j), (0, j)] = - t
+    lead0[lat(1, j), lat(0, j)] = -t
 
 # Then the lead to the right
 
 sym_lead1 = kwant.TranslationalSymmetry([lat.vec((1, 0))])
 lead1 = kwant.Builder(sym_lead1)
-lead1.default_site_group = lat
 
 for j in xrange(W):
-    lead1[(0, j)] = 4 * t
+    lead1[lat(0, j)] = 4 * t
 
     if j > 0:
-        lead1[(0, j), (0, j-1)] = - t
+        lead1[lat(0, j), lat(0, j - 1)] = -t
 
-    lead1[(1, j), (0, j)] = - t
+    lead1[lat(1, j), lat(0, j)] = -t
 
 # Then attach the leads to the system
 
@@ -80,7 +77,7 @@ kwant.plot(sys)
 
 # Finalize the system
 
-fsys = sys.finalized()
+sys = sys.finalized()
 
 # Now that we have the system, we can compute conductance
 
@@ -90,7 +87,7 @@ for ie in xrange(100):
     energy = ie * 0.01
 
     # compute the scattering matrix at energy energy
-    smatrix = kwant.solvers.sparse.solve(fsys, energy)
+    smatrix = kwant.solve(sys, energy)
 
     # compute the transmission probability from lead 0 to
     # lead 1
@@ -99,9 +96,10 @@ for ie in xrange(100):
 
 # Use matplotlib to write output
 # We should see conductance steps
-import pylab
+from matplotlib import pyplot
 
-pylab.plot(energies, data)
-pylab.xlabel("energy [in units of t]")
-pylab.ylabel("conductance [in units of e^2/h]")
-pylab.show()
+pyplot.figure()
+pyplot.plot(energies, data)
+pyplot.xlabel("energy [in units of t]")
+pyplot.ylabel("conductance [in units of e^2/h]")
+pyplot.show()
diff --git a/examples/tutorial1b.py b/tutorial/1-quantum_wire_revisited.py
similarity index 75%
rename from examples/tutorial1b.py
rename to tutorial/1-quantum_wire_revisited.py
index 0b8ee7fe48a6d316684ca25e8e7c5de39b8373cf..8f3202a5b4c098d887dd9756ef6fd77935f6dec6 100644
--- a/examples/tutorial1b.py
+++ b/tutorial/1-quantum_wire_revisited.py
@@ -13,7 +13,8 @@
 import kwant
 
 # For plotting
-import pylab
+from matplotlib import pyplot
+
 
 def make_system(a=1, t=1.0, W=10, L=30):
     # Start with an empty tight-binding system and a single square lattice.
@@ -21,10 +22,9 @@ def make_system(a=1, t=1.0, W=10, L=30):
     lat = kwant.lattice.Square(a)
 
     sys = kwant.Builder()
-    sys.default_site_group = lat
 
     #### Define the scattering region. ####
-    sys[((x, y) for x in range(L) for y in range(W))] = 4 * t
+    sys[(lat(x, y) for x in range(L) for y in range(W))] = 4 * t
     for hopping in lat.nearest:
         sys[sys.possible_hoppings(*hopping)] = -t
 
@@ -34,11 +34,10 @@ def make_system(a=1, t=1.0, W=10, L=30):
     # realspace vector)
     sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead0 = kwant.Builder(sym_lead0)
-    lead0.default_site_group = lat
 
-    lead0[((0, j) for j in xrange(W))] = 4 * t
+    lead0[(lat(0, j) for j in xrange(W))] = 4 * t
     for hopping in lat.nearest:
-        lead0[lead0.possible_hoppings(*hopping)] = - t
+        lead0[lead0.possible_hoppings(*hopping)] = -t
 
     # ... then the lead to the right.  We use a method that returns a copy of
     # `lead0` with its direction reversed.
@@ -50,17 +49,19 @@ def make_system(a=1, t=1.0, W=10, L=30):
 
     return sys
 
-def plot_conductance(fsys, energies):
+
+def plot_conductance(sys, energies):
     # Compute conductance
     data = []
     for energy in energies:
-        smatrix = kwant.solvers.sparse.solve(fsys, energy)
+        smatrix = kwant.solve(sys, energy)
         data.append(smatrix.transmission(1, 0))
 
-    pylab.plot(energies, data)
-    pylab.xlabel("energy [in units of t]")
-    pylab.ylabel("conductance [in units of e^2/h]")
-    pylab.show()
+    pyplot.figure()
+    pyplot.plot(energies, data)
+    pyplot.xlabel("energy [in units of t]")
+    pyplot.ylabel("conductance [in units of e^2/h]")
+    pyplot.show()
 
 
 def main():
@@ -70,10 +71,10 @@ def main():
     kwant.plot(sys)
 
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
     # We should see conductance steps.
-    plot_conductance(fsys, energies=[0.01 * i for i in xrange(100)])
+    plot_conductance(sys, energies=[0.01 * i for i in xrange(100)])
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/examples/tutorial2c.py b/tutorial/2-ab_ring.py
similarity index 76%
rename from examples/tutorial2c.py
rename to tutorial/2-ab_ring.py
index 06a5d6be1e009fee145de1b2d9964148078fef1e..16559d6ab362b0ff93769cac74866ca86e6f2ddd 100644
--- a/examples/tutorial2c.py
+++ b/tutorial/2-ab_ring.py
@@ -15,7 +15,8 @@ from math import pi
 import kwant
 
 # For plotting
-import pylab
+from matplotlib import pyplot
+
 
 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.
@@ -29,13 +30,13 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
     # Now, we aim for a more complex shape, namely a ring (or annulus)
     def ring(pos):
         (x, y) = pos
-        rsq = x**2 + y**2
-        return ( r1**2 < rsq < r2**2)
+        rsq = x ** 2 + y ** 2
+        return (r1 ** 2 < rsq < r2 ** 2)
 
     # and add the corresponding lattice points using the `shape`-function
-    sys[lat.shape(ring, (0, r1+1))] = 4 * t
+    sys[lat.shape(ring, (0, r1 + 1))] = 4 * t
     for hopping in lat.nearest:
-        sys[sys.possible_hoppings(*hopping)] = - t
+        sys[sys.possible_hoppings(*hopping)] = -t
 
     # 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
@@ -51,10 +52,10 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
 
         # possible_hoppings with the argument (1, 0) below
         # returns hoppings ordered as ((i+1, j), (i, j))
-        return iy0 < 0 and ix0 == 1 # ix1 == 0 then implied
+        return iy0 < 0 and ix0 == 1  # ix1 == 0 then implied
 
     # Modify only those hopings in x-direction that cross the branch cut
-    sys[(hop for hop in sys.possible_hoppings((1,0), lat, lat)
+    sys[(hop for hop in sys.possible_hoppings((1, 0), lat, lat)
          if crosses_branchcut(hop))] = fluxphase
 
     #### Define the leads. ####
@@ -64,11 +65,11 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
 
     def lead_shape(pos):
         (x, y) = pos
-        return (-1 < x < 1) and (-W/2 < y < W/2)
+        return (-1 < x < 1) and (-W / 2 < y < W / 2)
 
     lead0[lat.shape(lead_shape, (0, 0))] = 4 * t
     for hopping in lat.nearest:
-        lead0[lead0.possible_hoppings(*hopping)] = - t
+        lead0[lead0.possible_hoppings(*hopping)] = -t
 
     # Then the lead to the right
     # [again, obtained using reversed()]
@@ -81,23 +82,24 @@ def make_system(a=1, t=1.0, W=10, r1=10, r2=20):
     return sys
 
 
-def plot_conductance(fsys, energy, fluxes):
+def plot_conductance(sys, energy, fluxes):
     # compute conductance
     # global variable phi controls the flux
     global phi
 
-    normalized_fluxes = [flux/(2 * pi) for flux in fluxes]
+    normalized_fluxes = [flux / (2 * pi) for flux in fluxes]
     data = []
     for flux in fluxes:
         phi = flux
 
-        smatrix = kwant.solve(fsys, energy)
+        smatrix = kwant.solve(sys, energy)
         data.append(smatrix.transmission(1, 0))
 
-    pylab.plot(normalized_fluxes, data)
-    pylab.xlabel("flux [in units of the flux quantum]")
-    pylab.ylabel("conductance [in units of e^2/h]")
-    pylab.show()
+    pyplot.figure()
+    pyplot.plot(normalized_fluxes, data)
+    pyplot.xlabel("flux [in units of the flux quantum]")
+    pyplot.ylabel("conductance [in units of e^2/h]")
+    pyplot.show()
 
 
 def main():
@@ -107,10 +109,10 @@ def main():
     kwant.plot(sys)
 
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
     # We should see a conductance that is periodic with the flux quantum
-    plot_conductance(fsys, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
+    plot_conductance(sys, energy=0.15, fluxes=[0.01 * i * 3 * 2 * pi
                                                 for i in xrange(100)])
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/examples/tutorial2b.py b/tutorial/2-quantum_well.py
similarity index 80%
rename from examples/tutorial2b.py
rename to tutorial/2-quantum_well.py
index 94b778b8117bf29c73a3761f9aa95b8e74f572e2..d51ba0077dc47d0e5d3921caa895808ef88bf80b 100644
--- a/examples/tutorial2b.py
+++ b/tutorial/2-quantum_well.py
@@ -9,7 +9,7 @@
 import kwant
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
 # global variable governing the behavior of potential() in
 # make_system()
@@ -21,7 +21,6 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
     lat = kwant.lattice.Square(a)
 
     sys = kwant.Builder()
-    sys.default_site_group = lat
 
     #### Define the scattering region. ####
     # Potential profile
@@ -36,7 +35,7 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
     def onsite(site):
         return 4 * t + potential(site)
 
-    sys[((x, y) for x in range(L) for y in range(W))] = onsite
+    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
 
@@ -46,11 +45,10 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
     # realspace vector)
     sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead0 = kwant.Builder(sym_lead0)
-    lead0.default_site_group = lat
 
-    lead0[((0, j) for j in xrange(W))] = 4 * t
+    lead0[(lat(0, j) for j in xrange(W))] = 4 * t
     for hopping in lat.nearest:
-        lead0[lead0.possible_hoppings(*hopping)] = - t
+        lead0[lead0.possible_hoppings(*hopping)] = -t
 
     # ... then the lead to the right.  We use a method that returns a copy of
     # `lead0` with its direction reversed.
@@ -62,7 +60,8 @@ def make_system(a=1, t=1.0, W=10, L=30, L_well=10):
 
     return sys
 
-def plot_conductance(fsys, energy, welldepths):
+
+def plot_conductance(sys, energy, welldepths):
     # We specify that we want to not only read, but also write to a
     # global variable.
     global pot
@@ -73,13 +72,14 @@ def plot_conductance(fsys, energy, welldepths):
         # Set the global variable that defines the potential well depth
         pot = -welldepth
 
-        smatrix = kwant.solve(fsys, energy)
+        smatrix = kwant.solve(sys, energy)
         data.append(smatrix.transmission(1, 0))
 
-    pylab.plot(welldepths, data)
-    pylab.xlabel("well depth [in units of t]")
-    pylab.ylabel("conductance [in units of e^2/h]")
-    pylab.show()
+    pyplot.figure()
+    pyplot.plot(welldepths, data)
+    pyplot.xlabel("well depth [in units of t]")
+    pyplot.ylabel("conductance [in units of e^2/h]")
+    pyplot.show()
 
 
 def main():
@@ -89,10 +89,10 @@ def main():
     kwant.plot(sys)
 
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
     # We should see conductance steps.
-    plot_conductance(fsys, energy=0.2,
+    plot_conductance(sys, energy=0.2,
                      welldepths=[0.01 * i for i in xrange(100)])
 
 
diff --git a/examples/tutorial2a.py b/tutorial/2-spin_orbit.py
similarity index 71%
rename from examples/tutorial2a.py
rename to tutorial/2-spin_orbit.py
index d88a1746374233753d31f794055303446f74707f..009132b21f00674856445203deaa7ad692a6397e 100644
--- a/examples/tutorial2a.py
+++ b/tutorial/2-spin_orbit.py
@@ -13,7 +13,7 @@
 import kwant
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
 # For matrix support
 import numpy
@@ -31,30 +31,28 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
     lat = kwant.lattice.Square(a)
 
     sys = kwant.Builder()
-    sys.default_site_group = lat
 
     #### Define the scattering region. ####
-    sys[((x, y) for x in range(L) for y in range(W))] = 4 * t * sigma_0 + \
+    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
-    sys[sys.possible_hoppings((1, 0), lat, lat)] = - t * sigma_0 - \
+    sys[sys.possible_hoppings((1, 0), lat, lat)] = -t * sigma_0 - \
         1j * alpha * sigma_y
     # hoppings in y-directions
-    sys[sys.possible_hoppings((0, 1), lat, lat)] = - t * sigma_0 + \
+    sys[sys.possible_hoppings((0, 1), lat, lat)] = -t * sigma_0 + \
         1j * alpha * sigma_x
 
     #### Define the leads. ####
     # left lead
     sym_lead0 = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead0 = kwant.Builder(sym_lead0)
-    lead0.default_site_group = lat
 
-    lead0[((0, j) for j in xrange(W))] = 4 * t * sigma_0 + e_z * sigma_z
+    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 - \
+    lead0[lead0.possible_hoppings((1, 0), lat, lat)] = -t * sigma_0 - \
         1j * alpha * sigma_y
     # hoppings in y-directions
-    lead0[lead0.possible_hoppings((0, 1), lat, lat)] = - t * sigma_0 + \
+    lead0[lead0.possible_hoppings((0, 1), lat, lat)] = -t * sigma_0 + \
         1j * alpha * sigma_x
 
     # Then the lead to the right
@@ -67,17 +65,18 @@ def make_system(a=1, t=1.0, alpha=0.5, e_z=0.08, W=10, L=30):
 
     return sys
 
-def plot_conductance(fsys, energies):
+def plot_conductance(sys, energies):
     # Compute conductance
     data = []
     for energy in energies:
-        smatrix = kwant.solve(fsys, energy)
+        smatrix = kwant.solve(sys, energy)
         data.append(smatrix.transmission(1, 0))
 
-    pylab.plot(energies, data)
-    pylab.xlabel("energy [in units of t]")
-    pylab.ylabel("conductance [in units of e^2/h]")
-    pylab.show()
+    pyplot.figure()
+    pyplot.plot(energies, data)
+    pyplot.xlabel("energy [in units of t]")
+    pyplot.ylabel("conductance [in units of e^2/h]")
+    pyplot.show()
 
 
 def main():
@@ -87,10 +86,10 @@ def main():
     kwant.plot(sys)
 
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
     # We should see non-monotonic conductance steps.
-    plot_conductance(fsys, energies=[0.01 * i - 0.3 for i in xrange(100)])
+    plot_conductance(sys, energies=[0.01 * i - 0.3 for i in xrange(100)])
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/examples/tutorial3a.py b/tutorial/3-band_structure.py
similarity index 64%
rename from examples/tutorial3a.py
rename to tutorial/3-band_structure.py
index 3b9e4dcb31b572ebc2a42823cc095c6dc1152456..8a17633c104dd21b1d94ef8f802111b066369675 100644
--- a/examples/tutorial3a.py
+++ b/tutorial/3-band_structure.py
@@ -8,11 +8,10 @@
 
 import kwant
 
-import numpy as np
 from math import pi
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
 
 def make_lead(a=1, t=1.0, W=10):
@@ -21,39 +20,39 @@ def make_lead(a=1, t=1.0, W=10):
 
     sym_lead = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead = kwant.Builder(sym_lead)
-    lead.default_site_group = lat
 
     # build up one unit cell of the lead, and add the hoppings
     # to the next unit cell
     for j in xrange(W):
-        lead[(0, j)] = 4 * t
+        lead[lat(0, j)] = 4 * t
 
         if j > 0:
-            lead[(0, j), (0, j-1)] = - t
+            lead[lat(0, j), lat(0, j - 1)] = -t
 
-        lead[(1, j), (0, j)] = - t
+        lead[lat(1, j), lat(0, j)] = -t
 
     return lead
 
 
-def plot_bandstructure(flead, momenta):
+def plot_bandstructure(lead, momenta):
     # Use the method ``energies`` of the finalized lead to compute
     # the bandstructure
-    energy_list = [flead.energies(k) for k in momenta]
+    energy_list = [lead.energies(k) for k in momenta]
 
-    pylab.plot(momenta, energy_list)
-    pylab.xlabel("momentum [in units of (lattice constant)^-1]")
-    pylab.ylabel("energy [in units of t]")
-    pylab.show()
+    pyplot.figure()
+    pyplot.plot(momenta, energy_list)
+    pyplot.xlabel("momentum [in units of (lattice constant)^-1]")
+    pyplot.ylabel("energy [in units of t]")
+    pyplot.show()
 
 
 def main():
-    flead = make_lead().finalized()
+    lead = make_lead().finalized()
 
     # list of momenta at which the bands should be computed
-    momenta = np.arange(-pi, pi + .01, 0.02 * pi)
+    momenta = [-pi + 0.02 * pi * i for i in xrange(101)]
 
-    plot_bandstructure(flead, momenta)
+    plot_bandstructure(lead, momenta)
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/examples/tutorial3b.py b/tutorial/3-closed_system.py
similarity index 81%
rename from examples/tutorial3b.py
rename to tutorial/3-closed_system.py
index 273c8d4b423d071000f58300e9fea41b67e46768..f09603abff97a85a50c96c883b69effbee9678e4 100644
--- a/examples/tutorial3b.py
+++ b/tutorial/3-closed_system.py
@@ -16,7 +16,8 @@ import kwant
 import scipy.linalg as la
 
 # For plotting
-import pylab
+from matplotlib import pyplot
+
 
 def make_system(a=1, t=1.0, r=10):
     # Start with an empty tight-binding system and a single square lattice.
@@ -29,8 +30,8 @@ def make_system(a=1, t=1.0, r=10):
     # Define the quantum dot
     def circle(pos):
         (x, y) = pos
-        rsq = x**2 + y**2
-        return rsq < r**2
+        rsq = x ** 2 + y ** 2
+        return rsq < r ** 2
 
     def hopx(site1, site2):
         # The magnetic field is controlled by the global variable B
@@ -41,13 +42,13 @@ def make_system(a=1, t=1.0, r=10):
     # hoppings in x-direction
     sys[sys.possible_hoppings((1, 0), lat, lat)] = hopx
     # hoppings in y-directions
-    sys[sys.possible_hoppings((0, 1), lat, lat)] = - t
+    sys[sys.possible_hoppings((0, 1), lat, lat)] = -t
 
     # It's a closed system for a change, so no leads
     return sys
 
 
-def plot_spectrum(fsys, Bfields):
+def plot_spectrum(sys, Bfields):
     # global variable B controls the magnetic field
     global B
 
@@ -61,17 +62,18 @@ def plot_spectrum(fsys, Bfields):
         B = Bfield
 
         # Obtain the Hamiltonian as a dense matrix
-        ham_mat = fsys.hamiltonian_submatrix()[0]
+        ham_mat = sys.hamiltonian_submatrix()[0]
 
         ev = la.eigh(ham_mat, eigvals_only=True)
 
         # we only plot the 15 lowest eigenvalues
         energies.append(ev[:15])
 
-    pylab.plot(Bfields, energies)
-    pylab.xlabel("magnetic field [some arbitrary units]")
-    pylab.ylabel("energy [in units of t]")
-    pylab.show()
+    pyplot.figure()
+    pyplot.plot(Bfields, energies)
+    pyplot.xlabel("magnetic field [some arbitrary units]")
+    pyplot.ylabel("energy [in units of t]")
+    pyplot.show()
 
 
 def main():
@@ -81,11 +83,11 @@ def main():
     kwant.plot(sys)
 
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
     # We should observe energy levels that flow towards Landau
     # level energies with increasing magnetic field
-    plot_spectrum(fsys, [iB * 0.002 for iB in xrange(100)])
+    plot_spectrum(sys, [iB * 0.002 for iB in xrange(100)])
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/examples/tutorial4.py b/tutorial/4-graphene.py
similarity index 71%
rename from examples/tutorial4.py
rename to tutorial/4-graphene.py
index a31f00dde0a5965b07619fc5a4ee4da6818ec391..95782a819d213beec71f6f9ef89fc02cf28143b6 100644
--- a/examples/tutorial4.py
+++ b/tutorial/4-graphene.py
@@ -7,9 +7,8 @@
 #  - Application of all the aspects of tutorials 1-3 to a more complicated
 #    lattice, namely graphene
 
-from __future__ import division # so that 1/2 == 0.5, and not 0
-from math import pi, sqrt
-import numpy as np
+from __future__ import division  # so that 1/2 == 0.5, and not 0
+from math import pi, sqrt, tanh
 
 import kwant
 
@@ -17,13 +16,13 @@ import kwant
 import scipy.sparse.linalg as sla
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
 
 # Define the graphene lattice
-sin_30, cos_30 = (1/2, np.sqrt(3)/2)
+sin_30, cos_30 = (1 / 2, sqrt(3) / 2)
 graphene = kwant.make_lattice([(1, 0), (sin_30, cos_30)],
-                              [(0, 0), (0, 1/np.sqrt(3))])
+                              [(0, 0), (0, 1 / sqrt(3))])
 a, b = graphene.sublattices
 
 
@@ -33,27 +32,27 @@ def make_system(r=10, w=2.0, pot=0.1):
     # circular scattering region
     def circle(pos):
         x, y = pos
-        return x**2 + y**2 < r**2
+        return x ** 2 + y ** 2 < r ** 2
 
-    sys= kwant.Builder()
+    sys = kwant.Builder()
 
     # w: width and pot: potential maximum of the p-n junction
     def potential(site):
         (x, y) = site.pos
         d = y * cos_30 + x * sin_30
-        return pot * np.tanh(d / w)
+        return pot * tanh(d / w)
 
-    sys[graphene.shape(circle, (0,0))] = potential
+    sys[graphene.shape(circle, (0, 0))] = potential
 
     # specify the hoppings of the graphene lattice in the
     # format expected by possibe_hoppings()
     hoppings = (((0, 0), b, a), ((0, 1), b, a), ((-1, 1), b, a))
     for hopping in hoppings:
-        sys[sys.possible_hoppings(*hopping)] = - 1
+        sys[sys.possible_hoppings(*hopping)] = -1
 
     # Modify the scattering region
-    del sys[a(0,0)]
-    sys[a(-2,1), b(2, 2)] = -1
+    del sys[a(0, 0)]
+    sys[a(-2, 1), b(2, 2)] = -1
 
     #### Define the leads. ####
     # left lead
@@ -64,9 +63,9 @@ def make_system(r=10, w=2.0, pot=0.1):
         return (-1 < x < 1) and (-0.4 * r < y < 0.4 * r)
 
     lead0 = kwant.Builder(sym0)
-    lead0[graphene.shape(lead0_shape, (0,0))] = - pot
+    lead0[graphene.shape(lead0_shape, (0, 0))] = -pot
     for hopping in hoppings:
-        lead0[lead0.possible_hoppings(*hopping)] = - 1
+        lead0[lead0.possible_hoppings(*hopping)] = -1
 
     # The second lead, going to the top right
     sym1 = kwant.TranslationalSymmetry([graphene.vec((0, 1))])
@@ -78,9 +77,9 @@ def make_system(r=10, w=2.0, pot=0.1):
         return (-1 < u < 1) and (-0.4 * r < v < 0.4 * r)
 
     lead1 = kwant.Builder(sym1)
-    lead1[graphene.shape(lead1_shape, (0,0))] = pot
+    lead1[graphene.shape(lead1_shape, (0, 0))] = pot
     for hopping in hoppings:
-        lead1[lead1.possible_hoppings(*hopping)] = - 1
+        lead1[lead1.possible_hoppings(*hopping)] = -1
 
     return sys, [lead0, lead1]
 
@@ -93,23 +92,24 @@ def compute_evs(sys):
         # This requires scipy version >= 0.9.0
         # Failure (i.e. insufficient scipy version) is not critical
         # for the remainder of the tutorial, hence the try-block
-        evs = scipy.sparse.linalg.eigs(sparse_mat, 2)[0]
+        evs = sla.eigs(sparse_mat, 2)[0]
         print evs
     except:
         pass
 
 
-def plot_conductance(fsys, energies):
+def plot_conductance(sys, energies):
     # Compute transmission as a function of energy
     data = []
     for energy in energies:
-        smatrix = kwant.solve(fsys, energy)
+        smatrix = kwant.solve(sys, energy)
         data.append(smatrix.transmission(0, 1))
 
-    pylab.plot(energies, data)
-    pylab.xlabel("energy [in units of t]")
-    pylab.ylabel("conductance [in units of e^2/h]")
-    pylab.show()
+    pyplot.figure()
+    pyplot.plot(energies, data)
+    pyplot.xlabel("energy [in units of t]")
+    pyplot.ylabel("conductance [in units of e^2/h]")
+    pyplot.show()
 
 
 def plot_bandstructure(flead, momenta):
@@ -117,10 +117,11 @@ def plot_bandstructure(flead, momenta):
     # the bandstructure
     energy_list = [flead.energies(k) for k in momenta]
 
-    pylab.plot(momenta, energy_list)
-    pylab.xlabel("momentum [in untis of (lattice constant)^-1]")
-    pylab.ylabel("energy [in units of t]")
-    pylab.show()
+    pyplot.figure()
+    pyplot.plot(momenta, energy_list)
+    pyplot.xlabel("momentum [in untis of (lattice constant)^-1]")
+    pyplot.ylabel("energy [in units of t]")
+    pyplot.show()
 
 
 def main():
@@ -148,15 +149,15 @@ def main():
     kwant.plot(sys, symbols=plotter_symbols)
 
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
     # Compute the band structure of lead 0.
-    momenta = np.arange(-pi, pi + .01, 0.1 * pi)
-    plot_bandstructure(fsys.leads[0], momenta)
+    momenta = [-pi + 0.02 * pi * i for i in xrange(101)]
+    plot_bandstructure(sys.leads[0], momenta)
 
     # Plot conductance.
-    energies = np.arange(-2 * pot, 2 * pot, pot / 10.5)
-    plot_conductance(fsys, energies)
+    energies = [-2 * pot + 4. / 50. * pot * i for i in xrange(51)]
+    plot_conductance(sys, energies)
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/examples/tutorial5a.py b/tutorial/5-superconductor_band_structure.py
similarity index 64%
rename from examples/tutorial5a.py
rename to tutorial/5-superconductor_band_structure.py
index 60bf12ef9b413df2d6f312c8dba9b3e6fda300ea..27e39d3d298e9bf1cbed44b8cdbe940481160d42 100644
--- a/examples/tutorial5a.py
+++ b/tutorial/5-superconductor_band_structure.py
@@ -16,10 +16,11 @@ import numpy as np
 from math import pi
 
 # For plotting
-import pylab
+from matplotlib import pyplot
+
+tau_x = np.array([[0, 1], [1, 0]])
+tau_z = np.array([[1, 0], [0, -1]])
 
-tau_x = np.array([[0,1],[1,0]])
-tau_z = np.array([[1,0],[0,-1]])
 
 def make_lead(a=1, t=1.0, mu=0.7, Delta=0.1, W=10):
     # Start with an empty lead with a single square lattice
@@ -27,41 +28,41 @@ def make_lead(a=1, t=1.0, mu=0.7, Delta=0.1, W=10):
 
     sym_lead = kwant.TranslationalSymmetry([lat.vec((-1, 0))])
     lead = kwant.Builder(sym_lead)
-    lead.default_site_group = lat
 
     # build up one unit cell of the lead, and add the hoppings
     # to the next unit cell
     for j in xrange(W):
-        lead[(0, j)] = (4 * t - mu) * tau_z + Delta * tau_x
+        lead[lat(0, j)] = (4 * t - mu) * tau_z + Delta * tau_x
 
         if j > 0:
-            lead[(0, j), (0, j-1)] = - t * tau_z
+            lead[lat(0, j), lat(0, j - 1)] = -t * tau_z
 
-        lead[(1, j), (0, j)] = - t * tau_z
+        lead[lat(1, j), lat(0, j)] = -t * tau_z
 
     return lead
 
 
-def plot_bandstructure(flead, momenta):
+def plot_bandstructure(lead, momenta):
     # Use the method ``energies`` of the finalized lead to compute
     # the bandstructure
-    energy_list = [flead.energies(k) for k in momenta]
+    energy_list = [lead.energies(k) for k in momenta]
 
-    pylab.plot(momenta, energy_list)
-    pylab.xlabel("momentum [in untis of (lattice constant)^-1]")
-    pylab.ylabel("energy [in units of t]")
-    pylab.ylim([-0.8, 0.8])
-    pylab.show()
+    pyplot.figure()
+    pyplot.plot(momenta, energy_list)
+    pyplot.xlabel("momentum [in untis of (lattice constant)^-1]")
+    pyplot.ylabel("energy [in units of t]")
+    pyplot.ylim([-0.8, 0.8])
+    pyplot.show()
 
 
 def main():
     # Make system and finalize it right away.
-    flead = make_lead().finalized()
+    lead = make_lead().finalized()
 
     # list of momenta at which the bands should be computed
-    momenta = np.arange(-1.5, 1.5 + .0001, 0.002 * pi)
+    momenta = np.linspace(-1.5, 1.5, 201)
 
-    plot_bandstructure(flead, momenta)
+    plot_bandstructure(lead, momenta)
 
 
 # Call the main function if the script gets executed (as opposed to imported).
diff --git a/examples/tutorial5b.py b/tutorial/5-superconductor_transport.py
similarity index 80%
rename from examples/tutorial5b.py
rename to tutorial/5-superconductor_transport.py
index 2771599c90a68007e4e5d1008f19583cb8de755b..977fb1e121960ced4562f25725b788b1103cb77b 100644
--- a/examples/tutorial5b.py
+++ b/tutorial/5-superconductor_transport.py
@@ -10,10 +10,8 @@
 import kwant
 
 # For plotting
-import pylab
+from matplotlib import pyplot
 
-# For matrix support
-import numpy
 
 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):
@@ -35,14 +33,14 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
          for y in range(W))] = mu - 4 * t - barrier
 
     # hoppings in x and y-directions, for both electrons and holes
-    sys[sys.possible_hoppings((1, 0), lat_e, lat_e)] = - t
-    sys[sys.possible_hoppings((0, 1), lat_e, lat_e)] = - t
+    sys[sys.possible_hoppings((1, 0), lat_e, lat_e)] = -t
+    sys[sys.possible_hoppings((0, 1), lat_e, lat_e)] = -t
     sys[sys.possible_hoppings((1, 0), lat_h, lat_h)] = t
     sys[sys.possible_hoppings((0, 1), lat_h, lat_h)] = t
 
     # Superconducting order parameter enters as hopping between
     # electrons and holes
-    sys[((lat_e(x,y), lat_h(x, y)) for x in range(Deltapos, L)
+    sys[((lat_e(x, y), lat_h(x, y)) for x in range(Deltapos, L)
          for y in range(W))] = Delta
 
     #### Define the leads. ####
@@ -52,8 +50,8 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
 
     lead0[(lat_e(0, j) for j in xrange(W))] = 4 * t - mu
     # hoppings in x and y-direction
-    lead0[lead0.possible_hoppings((1, 0), lat_e, lat_e)] = - t
-    lead0[lead0.possible_hoppings((0, 1), lat_e, lat_e)] = - t
+    lead0[lead0.possible_hoppings((1, 0), lat_e, lat_e)] = -t
+    lead0[lead0.possible_hoppings((0, 1), lat_e, lat_e)] = -t
 
     # left hole lead
     sym_lead1 = kwant.TranslationalSymmetry([lat_h.vec((-1, 0))])
@@ -73,8 +71,8 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
     lead2[(lat_e(0, j) for j in xrange(W))] = 4 * t - mu
     lead2[(lat_h(0, j) for j in xrange(W))] = mu - 4 * t
     # hoppings in x and y-direction
-    lead2[lead2.possible_hoppings((1, 0), lat_e, lat_e)] = - t
-    lead2[lead2.possible_hoppings((0, 1), lat_e, lat_e)] = - t
+    lead2[lead2.possible_hoppings((1, 0), lat_e, lat_e)] = -t
+    lead2[lead2.possible_hoppings((0, 1), lat_e, lat_e)] = -t
     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
@@ -86,20 +84,22 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
 
     return sys
 
-def plot_conductance(fsys, energies):
+
+def plot_conductance(sys, energies):
     # Compute conductance
     data = []
     for energy in energies:
-        smatrix = kwant.solve(fsys, energy)
+        smatrix = kwant.solve(sys, energy)
         # Conductance is N - R_ee + R_he
         data.append(smatrix.submatrix(0, 0).shape[0] -
                     smatrix.transmission(0, 0) +
                     smatrix.transmission(1, 0))
 
-    pylab.plot(energies, data)
-    pylab.xlabel("energy [in units of t]")
-    pylab.ylabel("conductance [in units of e^2/h]")
-    pylab.show()
+    pyplot.figure()
+    pyplot.plot(energies, data)
+    pyplot.xlabel("energy [in units of t]")
+    pyplot.ylabel("conductance [in units of e^2/h]")
+    pyplot.show()
 
 
 def main():
@@ -109,9 +109,9 @@ def main():
     kwant.plot(sys)
 
     # Finalize the system.
-    fsys = sys.finalized()
+    sys = sys.finalized()
 
-    plot_conductance(fsys, energies=[0.002 * i for i in xrange(100)])
+    plot_conductance(sys, energies=[0.002 * i for i in xrange(100)])
 
 
 # Call the main function if the script gets executed (as opposed to imported).