diff --git a/.gitignore b/.gitignore
index 7357693774aea277987ca67c03ecd729581debb1..45e8fc9d83981c9d2121e252141e59f6d826b1f2 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,13 +8,14 @@
 /build
 /dist
 /doc/build
-/doc/source/tutorial/*.py
 /doc/source/reference/generated/
-/doc/source/images/*.png
-/doc/source/images/*.pdf
-/doc/source/images/.*_flag
-/doc/source/images/[a-zA-Z]*.py
-/doc/source/images/*.txt
+/doc/source/code/include/*.py
+/doc/source/code/figure/*.png
+/doc/source/code/figure/*.pdf
+/doc/source/code/figure/.*_flag
+/doc/source/code/figure/[a-zA-Z]*.py
+/doc/source/code/figure/*.txt
+/doc/source/code/download/
 /build.conf
 /kwant.egg-info/
 /MANIFEST.in
diff --git a/doc/Makefile b/doc/Makefile
index 8df4d59b00d3f81b876e8a3249cca7ddb5eab9c8..b9fd88395b2ca3a84cd4de24648bbd49ff8914c7 100644
--- a/doc/Makefile
+++ b/doc/Makefile
@@ -1,6 +1,6 @@
 # Makefile for Sphinx documentation
 
-# Copyright 2011-2013 Kwant authors.
+# Copyright 2011-2017 Kwant authors.
 #
 # This file is part of Kwant.  It is subject to the license terms in the file
 # LICENSE.rst found in the top-level directory of this distribution and at
@@ -21,21 +21,23 @@ ALLSPHINXOPTS   = -d $(BUILDDIR)/doctrees $(PAPEROPT_$(PAPER)) $(SPHINXOPTS) sou
 
 # We convert all SVG files to PDF for LaTeX output.  For HTML output, we don't
 # create PNGs but rather use the SVG files directly.
-IMAGESOURCES    = $(shell find source -name "*.svg")
-GENERATEDPDF    = $(patsubst %.svg,%.pdf,$(IMAGESOURCES))
+FIGURESOURCES    = $(shell find source -name "*.svg")
+GENERATEDPDF    = $(patsubst %.svg,%.pdf,$(FIGURESOURCES))
 
-# Image generation from patched tutorial scripts
+# Figure generation from patched tutorial scripts
 #
 # As make does not support the generation of multiple targets by a single
 # invocation of a (non-implicit) rule, we use a trick: We pretend to be
-# generating a single (empty) flag file per invocation.  The image files are
+# generating a single (empty) flag file per invocation.  The figure files are
 # generated as well, but only as side-effects.  Each flag file is used to
-# remember the time at which the corresponding image-generating script was run.
+# remember the time at which the corresponding figure-generating script was run.
 # This works perfectly unless the actual output files are deleted without
 # deleting the corresponding flag file.
-SCRIPTS = $(patsubst source/images/%.diff,%,$(wildcard source/images/*.py.diff))
-FLAGS = $(patsubst %.py,source/images/.%_flag,$(SCRIPTS))
-INCLUDES = $(patsubst %,source/tutorial/%,$(SCRIPTS))
+FIGSCRIPTS = $(patsubst %.diff,%,$(notdir $(wildcard source/code/figure/*.py.diff)))
+FIGURES = $(patsubst %.py,source/code/figure/.%_flag,$(FIGSCRIPTS))
+SCRIPTS = $(sort $(FIGSCRIPTS) $(notdir $(wildcard source/code/include/*.py)))
+INCLUDES = $(patsubst %,source/code/include/%,$(SCRIPTS))
+DOWNLOADS = $(patsubst %,source/code/download/%,$(SCRIPTS))
 
 .PHONY: help clean realclean html dirhtml pickle json htmlhelp qthelp latex changes linkcheck doctest
 
@@ -59,39 +61,40 @@ clean:
 	-rm -rf source/reference/generated
 
 realclean: clean
-	-rm -f $(FLAGS)
-	-rm -f $(INCLUDES)
-	-rm -f $(patsubst %,source/images/%,$(SCRIPTS))
-	-rm -f $(patsubst %.py,source/images/%_*.png,$(SCRIPTS))
-	-rm -f $(patsubst %.py,source/images/%_*.pdf,$(SCRIPTS))
-
-html:	$(FLAGS) $(INCLUDES)
+	-rm -f $(FIGURES)
+	-rm -f $(patsubst %,source/code/include/%,$(FIGSCRIPTS))
+	-rm -f $(DOWNLOADS)
+	-rm -f $(patsubst %,source/code/figure/%,$(FIGSCRIPTS))
+	-rm -f $(patsubst %.py,source/code/figure/%_*.png,$(FIGSCRIPTS))
+	-rm -f $(patsubst %.py,source/code/figure/%_*.pdf,$(FIGSCRIPTS))
+
+html:	$(FIGURES) $(INCLUDES) $(DOWNLOADS)
 	$(SPHINXBUILD) -b html $(ALLSPHINXOPTS) $(BUILDDIR)/html
 	@echo
 	@echo "Build finished. The HTML pages are in $(BUILDDIR)/html."
 
-dirhtml: $(FLAGS) $(INCLUDES)
+dirhtml: $(FIGURES) $(INCLUDES) $(DOWNLOADS)
 	$(SPHINXBUILD) -b dirhtml $(ALLSPHINXOPTS) $(BUILDDIR)/dirhtml
 	@echo
 	@echo "Build finished. The HTML pages are in $(BUILDDIR)/dirhtml."
 
-pickle: $(FLAGS) $(INCLUDES)
+pickle: $(FIGURES) $(INCLUDES) $(DOWNLOADS)
 	$(SPHINXBUILD) -b pickle $(ALLSPHINXOPTS) $(BUILDDIR)/pickle
 	@echo
 	@echo "Build finished; now you can process the pickle files."
 
-json:   $(FLAGS) $(INCLUDES)
+json:   $(FIGURES) $(INCLUDES) $(DOWNLOADS)
 	$(SPHINXBUILD) -b json $(ALLSPHINXOPTS) $(BUILDDIR)/json
 	@echo
 	@echo "Build finished; now you can process the JSON files."
 
-htmlhelp: $(FLAGS) $(INCLUDES)
+htmlhelp: $(FIGURES) $(INCLUDES) $(DOWNLOADS)
 	$(SPHINXBUILD) -b htmlhelp $(ALLSPHINXOPTS) $(BUILDDIR)/htmlhelp
 	@echo
 	@echo "Build finished; now you can run HTML Help Workshop with the" \
 	      ".hhp project file in $(BUILDDIR)/htmlhelp."
 
-qthelp: $(FLAGS) $(INCLUDES)
+qthelp: $(FIGURES) $(INCLUDES) $(DOWNLOADS)
 	$(SPHINXBUILD) -b qthelp $(ALLSPHINXOPTS) $(BUILDDIR)/qthelp
 	@echo
 	@echo "Build finished; now you can run "qcollectiongenerator" with the" \
@@ -100,7 +103,7 @@ qthelp: $(FLAGS) $(INCLUDES)
 	@echo "To view the help file:"
 	@echo "# assistant -collectionFile $(BUILDDIR)/qthelp/kwant.qhc"
 
-latex:  $(GENERATEDPDF) $(FLAGS) $(INCLUDES)
+latex:  $(GENERATEDPDF) $(FIGURES) $(INCLUDES) $(DOWNLOADS)
 	$(SPHINXBUILD) -b latex $(ALLSPHINXOPTS) $(BUILDDIR)/latex
 	@echo
 	@echo "Build finished; the LaTeX files are in $(BUILDDIR)/latex."
@@ -129,16 +132,20 @@ doctest:
 # Make tutorial scripts by extracting the (complete!) context of the "patches".
 # We make sure not to use 'wiggle' here.
 .SECONDARY:
-source/tutorial/%.py: source/images/%.py.diff
+source/code/include/%.py: source/code/figure/%.py.diff
 	@sed -n '/^[- ]/ s/^.//p' <$< >$@
 	@touch -r $< $@
 
-# Make the image generation scripts by patching tutorial scripts.  If the
+source/code/download/%.py: source/code/include/%.py
+	@mkdir -p source/code/download
+	@grep -v '^#HIDDEN' <$< >$@
+
+# Make the figure generation scripts by patching tutorial scripts.  If the
 # tutorial scripts haven't been modified, don't patch but directly extract the
-# image generation scripts.  This means that 'wiggle' is only needed when the
+# figure generation scripts.  This means that 'wiggle' is only needed when the
 # tutorial scripts have been modified.
 .SECONDARY:
-source/images/%.py: source/tutorial/%.py
+source/code/figure/%.py: source/code/include/%.py
 	@if [ $< -nt $@.diff ]; then \
 	    cp $< $@; \
 	    rm -f $@.porig; \
@@ -153,20 +160,20 @@ source/images/%.py: source/tutorial/%.py
 	    touch -r $@.diff $@; \
 	fi
 
-# Make the image generation scripts also depend on the diffs.
+# Make the figure generation scripts also depend on the diffs.
 define makedep
-source/images/$(1): source/images/$(1).diff
+source/code/figure/$(1): source/code/figure/$(1).diff
 endef
-$(foreach name,$(SCRIPTS),$(eval $(call makedep,$(name))))
+$(foreach name,$(FIGSCRIPTS),$(eval $(call makedep,$(name))))
 
-# Run an image generation script.  When successful, and if the script is newer
+# Run an figure generation script.  When successful, and if the script is newer
 # than the corresponding diff, recreate the latter.  Note that the
-# corresponding tutorial script cannot be newer, since if it is, the image
+# corresponding tutorial script cannot be newer, since if it is, the figure
 # generation script is generated from it by patching.
-source/images/.%_flag: source/images/%.py
+source/code/figure/.%_flag: source/code/figure/%.py
 	cd $(dir $<) && python3 $(notdir $<)
 	@if [ ! -f $<.diff -o $< -nt $<.diff ]; then \
-	    wiggle --diff --lines source/tutorial/$(notdir $<) $< >$<.diff; \
+	    wiggle --diff --lines source/code/include/$(notdir $<) $< >$<.diff; \
 	    touch -r $< $<.diff; \
 	fi
 	@rm -f $<.porig
diff --git a/doc/source/images/_defs.py b/doc/source/code/figure/_defs.py
similarity index 94%
rename from doc/source/images/_defs.py
rename to doc/source/code/figure/_defs.py
index 42c864b71ac6e66e0b387845fea83c5a91327b60..985e2e9a3cc8f0057e53ea9726ff252975ad7343 100644
--- a/doc/source/images/_defs.py
+++ b/doc/source/code/figure/_defs.py
@@ -9,7 +9,7 @@ matplotlib.use('Agg')
 ################################################################
 import sys
 from distutils.util import get_platform
-sys.path.insert(0, "../../../build/lib.{0}-{1}.{2}".format(
+sys.path.insert(0, "../../../../build/lib.{0}-{1}.{2}".format(
         get_platform(), *sys.version_info[:2]))
 
 ################################################################
diff --git a/doc/source/images/ab_ring.py.diff b/doc/source/code/figure/ab_ring.py.diff
similarity index 100%
rename from doc/source/images/ab_ring.py.diff
rename to doc/source/code/figure/ab_ring.py.diff
diff --git a/doc/source/images/ab_ring_sketch.svg b/doc/source/code/figure/ab_ring_sketch.svg
similarity index 100%
rename from doc/source/images/ab_ring_sketch.svg
rename to doc/source/code/figure/ab_ring_sketch.svg
diff --git a/doc/source/images/ab_ring_sketch2.svg b/doc/source/code/figure/ab_ring_sketch2.svg
similarity index 100%
rename from doc/source/images/ab_ring_sketch2.svg
rename to doc/source/code/figure/ab_ring_sketch2.svg
diff --git a/doc/source/images/band_structure.py.diff b/doc/source/code/figure/band_structure.py.diff
similarity index 100%
rename from doc/source/images/band_structure.py.diff
rename to doc/source/code/figure/band_structure.py.diff
diff --git a/doc/source/images/closed_system.py.diff b/doc/source/code/figure/closed_system.py.diff
similarity index 100%
rename from doc/source/images/closed_system.py.diff
rename to doc/source/code/figure/closed_system.py.diff
diff --git a/doc/source/images/discretize.py.diff b/doc/source/code/figure/discretize.py.diff
similarity index 100%
rename from doc/source/images/discretize.py.diff
rename to doc/source/code/figure/discretize.py.diff
diff --git a/doc/source/images/faq.py.diff b/doc/source/code/figure/faq.py.diff
similarity index 100%
rename from doc/source/images/faq.py.diff
rename to doc/source/code/figure/faq.py.diff
diff --git a/doc/source/images/graphene.py.diff b/doc/source/code/figure/graphene.py.diff
similarity index 100%
rename from doc/source/images/graphene.py.diff
rename to doc/source/code/figure/graphene.py.diff
diff --git a/doc/source/images/kernel_polynomial_method.py.diff b/doc/source/code/figure/kernel_polynomial_method.py.diff
similarity index 100%
rename from doc/source/images/kernel_polynomial_method.py.diff
rename to doc/source/code/figure/kernel_polynomial_method.py.diff
diff --git a/doc/source/images/magnetic_texture.py.diff b/doc/source/code/figure/magnetic_texture.py.diff
similarity index 100%
rename from doc/source/images/magnetic_texture.py.diff
rename to doc/source/code/figure/magnetic_texture.py.diff
diff --git a/doc/source/images/plot_graphene.py.diff b/doc/source/code/figure/plot_graphene.py.diff
similarity index 100%
rename from doc/source/images/plot_graphene.py.diff
rename to doc/source/code/figure/plot_graphene.py.diff
diff --git a/doc/source/images/plot_qahe.py.diff b/doc/source/code/figure/plot_qahe.py.diff
similarity index 100%
rename from doc/source/images/plot_qahe.py.diff
rename to doc/source/code/figure/plot_qahe.py.diff
diff --git a/doc/source/images/plot_zincblende.py.diff b/doc/source/code/figure/plot_zincblende.py.diff
similarity index 100%
rename from doc/source/images/plot_zincblende.py.diff
rename to doc/source/code/figure/plot_zincblende.py.diff
diff --git a/doc/source/images/quantum_well.py.diff b/doc/source/code/figure/quantum_well.py.diff
similarity index 100%
rename from doc/source/images/quantum_well.py.diff
rename to doc/source/code/figure/quantum_well.py.diff
diff --git a/doc/source/images/quantum_wire.py.diff b/doc/source/code/figure/quantum_wire.py.diff
similarity index 100%
rename from doc/source/images/quantum_wire.py.diff
rename to doc/source/code/figure/quantum_wire.py.diff
diff --git a/doc/source/images/spin_orbit.py.diff b/doc/source/code/figure/spin_orbit.py.diff
similarity index 100%
rename from doc/source/images/spin_orbit.py.diff
rename to doc/source/code/figure/spin_orbit.py.diff
diff --git a/doc/source/images/superconductor.py.diff b/doc/source/code/figure/superconductor.py.diff
similarity index 100%
rename from doc/source/images/superconductor.py.diff
rename to doc/source/code/figure/superconductor.py.diff
diff --git a/doc/source/images/superconductor_transport_sketch.svg b/doc/source/code/figure/superconductor_transport_sketch.svg
similarity index 100%
rename from doc/source/images/superconductor_transport_sketch.svg
rename to doc/source/code/figure/superconductor_transport_sketch.svg
diff --git a/doc/source/tutorial/quantum_wire_revisited.py b/doc/source/code/include/quantum_wire_revisited.py
similarity index 100%
rename from doc/source/tutorial/quantum_wire_revisited.py
rename to doc/source/code/include/quantum_wire_revisited.py
diff --git a/doc/source/pre/whatsnew/1.3.rst b/doc/source/pre/whatsnew/1.3.rst
index 7380e4e14a26a096a8cb8a7b79a3e704ce8d6c68..0f4e250fdcf608271da94544e6396cddd4638ea8 100644
--- a/doc/source/pre/whatsnew/1.3.rst
+++ b/doc/source/pre/whatsnew/1.3.rst
@@ -28,11 +28,11 @@ tight-binding band structures or construct systems with different/lower
 symmetry. For example in just a few lines we can construct a two-band model that exhibits
 a quantum anomalous spin Hall phase:
 
-.. literalinclude:: ../../tutorial/plot_qahe.py
+.. literalinclude:: ../../code/include/plot_qahe.py
     :start-after: HIDDEN_BEGIN_model
     :end-before: HIDDEN_END_model
 
-From: :download:`QAHE example script <../../tutorial/plot_qahe.py>`
+From: :download:`QAHE example script <../../code/download/plot_qahe.py>`
 
 See the tutorial: :doc:`../../tutorial/discretize`
 
@@ -71,13 +71,13 @@ The example below shows edge states of a quantum anomalous Hall phase
 of the two-band model shown in the `above section
 <#tools-for-continuum-hamiltonians>`_:
 
-.. literalinclude:: ../../tutorial/plot_qahe.py
+.. literalinclude:: ../../code/include/plot_qahe.py
     :start-after: HIDDEN_BEGIN_current
     :end-before: HIDDEN_END_current
 
-.. image:: ../../images/plot_qahe_current.*
+.. image:: ../../code/figure/plot_qahe_current.*
 
-From: :download:`QAHE example script <../../tutorial/plot_qahe.py>`
+From: :download:`QAHE example script <../../code/download/plot_qahe.py>`
 
 Scattering states with discrete symmetries and conservation laws
 ----------------------------------------------------------------
diff --git a/doc/source/tutorial/discretize.rst b/doc/source/tutorial/discretize.rst
index a36149b62defc48cc2a36aee9fea80172e3e7265..c4a5b6a864db1843470365510ee8c73d44fe5d78 100644
--- a/doc/source/tutorial/discretize.rst
+++ b/doc/source/tutorial/discretize.rst
@@ -18,7 +18,7 @@ continuum models and for discretizing them into tight-binding models.
 
 .. seealso::
     The complete source code of this tutorial can be found in
-    :download:`tutorial/discretize.py <../../../tutorial/discretize.py>`
+    :download:`discretize.py </code/download/discretize.py>`
 
 
 .. _tutorial_discretizer_introduction:
@@ -60,7 +60,7 @@ symmetry that serves as a template.
 (We will see how to use the template to build systems with a particular
 shape later).
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_symbolic_discretization
     :end-before: #HIDDEN_END_symbolic_discretization
 
@@ -70,7 +70,7 @@ discretization process.
 
 The builder produced by ``discretize`` may be printed to show the source code of its onsite and hopping functions (this is a special feature of builders returned by ``discretize``):
 
-.. literalinclude:: ../images/discretizer_intro_verbose.txt
+.. literalinclude:: /code/figure/discretizer_intro_verbose.txt
 
 .. specialnote:: Technical details
 
@@ -134,7 +134,7 @@ where :math:`V(x, y)` is some arbitrary potential.
 First, use ``discretize`` to obtain a
 builder that we will use as a template:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_template
     :end-before: #HIDDEN_END_template
 
@@ -142,18 +142,18 @@ We now use this system with the `~kwant.builder.Builder.fill`
 method of `~kwant.builder.Builder` to construct the system we
 want to investigate:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_fill
     :end-before: #HIDDEN_END_fill
 
 After finalizing this system, we can plot one of the system's
 energy eigenstates:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_plot_eigenstate
     :end-before: #HIDDEN_END_plot_eigenstate
 
-.. image:: ../images/discretizer_gs.*
+.. image:: /code/figure/discretizer_gs.*
 
 Note in the above that we provided the function ``V`` to
 ``syst.hamiltonian_submatrix`` using ``params=dict(V=potential)``, rather than
@@ -171,32 +171,32 @@ model [1]_ [2]_, one can provide matrix input to `~kwant.continuum.discretize`
 using ``identity`` and ``kron``. For example, the definition of the BHZ model can be
 written succinctly as:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_define_qsh
     :end-before: #HIDDEN_END_define_qsh
 
 We can then make a ribbon out of this template system:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_define_qsh_build
     :end-before: #HIDDEN_END_define_qsh_build
 
 and plot its dispersion using `kwant.plotter.bands`:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_plot_qsh_band
     :end-before: #HIDDEN_END_plot_qsh_band
 
-.. image:: ../images/discretizer_qsh_band.*
+.. image:: /code/figure/discretizer_qsh_band.*
 
 In the above we see the edge states of the quantum spin Hall effect, which
 we can visualize using `kwant.plotter.map`:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_plot_qsh_wf
     :end-before: #HIDDEN_END_plot_qsh_wf
 
-.. image:: ../images/discretizer_qsh_wf.*
+.. image:: /code/figure/discretizer_qsh_wf.*
 
 
 Limitations of discretization
@@ -233,28 +233,28 @@ example. Let us start from the continuum Hamiltonian
 We start by defining this model as a string and setting the value of the
 :math:`α` parameter:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_ls_def
     :end-before: #HIDDEN_END_ls_def
 
 Now we can use `kwant.continuum.lambdify` to obtain a function that computes
 :math:`H(k)`:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_ls_hk_cont
     :end-before: #HIDDEN_END_ls_hk_cont
 
 We can also construct a discretized approximation using
 `kwant.continuum.discretize`, in a similar manner to previous examples:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_ls_hk_tb
     :end-before: #HIDDEN_END_ls_hk_tb
 
 Below we can see the continuum and tight-binding dispersions for two
 different values of the discretization grid spacing :math:`a`:
 
-.. image:: ../images/discretizer_lattice_spacing.*
+.. image:: /code/figure/discretizer_lattice_spacing.*
 
 
 We clearly see that the smaller grid spacing is, the better we approximate
@@ -279,20 +279,20 @@ It is possible to use ``identity`` (for identity matrix), ``kron`` (for Kronecke
 expressions involving matrices. Matrices can also be provided explicitly using
 square ``[]`` brackets. For example, all following expressions are equivalent:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_subs_1
     :end-before: #HIDDEN_END_subs_1
 
-.. literalinclude:: ../images/discretizer_subs_1.txt
+.. literalinclude:: /code/figure/discretizer_subs_1.txt
 
 We can use the ``locals`` keyword parameter to substitute expressions
 and numerical values:
 
-.. literalinclude:: discretize.py
+.. literalinclude:: /code/include/discretize.py
     :start-after: #HIDDEN_BEGIN_subs_2
     :end-before: #HIDDEN_END_subs_2
 
-.. literalinclude:: ../images/discretizer_subs_2.txt
+.. literalinclude:: /code/figure/discretizer_subs_2.txt
 
 Symbolic expressions obtained in this way can be directly passed to all
 ``discretizer`` functions.
diff --git a/doc/source/tutorial/faq.rst b/doc/source/tutorial/faq.rst
index eccec8a8993eead562c9c368e85d2a230a071ab8..3833f4743bdd9450f051a426c3ad745ee9005fa1 100644
--- a/doc/source/tutorial/faq.rst
+++ b/doc/source/tutorial/faq.rst
@@ -48,11 +48,11 @@ combination of family and tag uniquely defines a site.
 
 For example let us create an empty tight binding system and add two sites:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_site
     :end-before: #HIDDEN_END_site
 
-.. image:: ../images/faq_site.*
+.. image:: /code/figure/faq_site.*
 
 In the above snippet we added 2 sites: ``lat(1, 0)`` and ``lat(0, 1)``. Both
 of these sites belong to the same family, ``lat``, but have different tags:
@@ -79,11 +79,11 @@ tuples, for example lists, are not treated as hoppings.
 Starting from the example code from `What is a site?`_, we can add a hopping
 to our system in the following way:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_hopping
     :end-before: #HIDDEN_END_hopping
 
-.. image:: ../images/faq_hopping.*
+.. image:: /code/figure/faq_hopping.*
 
 Visually, a hopping is represented as a line that joins two sites.
 
@@ -139,16 +139,16 @@ Let's create two monatomic lattices (``lat_a`` and ``lat_b``).  ``(1, 0)``
 and ``(0, 1)`` will be the primitive vectors and ``(0, 0)`` and ``(0.5, 0.5)``
 the origins of the two lattices:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_lattice_monatomic
     :end-before: #HIDDEN_END_lattice_monatomic
 
-.. image:: ../images/faq_lattice.*
+.. image:: /code/figure/faq_lattice.*
 
 We can also create a ``Polyatomic`` lattice with the same primitive vectors and
 two sites in the basis:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_lattice_polyatomic
     :end-before: #HIDDEN_END_lattice_polyatomic
 
@@ -169,17 +169,17 @@ When plotting, how to color the different sublattices differently?
 ==================================================================
 In the following example we shall use a kagome lattice, which has three sublattices.
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_colors1
     :end-before: #HIDDEN_END_colors1
 
 As we can see below, we create a new plotting function that assigns a color for each family, and a different size for the hoppings depending on the family of the two sites. Finally we add sites and hoppings to our system and plot it with the new function.
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_colors2
     :end-before: #HIDDEN_END_colors2
 
-.. image:: ../images/faq_colors.*
+.. image:: /code/figure/faq_colors.*
 
 
 How to create many similar hoppings in one go?
@@ -197,31 +197,31 @@ integers).
 
 The following example shows how this can be used:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_direction1
     :end-before: #HIDDEN_END_direction1
 
-.. image:: ../images/faq_direction1.*
+.. image:: /code/figure/faq_direction1.*
 
 Note that ``HoppingKind`` only works with site families so you cannot use
 them directly with ``Polyatomic`` lattices; you have to explicitly specify
 the sublattices when creating a ``HoppingKind``:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_direction2
     :end-before: #HIDDEN_END_direction2
 
 Here, we want the hoppings between the sites from sublattice b with a direction of (0,1) in the lattice coordinates.
 
-.. image:: ../images/faq_direction2.*
+.. image:: /code/figure/faq_direction2.*
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_direction3
     :end-before: #HIDDEN_END_direction3
 
 Here, we create hoppings between the sites of the same lattice coordinates but from different families.
 
-.. image:: ../images/faq_direction3.*
+.. image:: /code/figure/faq_direction3.*
 
 
 How to set the hoppings between adjacent sites?
@@ -230,32 +230,32 @@ How to set the hoppings between adjacent sites?
 that returns a list of ``HoppingKind`` instances that connect sites with their
 (n-nearest) neighors:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_adjacent1
     :end-before: #HIDDEN_END_adjacent1
 
-.. image:: ../images/faq_adjacent1.*
-.. image:: ../images/faq_adjacent2.*
+.. image:: /code/figure/faq_adjacent1.*
+.. image:: /code/figure/faq_adjacent2.*
 
 As we can see in the figure above, ``lat.neighbors()`` (on the left) returns the hoppings between the first nearest neighbors and ``lat.neighbors(2)`` (on the right) returns the hoppings between the second nearest neighbors.
 
 When using a ``Polyatomic`` lattice ``neighbors()`` knows about the different
 sublattices:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_adjacent2
     :end-before: #HIDDEN_END_adjacent2
 
-.. image:: ../images/faq_adjacent3.*
+.. image:: /code/figure/faq_adjacent3.*
 
 However, if we use the ``neighbors()`` method of a single sublattice, we will
 only get the neighbors *on that sublattice*:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_adjacent3
     :end-before: #HIDDEN_END_adjacent3
 
-.. image:: ../images/faq_adjacent4.*
+.. image:: /code/figure/faq_adjacent4.*
 
 Note in the above that there are *only* hoppings between the blue sites. This
 is an artifact of the visualisation: the red and green sites just happen to lie
@@ -268,12 +268,12 @@ To make a hole in the system, use ``del syst[site]``, just like with any other
 mapping. In the following example we remove all sites inside some "hole"
 region:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_hole
     :end-before: #HIDDEN_END_hole
 
-.. image:: ../images/faq_hole1.*
-.. image:: ../images/faq_hole2.*
+.. image:: /code/figure/faq_hole1.*
+.. image:: /code/figure/faq_hole2.*
 
 ``del syst[site]`` also works after hoppings have been added to the system.
 If a site is deleted, then all the hoppings to/from that site are also deleted.
@@ -287,7 +287,7 @@ what is a builder?`_ if you do not know the difference).
 
 We can access the sites of a ``Builder`` by using its `~kwant.builder.Builder.sites` method:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_sites1
     :end-before: #HIDDEN_END_sites1
 
@@ -300,7 +300,7 @@ well be returned in a different order.
 After finalization, when we are dealing with a ``System``, the sites themselves
 are stored in a list, which can be accessed via the ``sites`` attribute:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_sites2
     :end-before: #HIDDEN_END_sites2
 
@@ -310,7 +310,7 @@ the system Hamiltonian, system wavefunctions etc. (see `How does Kwant order com
 ``System`` also contains the inverse mapping, ``id_by_site`` which gives us
 the index of a given site within the system:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_sites3
     :end-before: #HIDDEN_END_sites3
 
@@ -322,19 +322,19 @@ which we want to connect to leads that contain sites from a square lattice.
 
 First we construct the central system:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_different_lattice1
     :end-before: #HIDDEN_END_different_lattice1
 
-.. image:: ../images/faq_different_lattice1.*
+.. image:: /code/figure/faq_different_lattice1.*
 
 and the lead:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_different_lattice2
     :end-before: #HIDDEN_END_different_lattice2
 
-.. image:: ../images/faq_different_lattice2.*
+.. image:: /code/figure/faq_different_lattice2.*
 
 We cannot simply use `~kwant.builder.Builder.attach_lead` to attach this lead to the
 system with the honeycomb lattice because Kwant does not know how sites from
@@ -343,19 +343,19 @@ these two lattices should be connected.
 We must first add a layer of sites from the square lattice to the system and manually
 add the hoppings from these sites to the sites from the honeycomb lattice:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_different_lattice3
     :end-before: #HIDDEN_END_different_lattice3
 
-.. image:: ../images/faq_different_lattice3.*
+.. image:: /code/figure/faq_different_lattice3.*
 
 ``attach_lead()`` will now be able to attach the lead:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_different_lattice4
     :end-before: #HIDDEN_END_different_lattice4
 
-.. image:: ../images/faq_different_lattice4.*
+.. image:: /code/figure/faq_different_lattice4.*
 
 
 How to cut a finite system out of a system with translational symmetries?
@@ -369,7 +369,7 @@ will be repeated in the desired shape to create the final system.
 
 For example, say we want to create a simple model on a cubic lattice:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_fill1
     :end-before: #HIDDEN_END_fill1
 
@@ -377,21 +377,21 @@ We have now created our "template" ``Builder`` which has 3 translational
 symmetries. Next we will fill a system with no translational symmetries with
 sites and hoppings from the template inside a cuboid:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_fill2
     :end-before: #HIDDEN_END_fill2
 
-.. image:: ../images/faq_fill2.*
+.. image:: /code/figure/faq_fill2.*
 
 We can then use the original template to create a lead, which has 1 translational
 symmetry. We can then use this lead as a template to fill another section of
 the system with a cylinder of sites and hoppings:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_fill3
     :end-before: #HIDDEN_END_fill3
 
-.. image:: ../images/faq_fill3.*
+.. image:: /code/figure/faq_fill3.*
 
 
 How does Kwant order the propagating modes of a lead?
@@ -402,7 +402,7 @@ achieved with the `~kwant.system.InfiniteSystem.modes` method, which returns a
 pair of objects, the first of which contains the propagating modes of the
 system in a `~kwant.physics.PropagatingModes` object:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_pm
     :end-before: #HIDDEN_END_pm
 
@@ -411,7 +411,7 @@ modes at the requested energy (2.5 in this example).  In order to understand
 the order in which these quantities are returned it is often useful to look at
 the a section of the band structure for the system in question:
 
-.. image:: ../images/faq_pm1.*
+.. image:: /code/figure/faq_pm1.*
 
 On the above band structure we have labelled the 4 modes in the order
 that they appear in the output of ``modes()`` at energy 2.5. Note that
@@ -425,7 +425,7 @@ the modes are sorted in the following way:
 For more complicated systems and band structures this can lead to some
 possibly unintuitive orderings:
 
-.. image:: ../images/faq_pm2.*
+.. image:: /code/figure/faq_pm2.*
 
 
 How does Kwant order scattering states?
@@ -447,10 +447,10 @@ When all the site families present in a system have only 1 degree of freedom
 per site (i.e.  the all the onsites are scalars) then the index into a
 wavefunction defined over the system is exactly the site index:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_ord1
     :end-before: #HIDDEN_END_ord1
-.. literalinclude:: ../images/faq_ord1.txt
+.. literalinclude:: /code/figure/faq_ord1.txt
 
 We see that the wavefunction on a single site is a single complex number, as
 expected.
@@ -462,10 +462,10 @@ to one another.  In the case where all site families in the system have the
 wavefunction into a matrix, where the row number indexes the site, and the
 column number indexes the degree of freedom on that site:
 
-.. literalinclude:: faq.py
+.. literalinclude:: /code/include/faq.py
     :start-after: #HIDDEN_BEGIN_ord2
     :end-before: #HIDDEN_END_ord2
-.. literalinclude:: ../images/faq_ord2.txt
+.. literalinclude:: /code/figure/faq_ord2.txt
 
 We see that the wavefunction on a single site is a *vector* of 2 complex numbers,
 as we expect.
diff --git a/doc/source/tutorial/first_steps.rst b/doc/source/tutorial/first_steps.rst
index fee31fce1de21d9130b6925dcc94dbe8e0420185..9a74435f0ef9fb3784e6b638893c1e7dc2545900 100644
--- a/doc/source/tutorial/first_steps.rst
+++ b/doc/source/tutorial/first_steps.rst
@@ -62,11 +62,11 @@ Transport through a quantum wire
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/quantum_wire.py <../../../tutorial/quantum_wire.py>`
+    :download:`quantum_wire.py </code/download/quantum_wire.py>`
 
 In order to use Kwant, we need to import it:
 
-.. literalinclude:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_dwhx
     :end-before: #HIDDEN_END_dwhx
 
@@ -76,7 +76,7 @@ The first step is now the definition of the system with scattering region and
 leads. For this we make use of the `~kwant.builder.Builder` type that allows to
 define a system in a convenient way. We need to create an instance of it:
 
-.. literalinclude:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_goiq
     :end-before: #HIDDEN_END_goiq
 
@@ -92,7 +92,7 @@ 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:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_suwo
     :end-before: #HIDDEN_END_suwo
 
@@ -111,7 +111,7 @@ 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:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_zfvr
     :end-before: #HIDDEN_END_zfvr
 
@@ -140,7 +140,7 @@ 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:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_xcmc
     :end-before: #HIDDEN_END_xcmc
 
@@ -155,7 +155,7 @@ 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:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_ndez
     :end-before: #HIDDEN_END_ndez
 
@@ -163,7 +163,7 @@ 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.  The
 isolated, infinite is attached at the correct position using
 
-.. literalinclude:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_fskr
     :end-before: #HIDDEN_END_fskr
 
@@ -175,7 +175,7 @@ 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:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_xhqc
     :end-before: #HIDDEN_END_xhqc
 
@@ -187,13 +187,13 @@ you do not need to worry about that.
 Now we have finished building our system! We plot it, to make sure we didn't
 make any mistakes:
 
-.. literalinclude:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_wsgh
     :end-before: #HIDDEN_END_wsgh
 
 This should bring up this picture:
 
-.. image:: /images/quantum_wire_syst.*
+.. image:: /code/figure/quantum_wire_syst.*
 
 The system is represented in the usual way for tight-binding systems:
 dots represent the lattice points `(i, j)`, and for every
@@ -203,14 +203,14 @@ fading color.
 
 In order to use our system for a transport calculation, we need to finalize it
 
-.. literalinclude:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_dngj
     :end-before: #HIDDEN_END_dngj
 
 Having successfully created a system, we now can immediately start to compute
 its conductance as a function of energy:
 
-.. literalinclude:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_buzn
     :end-before: #HIDDEN_END_buzn
 
@@ -227,13 +227,13 @@ 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:: quantum_wire.py
+.. literalinclude:: /code/include/quantum_wire.py
     :start-after: #HIDDEN_BEGIN_lliv
     :end-before: #HIDDEN_END_lliv
 
 This should yield the result
 
-.. image:: /images/quantum_wire_result.*
+.. image:: /code/figure/quantum_wire_result.*
 
 We see a conductance quantized in units of :math:`e^2/h`,
 increasing in steps as the energy is increased. The
@@ -333,7 +333,7 @@ Building the same system with less code
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/quantum_wire_revisited.py <../../../tutorial/quantum_wire_revisited.py>`
+    :download:`quantum_wire_revisited.py </code/download/quantum_wire_revisited.py>`
 
 Kwant allows for more than one way to build a system. The reason is that
 `~kwant.builder.Builder` is essentially just a container that can be filled in
@@ -350,14 +350,14 @@ 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:: quantum_wire_revisited.py
+.. literalinclude:: /code/include/quantum_wire_revisited.py
     :start-after: #HIDDEN_BEGIN_xkzy
     :end-before: #HIDDEN_END_xkzy
 
 Previously, the scattering region was build using two ``for``-loops.
 Instead, we now write:
 
-.. literalinclude:: quantum_wire_revisited.py
+.. literalinclude:: /code/include/quantum_wire_revisited.py
     :start-after: #HIDDEN_BEGIN_vvjt
     :end-before: #HIDDEN_END_vvjt
 
@@ -375,7 +375,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:: quantum_wire_revisited.py
+.. literalinclude:: /code/include/quantum_wire_revisited.py
     :start-after: #HIDDEN_BEGIN_nooi
     :end-before: #HIDDEN_END_nooi
 
@@ -397,7 +397,7 @@ values for all the nth-nearest neighbors at once, one can similarly use
 
 The left lead is constructed in an analogous way:
 
-.. literalinclude:: quantum_wire_revisited.py
+.. literalinclude:: /code/include/quantum_wire_revisited.py
     :start-after: #HIDDEN_BEGIN_iepx
     :end-before: #HIDDEN_END_iepx
 
@@ -408,20 +408,20 @@ symmetry vector.  Here, we only construct the left lead, and use the method
 of a lead pointing in the opposite direction.  Both leads are attached as
 before and the finished system returned:
 
-.. literalinclude:: quantum_wire_revisited.py
+.. literalinclude:: /code/include/quantum_wire_revisited.py
     :start-after: #HIDDEN_BEGIN_yxot
     :end-before: #HIDDEN_END_yxot
 
 The remainder of the script has been organized into two functions.  One for the
 plotting of the conductance.
 
-.. literalinclude:: quantum_wire_revisited.py
+.. literalinclude:: /code/include/quantum_wire_revisited.py
     :start-after: #HIDDEN_BEGIN_ayuk
     :end-before: #HIDDEN_END_ayuk
 
 And one ``main`` function.
 
-.. literalinclude:: quantum_wire_revisited.py
+.. literalinclude:: /code/include/quantum_wire_revisited.py
     :start-after: #HIDDEN_BEGIN_cjel
     :end-before: #HIDDEN_END_cjel
 
@@ -429,7 +429,7 @@ Finally, we use the following standard Python construct [#]_ to execute
 ``main`` if the program is used as a script (i.e. executed as
 ``python quantum_wire_revisited.py``):
 
-.. literalinclude:: quantum_wire_revisited.py
+.. literalinclude:: /code/include/quantum_wire_revisited.py
     :start-after: #HIDDEN_BEGIN_ypbj
     :end-before: #HIDDEN_END_ypbj
 
@@ -455,7 +455,7 @@ 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:: quantum_wire_revisited.py
+     .. literalinclude:: /code/include/quantum_wire_revisited.py
          :start-after: #HIDDEN_BEGIN_vvjt
          :end-before: #HIDDEN_END_vvjt
 
diff --git a/doc/source/tutorial/graphene.rst b/doc/source/tutorial/graphene.rst
index e063940c92a6676da4626e89d0d53a0c4177d1b9..b4281a80f370128c5c9860384c33ae6a20a9cb85 100644
--- a/doc/source/tutorial/graphene.rst
+++ b/doc/source/tutorial/graphene.rst
@@ -5,7 +5,7 @@ Beyond square lattices: graphene
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/graphene.py <../../../tutorial/graphene.py>`
+    :download:`graphene.py </code/download/graphene.py>`
 
 In the following example, we are going to calculate the
 conductance through a graphene quantum dot with a p-n junction
@@ -18,7 +18,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:: graphene.py
+.. literalinclude:: /code/include/graphene.py
     :start-after: #HIDDEN_BEGIN_hnla
     :end-before: #HIDDEN_END_hnla
 
@@ -31,7 +31,7 @@ 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:: graphene.py
+.. literalinclude:: /code/include/graphene.py
     :start-after: #HIDDEN_BEGIN_shzy
     :end-before: #HIDDEN_END_shzy
 
@@ -45,7 +45,7 @@ As a next step we add the hoppings, making use of
 `~kwant.builder.HoppingKind`. For illustration purposes we define
 the hoppings ourselves instead of using ``graphene.neighbors()``:
 
-.. literalinclude:: graphene.py
+.. literalinclude:: /code/include/graphene.py
     :start-after: #HIDDEN_BEGIN_hsmc
     :end-before: #HIDDEN_END_hsmc
 
@@ -63,7 +63,7 @@ respect to the two primitive vectors ``[(1, 0), (sin_30, cos_30)]``.
 
 Adding the hoppings however still works the same way:
 
-.. literalinclude:: graphene.py
+.. literalinclude:: /code/include/graphene.py
     :start-after: #HIDDEN_BEGIN_bfwb
     :end-before: #HIDDEN_END_bfwb
 
@@ -72,7 +72,7 @@ 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:: graphene.py
+.. literalinclude:: /code/include/graphene.py
     :start-after: #HIDDEN_BEGIN_efut
     :end-before: #HIDDEN_END_efut
 
@@ -81,7 +81,7 @@ is done by the sublattices `a` and `b`.
 
 The leads are defined almost as before:
 
-.. literalinclude:: graphene.py
+.. literalinclude:: /code/include/graphene.py
     :start-after: #HIDDEN_BEGIN_aakh
     :end-before: #HIDDEN_END_aakh
 
@@ -101,14 +101,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:: graphene.py
+.. literalinclude:: /code/include/graphene.py
     :start-after: #HIDDEN_BEGIN_kmmw
     :end-before: #HIDDEN_END_kmmw
 
 The computation of some eigenvalues of the closed system is done
 in the following piece of code:
 
-.. literalinclude:: graphene.py
+.. literalinclude:: /code/include/graphene.py
     :start-after: #HIDDEN_BEGIN_zydk
     :end-before: #HIDDEN_END_zydk
 
@@ -120,7 +120,7 @@ to the previous examples, and needs not be further explained here.
 
 Finally, in the `main` function we make and plot the system:
 
-.. literalinclude:: graphene.py
+.. literalinclude:: /code/include/graphene.py
     :start-after: #HIDDEN_BEGIN_itkk
     :end-before: #HIDDEN_END_itkk
 
@@ -137,11 +137,11 @@ independent on the overall figure size.
 
 Plotting the closed system gives this result:
 
-.. image:: ../images/graphene_syst1.*
+.. image:: /code/figure/graphene_syst1.*
 
 Computing the eigenvalues of largest magnitude,
 
-.. literalinclude:: graphene.py
+.. literalinclude:: /code/include/graphene.py
     :start-after: #HIDDEN_BEGIN_jmbi
     :end-before: #HIDDEN_END_jmbi
 
@@ -151,11 +151,11 @@ should yield two eigenvalues equal to ``[ 3.07869311,
 The remaining code of `main` attaches the leads to the system and plots it
 again:
 
-.. image:: ../images/graphene_syst2.*
+.. image:: /code/figure/graphene_syst2.*
 
 It computes the band structure of one of lead 0:
 
-.. image:: ../images/graphene_bs.*
+.. image:: /code/figure/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
@@ -163,7 +163,7 @@ zero energy, due to a potential in the leads).
 
 Finally the transmission through the system is computed,
 
-.. image:: ../images/graphene_result.*
+.. image:: /code/figure/graphene_result.*
 
 showing the typical resonance-like transmission probability through
 an open quantum dot
diff --git a/doc/source/tutorial/kpm.rst b/doc/source/tutorial/kpm.rst
index 8d79ce4f20734b0a964ad66cbce8611c0b2cf137..24f0c286de3e826825eb9d8af26bfa2274410124 100644
--- a/doc/source/tutorial/kpm.rst
+++ b/doc/source/tutorial/kpm.rst
@@ -24,7 +24,7 @@ vectors.  See notes on accuracy_ below for details.
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/kernel_polynomial_method.py <../../../tutorial/kernel_polynomial_method.py>`
+    :download:`kernel_polynomial_method.py </code/download/kernel_polynomial_method.py>`
 
 .. _accuracy:
 .. specialnote:: Performance and accuracy
@@ -82,14 +82,14 @@ to obtain the density of states of a graphene disk.
 
 We start by importing kwant and defining our system.
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_sys1
     :end-before: #HIDDEN_END_sys1
 
 After making a system we can then create a `~kwant.kpm.SpectralDensity`
 object that represents the density of states for this system.
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_kpm1
     :end-before: #HIDDEN_END_kpm1
 
@@ -97,7 +97,7 @@ The `~kwant.kpm.SpectralDensity` can then be called like a function to obtain a
 sequence of energies in the spectrum of the Hamiltonian, and the corresponding
 density of states at these energies.
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_kpm2
     :end-before: #HIDDEN_END_kpm2
 
@@ -106,11 +106,11 @@ not evenly distributed over the spectrum, see Ref. [1]_ for details), however
 it is also possible to provide an explicit sequence of energies at which to
 evaluate the density of states.
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_kpm3
     :end-before: #HIDDEN_END_kpm3
 
-.. image:: ../images/kpm_dos.*
+.. image:: /code/figure/kpm_dos.*
 
 In addition to being called like functions, `~kwant.kpm.SpectralDensity`
 objects also have a method `~kwant.kpm.SpectralDensity.integrate` which can be
@@ -118,22 +118,22 @@ used to integrate the density of states against some distribution function over
 the whole spectrum. If no distribution function is specified, then the uniform
 distribution is used:
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_int1
     :end-before: #HIDDEN_END_int1
 
-.. literalinclude:: ../images/kpm_normalization.txt
+.. literalinclude:: /code/figure/kpm_normalization.txt
 
 We see that the integral of the density of states is normalized to the total
 number of available states in the system. If we wish to calculate, say, the
 number of states populated in equilibrium, then we should integrate with
 respect to a Fermi-Dirac distribution:
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_int2
     :end-before: #HIDDEN_END_int2
 
-.. literalinclude:: ../images/kpm_total_states.txt
+.. literalinclude:: /code/figure/kpm_total_states.txt
 
 .. specialnote:: Stability and performance: spectral bounds
 
@@ -161,7 +161,7 @@ exactly is changed.
 The simplest way to obtain a more accurate solution is to use the
 ``add_moments`` method:
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_acc1
     :end-before: #HIDDEN_END_acc1
 
@@ -169,16 +169,16 @@ This will update the number of calculated moments and also the default
 number of sampling points such that the maximum distance between successive
 energy points is ``energy_resolution`` (see notes on accuracy_).
 
-.. image:: ../images/kpm_dos_acc.*
+.. image:: /code/figure/kpm_dos_acc.*
 
 Alternatively, you can directly increase the number of moments
 with ``add_moments``, or the number of random vectors with ``add_vectors``.
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_acc2
     :end-before: #HIDDEN_END_acc2
 
-.. image:: ../images/kpm_dos_r.*
+.. image:: /code/figure/kpm_dos_r.*
 
 
 .. _operator_spectral_density:
@@ -200,14 +200,14 @@ the spectral density of the given operator that is calculated.
 If an explicit matrix is provided then it must have the same
 shape as the system Hamiltonian.
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_op1
     :end-before: #HIDDEN_END_op1
 
 
 Or, to do the same calculation using `kwant.operator.Density`:
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_op2
     :end-before: #HIDDEN_END_op2
 
@@ -215,7 +215,7 @@ Using operators from `kwant.operator` allows us to calculate quantities
 such as the *local* density of states by telling the operator not to
 sum over all the sites of the system:
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_op3
     :end-before: #HIDDEN_END_op3
 
@@ -223,11 +223,11 @@ sum over all the sites of the system:
 which allows us to plot the local density of states at different
 point in the spectrum:
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_op4
     :end-before: #HIDDEN_END_op4
 
-.. image:: ../images/kpm_ldos.*
+.. image:: /code/figure/kpm_ldos.*
 
 This nicely illustrates the edge states of the graphene dot at zero
 energy, and the bulk states at higher energy.
@@ -248,7 +248,7 @@ To change how the random vectors are generated, you need only specify a
 function that takes the dimension of the Hilbert space as a single parameter,
 and which returns a vector in that Hilbert space:
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_fact1
     :end-before: #HIDDEN_END_fact1
 
@@ -275,7 +275,7 @@ first argument, and linear in its second argument. Below, we compare two
 methods for computing the local density of states, one using
 `kwant.operator.Density`, and the other using a custom function.
 
-.. literalinclude:: kernel_polynomial_method.py
+.. literalinclude:: /code/include/kernel_polynomial_method.py
     :start-after: #HIDDEN_BEGIN_blm
     :end-before: #HIDDEN_END_blm
 
diff --git a/doc/source/tutorial/operators.rst b/doc/source/tutorial/operators.rst
index c2c314acd3f9a89464ef2ea67d503fbe4d111b4e..c7e66feaf40efcd00339571bef2efa462d217838 100644
--- a/doc/source/tutorial/operators.rst
+++ b/doc/source/tutorial/operators.rst
@@ -13,7 +13,7 @@ texture.
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/magnetic_texture.py <../../../tutorial/magnetic_texture.py>`
+    :download:`magnetic_texture.py </code/download/magnetic_texture.py>`
 
 
 Introduction
@@ -47,26 +47,26 @@ of site :math:`i`, and :math:`r_i = \sqrt{x_i^2 + y_i^2}`.
 To define this model in Kwant we start as usual by defining functions that
 depend on the model parameters:
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_model
     :end-before: #HIDDEN_END_model
 
 and define our system as a square shape on a square lattice with two orbitals
 per site, with leads attached on the left and right:
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_syst
     :end-before: #HIDDEN_END_syst
 
 Below is a plot of a projection of :math:`\mathbf{m}_i` onto the x-y plane
 inside the scattering region. The z component is shown by the color scale:
 
-.. image:: ../images/mag_field_direction.*
+.. image:: /code/figure/mag_field_direction.*
 
 We will now be interested in analyzing the form of the scattering states
 that originate from the left lead:
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_wavefunction
     :end-before: #HIDDEN_END_wavefunction
 
@@ -82,7 +82,7 @@ When there are multiple degrees of freedom per site, however, one has to be
 more careful. In the present case with two (spin) degrees of freedom per site
 one could calculate the per-site density like:
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_ldos
     :end-before: #HIDDEN_END_ldos
 
@@ -91,14 +91,14 @@ local quantities we can meaningfully compute. For example, we may wish to
 calculate the local z-projected spin density. We could calculate
 this in the following way:
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_lsdz
     :end-before: #HIDDEN_END_lsdz
 
 If we wanted instead to calculate the local y-projected spin density, we would
 need to use an even more complicated expression:
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_lsdy
     :end-before: #HIDDEN_END_lsdy
 
@@ -107,7 +107,7 @@ book-keeping by providing a simple interface for defining operators that act on
 wavefunctions. To calculate the above quantities we would use the
 `~kwant.operator.Density` operator like so:
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_lden
     :end-before: #HIDDEN_END_lden
 
@@ -126,7 +126,7 @@ Below we can see colorplots of the above-calculated quantities. The array that
 is returned by evaluating a `~kwant.operator.Density` can be used directly with
 `kwant.plotter.map`:
 
-.. image:: ../images/spin_densities.*
+.. image:: /code/figure/spin_densities.*
 
 
 .. specialnote:: Technical Details
@@ -180,14 +180,14 @@ where :math:`\mathbf{H}_{ab}` is the hopping matrix from site :math:`b` to site
 :math:`a`.  For example, to calculate the local current and
 spin current:
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_current
     :end-before: #HIDDEN_END_current
 
 Evaluating a `~kwant.operator.Current` operator on a wavefunction returns a
 1D array of values that can be directly used with `kwant.plotter.current`:
 
-.. image:: ../images/spin_currents.*
+.. image:: /code/figure/spin_currents.*
 
 .. note::
 
@@ -262,7 +262,7 @@ scattering region.
 Doing this is as simple as passing a *function* when instantiating
 the `~kwant.operator.Current`, instead of a constant matrix:
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_following
     :end-before: #HIDDEN_END_following
 
@@ -286,7 +286,7 @@ Using this we can see that the spin current is essentially oriented along
 the direction of :math:`m_i` in the present regime where the onsite term
 in the Hamiltonian is dominant:
 
-.. image:: ../images/spin_current_comparison.*
+.. image:: /code/figure/spin_current_comparison.*
 
 .. note:: Although this example used exclusively `~kwant.operator.Current`,
           you can do the same with `~kwant.operator.Density`.
@@ -308,11 +308,11 @@ Density of states in a circle
 To calculate the density of states inside a circle of radius
 20 we can simply do:
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_density_cut
     :end-before: #HIDDEN_END_density_cut
 
-.. literalinclude:: ../images/circle_dos.txt
+.. literalinclude:: /code/figure/circle_dos.txt
 
 note that we also provide ``sum=True``, which means that evaluating the
 operator on a wavefunction will produce a single scalar. This is semantically
@@ -325,11 +325,11 @@ Current flowing through a cut
 Below we calculate the probability current and z-projected spin current near
 the interfaces with the left and right leads.
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_current_cut
     :end-before: #HIDDEN_END_current_cut
 
-.. literalinclude:: ../images/current_cut.txt
+.. literalinclude:: /code/figure/current_cut.txt
 
 We see that the probability current is conserved across the scattering region,
 but the z-projected spin current is not due to the fact that the Hamiltonian
@@ -360,8 +360,8 @@ This can be achieved with `~kwant.operator.Current.bind`:
              *different* set of parameters. This will almost certainly give
              incorrect results.
 
-.. literalinclude:: magnetic_texture.py
+.. literalinclude:: /code/include/magnetic_texture.py
     :start-after: #HIDDEN_BEGIN_bind
     :end-before: #HIDDEN_END_bind
 
-.. image:: ../images/bound_current.*
+.. image:: /code/figure/bound_current.*
diff --git a/doc/source/tutorial/plotting.rst b/doc/source/tutorial/plotting.rst
index 950ca0797645722428efab62caf0417cec6c1009..f8bb6586a1c2b7d6cd30cf084d983b5eef8b3ce5 100644
--- a/doc/source/tutorial/plotting.rst
+++ b/doc/source/tutorial/plotting.rst
@@ -12,13 +12,13 @@ these options can be used to achieve various very different objectives.
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/plot_graphene.py <../../../tutorial/plot_graphene.py>`
+    :download:`plot_graphene.py </code/download/plot_graphene.py>`
 
 We begin by first considering a circular graphene quantum dot (similar to what
 has been used in parts of the tutorial :ref:`tutorial-graphene`.)  In contrast
 to previous examples, we will also use hoppings beyond next-nearest neighbors:
 
-.. literalinclude:: plot_graphene.py
+.. literalinclude:: /code/include/plot_graphene.py
     :start-after: #HIDDEN_BEGIN_makesyst
     :end-before: #HIDDEN_END_makesyst
 
@@ -31,14 +31,14 @@ next-nearest-neighbor hopping [#]_.
 
 Of course, the system can be plotted simply with default settings:
 
-.. literalinclude:: plot_graphene.py
+.. literalinclude:: /code/include/plot_graphene.py
     :start-after: #HIDDEN_BEGIN_plotsyst1
     :end-before: #HIDDEN_END_plotsyst1
 
 However, due to the richer structure of the lattice, this results in a rather
 busy plot:
 
-.. image:: ../images/plot_graphene_syst1.*
+.. image:: /code/figure/plot_graphene_syst1.*
 
 A much clearer plot can be obtained by using different colors for both
 sublattices, and by having different line widths for different hoppings.  This
@@ -47,7 +47,7 @@ can be achieved by passing a function to the arguments of
 must be a function taking one site as argument, for hoppings a function taking
 the start end end site of hopping as arguments:
 
-.. literalinclude:: plot_graphene.py
+.. literalinclude:: /code/include/plot_graphene.py
     :start-after: #HIDDEN_BEGIN_plotsyst2
     :end-before: #HIDDEN_END_plotsyst2
 
@@ -55,7 +55,7 @@ Note that since we are using an unfinalized Builder, a ``site`` is really an
 instance of `~kwant.builder.Site`. With these adjustments we arrive at a plot
 that carries the same information, but is much easier to interpret:
 
-.. image:: ../images/plot_graphene_syst2.*
+.. image:: /code/figure/plot_graphene_syst2.*
 
 Apart from plotting the *system* itself, `~kwant.plotter.plot` can also be used
 to plot *data* living on the system.
@@ -67,7 +67,7 @@ nearest-neighbor hopping.  Computing the wave functions is done in the usual
 way (note that for a large-scale system, one would probably want to use sparse
 linear algebra):
 
-.. literalinclude:: plot_graphene.py
+.. literalinclude:: /code/include/plot_graphene.py
     :start-after: #HIDDEN_BEGIN_plotdata1
     :end-before: #HIDDEN_END_plotdata1
 
@@ -75,7 +75,7 @@ In most cases, to plot the wave function probability, one wouldn't use
 `~kwant.plotter.plot`, but rather `~kwant.plotter.map`. Here, we plot the
 `n`-th wave function using it:
 
-.. literalinclude:: plot_graphene.py
+.. literalinclude:: /code/include/plot_graphene.py
     :start-after: #HIDDEN_BEGIN_plotdata2
     :end-before: #HIDDEN_END_plotdata2
 
@@ -83,7 +83,7 @@ This results in a standard pseudocolor plot, showing in this case (``n=225``) a
 graphene edge state, i.e. a wave function mostly localized at the zigzag edges
 of the quantum dot.
 
-.. image:: ../images/plot_graphene_data1.*
+.. image:: /code/figure/plot_graphene_data1.*
 
 However although in general preferable, `~kwant.plotter.map` has a few
 deficiencies for this small system: For example, there are a few distortions at
@@ -99,7 +99,7 @@ the symbol shape depending on the sublattice. With a triangle pointing up and
 down on the respective sublattice, the symbols used by plot fill the space
 completely:
 
-.. literalinclude:: plot_graphene.py
+.. literalinclude:: /code/include/plot_graphene.py
     :start-after: #HIDDEN_BEGIN_plotdata3
     :end-before: #HIDDEN_END_plotdata3
 
@@ -116,7 +116,7 @@ is represented by an integer. In order to obtain the original
 
 With this we arrive at
 
-.. image:: ../images/plot_graphene_data2.*
+.. image:: /code/figure/plot_graphene_data2.*
 
 with the same information as `~kwant.plotter.map`, but with a cleaner look.
 
@@ -125,7 +125,7 @@ are best visible in a given plot. With `~kwant.plotter.plot` one can easily go
 beyond pseudocolor-like plots. For example, we can represent the wave function
 probability using the symbols itself:
 
-.. literalinclude:: plot_graphene.py
+.. literalinclude:: /code/include/plot_graphene.py
     :start-after: #HIDDEN_BEGIN_plotdata4
     :end-before: #HIDDEN_END_plotdata4
 
@@ -135,7 +135,7 @@ visible. The hoppings are also plotted in order to show the underlying lattice.
 
 With this, we arrive at
 
-.. image:: ../images/plot_graphene_data3.*
+.. image:: /code/figure/plot_graphene_data3.*
 
 which shows the edge state nature of the wave function most clearly.
 
@@ -150,7 +150,7 @@ which shows the edge state nature of the wave function most clearly.
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/plot_zincblende.py <../../../tutorial/plot_zincblende.py>`
+    :download:`plot_zincblende.py </code/download/plot_zincblende.py>`
 
 Zincblende is a very common crystal structure of semiconductors. It is a
 face-centered cubic crystal with two inequivalent atoms in the unit cell
@@ -159,7 +159,7 @@ structure, but two equivalent atoms per unit cell).
 
 It is very easily generated in Kwant with `kwant.lattice.general`:
 
-.. literalinclude:: plot_zincblende.py
+.. literalinclude:: /code/include/plot_zincblende.py
     :start-after: #HIDDEN_BEGIN_zincblende1
     :end-before: #HIDDEN_END_zincblende1
 
@@ -168,7 +168,7 @@ Note how we keep references to the two different sublattices for later use.
 A three-dimensional structure is created as easily as in two dimensions, by
 using the `~kwant.lattice.PolyatomicLattice.shape`-functionality:
 
-.. literalinclude:: plot_zincblende.py
+.. literalinclude:: /code/include/plot_zincblende.py
     :start-after: #HIDDEN_BEGIN_zincblende2
     :end-before: #HIDDEN_END_zincblende2
 
@@ -179,13 +179,13 @@ real calculation, several atomic orbitals would have to be considered).
 `~kwant.plotter.plot` can plot 3D systems just as easily as its two-dimensional
 counterparts:
 
-.. literalinclude:: plot_zincblende.py
+.. literalinclude:: /code/include/plot_zincblende.py
     :start-after: #HIDDEN_BEGIN_plot1
     :end-before: #HIDDEN_END_plot1
 
 resulting in
 
-.. image:: ../images/plot_zincblende_syst1.*
+.. image:: /code/figure/plot_zincblende_syst1.*
 
 You might notice that the standard options for plotting are quite different in
 3D than in 2D. For example, by default hoppings are not printed, but sites are
@@ -207,14 +207,14 @@ Also for 3D it is possible to customize the plot. For example, we
 can explicitly plot the hoppings as lines, and color sites differently
 depending on the sublattice:
 
-.. literalinclude:: plot_zincblende.py
+.. literalinclude:: /code/include/plot_zincblende.py
     :start-after: #HIDDEN_BEGIN_plot2
     :end-before: #HIDDEN_END_plot2
 
 which results in a 3D plot that allows to interactively (when plotted
 in a window) explore the crystal structure:
 
-.. image:: ../images/plot_zincblende_syst2.*
+.. image:: /code/figure/plot_zincblende_syst2.*
 
 Hence, a few lines of code using Kwant allow to explore all the different
 crystal lattices out there!
diff --git a/doc/source/tutorial/spectrum.rst b/doc/source/tutorial/spectrum.rst
index b119e3c79dc000d60f11e717c812e2cf8ef093ac..46843440019be37e7bc0d3e99124366eb6df68ee 100644
--- a/doc/source/tutorial/spectrum.rst
+++ b/doc/source/tutorial/spectrum.rst
@@ -6,7 +6,7 @@ Band structure calculations
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/band_structure.py <../../../tutorial/band_structure.py>`
+    :download:`band_structure.py </code/download/band_structure.py>`
 
 When doing transport simulations, one also often needs to know the band
 structure of the leads, i.e. the energies of the propagating plane waves in the
@@ -19,7 +19,7 @@ tight-binding wire.
 Computing band structures in Kwant is easy. Just define a lead in the
 usual way:
 
-.. literalinclude:: band_structure.py
+.. literalinclude:: /code/include/band_structure.py
     :start-after: #HIDDEN_BEGIN_zxip
     :end-before: #HIDDEN_END_zxip
 
@@ -42,13 +42,13 @@ do something more profound with the dispersion relation these energies may be
 calculated directly using `~kwant.physics.Bands`. For now we just plot the
 bandstructure:
 
-.. literalinclude:: band_structure.py
+.. literalinclude:: /code/include/band_structure.py
     :start-after: #HIDDEN_BEGIN_pejz
     :end-before: #HIDDEN_END_pejz
 
 This gives the result:
 
-.. image:: ../images/band_structure_result.*
+.. image:: /code/figure/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
@@ -61,7 +61,7 @@ Closed systems
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/closed_system.py <../../../tutorial/closed_system.py>`
+    :download:`closed_system.py </code/download/closed_system.py>`
 
 Although Kwant is (currently) mainly aimed towards transport problems, it
 can also easily be used to compute properties of closed systems -- after
@@ -74,14 +74,14 @@ To compute the eigenenergies and eigenstates, we will make use of the sparse
 linear algebra functionality of `SciPy <https://www.scipy.org>`_, which
 interfaces the ARPACK package:
 
-.. literalinclude:: closed_system.py
+.. literalinclude:: /code/include/closed_system.py
     :start-after: #HIDDEN_BEGIN_tibv
     :end-before: #HIDDEN_END_tibv
 
 We set up the system using the `shape`-function as in
 :ref:`tutorial-abring`, but do not add any leads:
 
-.. literalinclude:: closed_system.py
+.. literalinclude:: /code/include/closed_system.py
     :start-after: #HIDDEN_BEGIN_qlyd
     :end-before: #HIDDEN_END_qlyd
 
@@ -93,14 +93,14 @@ 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:: closed_system.py
+.. literalinclude:: /code/include/closed_system.py
     :start-after: #HIDDEN_BEGIN_yvri
     :end-before: #HIDDEN_END_yvri
 
 Note that we use sparse linear algebra to efficiently calculate only a
 few lowest eigenvalues. Finally, we obtain the result:
 
-.. image:: ../images/closed_system_result.*
+.. image:: /code/figure/closed_system_result.*
 
 At zero magnetic field several energy levels are degenerate (since our
 quantum dot is rather symmetric). These degeneracies are split
@@ -110,11 +110,11 @@ Landau level energies at higher magnetic fields [#]_.
 The eigenvectors are obtained very similarly, and can be plotted directly
 using `~kwant.plotter.map`:
 
-.. literalinclude:: closed_system.py
+.. literalinclude:: /code/include/closed_system.py
     :start-after: #HIDDEN_BEGIN_wave
     :end-before: #HIDDEN_END_wave
 
-.. image:: ../images/closed_system_eigenvector.*
+.. image:: /code/figure/closed_system_eigenvector.*
 
 The last two arguments to `~kwant.plotter.map` are optional.  The first prevents
 a colorbar from appearing.  The second, ``oversampling=1``, makes the image look
@@ -127,11 +127,11 @@ that they can have *non-zero* local current. We can calculate the local
 current due to a state by using `kwant.operator.Current` and plotting
 it using `kwant.plotter.current`:
 
-.. literalinclude:: closed_system.py
+.. literalinclude:: /code/include/closed_system.py
     :start-after: #HIDDEN_BEGIN_current
     :end-before: #HIDDEN_END_current
 
-.. image:: ../images/closed_system_current.*
+.. image:: /code/figure/closed_system_current.*
 
 .. specialnote:: Technical details
 
diff --git a/doc/source/tutorial/spin_potential_shape.rst b/doc/source/tutorial/spin_potential_shape.rst
index abe776ef24944dd99bdb43c0e22287a1c7f611f8..bcf95dacc7627cf3822474009480287802cc6889 100644
--- a/doc/source/tutorial/spin_potential_shape.rst
+++ b/doc/source/tutorial/spin_potential_shape.rst
@@ -11,7 +11,7 @@ Matrix structure of on-site and hopping elements
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/spin_orbit.py <../../../tutorial/spin_orbit.py>`
+    :download:`spin_orbit.py </code/download/spin_orbit.py>`
 
 We begin by extending the simple 2DEG-Hamiltonian by a Rashba spin-orbit
 coupling and a Zeeman splitting due to an external magnetic field:
@@ -42,14 +42,14 @@ use matrices in our program, we import the Tinyarray package.  (`NumPy
 <http://www.numpy.org/>`_ would work as well, but Tinyarray is much faster
 for small arrays.)
 
-.. literalinclude:: spin_orbit.py
+.. literalinclude:: /code/include/spin_orbit.py
     :start-after: #HIDDEN_BEGIN_xumz
     :end-before: #HIDDEN_END_xumz
 
 For convenience, we define the Pauli-matrices first (with :math:`\sigma_0` the
 unit matrix):
 
-.. literalinclude:: spin_orbit.py
+.. literalinclude:: /code/include/spin_orbit.py
     :start-after: #HIDDEN_BEGIN_hwbt
     :end-before: #HIDDEN_END_hwbt
 
@@ -57,7 +57,7 @@ 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:: spin_orbit.py
+.. literalinclude:: /code/include/spin_orbit.py
     :start-after: #HIDDEN_BEGIN_uxrm
     :end-before: #HIDDEN_END_uxrm
 
@@ -85,14 +85,14 @@ when specifying `(1, 0)` it is not necessary to specify `(-1, 0)`),
 
 The leads also allow for a matrix structure,
 
-.. literalinclude:: spin_orbit.py
+.. literalinclude:: /code/include/spin_orbit.py
     :start-after: #HIDDEN_BEGIN_yliu
     :end-before: #HIDDEN_END_yliu
 
 The remainder of the code is unchanged, and as a result we should obtain
 the following, clearly non-monotonic conductance steps:
 
-.. image:: ../images/spin_orbit_result.*
+.. image:: /code/figure/spin_orbit_result.*
 
 .. specialnote:: Technical details
 
@@ -129,7 +129,7 @@ Spatially dependent values through functions
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/quantum_well.py <../../../tutorial/quantum_well.py>`
+    :download:`quantum_well.py </code/download/quantum_well.py>`
 
 Up to now, all examples had position-independent matrix-elements
 (and thus translational invariance along the wire, which
@@ -148,7 +148,7 @@ 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:: quantum_well.py
+.. literalinclude:: /code/include/quantum_well.py
     :start-after: #HIDDEN_BEGIN_ehso
     :end-before: #HIDDEN_END_ehso
 
@@ -161,7 +161,7 @@ namespace of `make_system`.
 Kwant now allows us to pass a function as a value to
 `~kwant.builder.Builder`:
 
-.. literalinclude:: quantum_well.py
+.. literalinclude:: /code/include/quantum_well.py
     :start-after: #HIDDEN_BEGIN_coid
     :end-before: #HIDDEN_END_coid
 
@@ -180,7 +180,7 @@ of the lead -- this should be kept in mind.
 
 Finally, we compute the transmission probability:
 
-.. literalinclude:: quantum_well.py
+.. literalinclude:: /code/include/quantum_well.py
     :start-after: #HIDDEN_BEGIN_sqvr
     :end-before: #HIDDEN_END_sqvr
 
@@ -189,7 +189,7 @@ additional arguments to the functions that provide the Hamiltonian matrix
 elements.  In this example we are able to solve the system for different depths
 of the potential well by passing the potential value. We obtain the result:
 
-.. image:: ../images/quantum_well_result.*
+.. image:: /code/figure/quantum_well_result.*
 
 Starting from no potential (well depth = 0), we observe the typical
 oscillatory transmission behavior through resonances in the quantum well.
@@ -216,13 +216,13 @@ Nontrivial shapes
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/ab_ring.py <../../../tutorial/ab_ring.py>`
+    :download:`ab_ring.py </code/download/ab_ring.py>`
 
 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/ab_ring_sketch.*
+.. image:: /code/figure/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
@@ -236,7 +236,7 @@ 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:: ab_ring.py
+.. literalinclude:: /code/include/ab_ring.py
     :start-after: #HIDDEN_BEGIN_eusz
     :end-before: #HIDDEN_END_eusz
 
@@ -247,7 +247,7 @@ 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:: ab_ring.py
+.. literalinclude:: /code/include/ab_ring.py
     :start-after: #HIDDEN_BEGIN_lcak
     :end-before: #HIDDEN_END_lcak
 
@@ -262,7 +262,7 @@ along the branch cut in the lower arm of the ring. For this we select
 all hoppings in x-direction that are of the form `(lat(1, j), lat(0, j))`
 with ``j<0``:
 
-.. literalinclude:: ab_ring.py
+.. literalinclude:: /code/include/ab_ring.py
     :start-after: #HIDDEN_BEGIN_lvkt
     :end-before: #HIDDEN_END_lvkt
 
@@ -277,7 +277,7 @@ by the parameter `phi`.
 
 For the leads, we can also use the ``lat.shape``-functionality:
 
-.. literalinclude:: ab_ring.py
+.. literalinclude:: /code/include/ab_ring.py
     :start-after: #HIDDEN_BEGIN_qwgr
     :end-before: #HIDDEN_END_qwgr
 
@@ -288,7 +288,7 @@ no restriction on ``x`` in ``lead_shape``) [#]_.
 
 Attaching the leads is done as before:
 
-.. literalinclude:: ab_ring.py
+.. literalinclude:: /code/include/ab_ring.py
     :start-after: #HIDDEN_BEGIN_skbz
     :end-before: #HIDDEN_END_skbz
 
@@ -303,16 +303,16 @@ 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/ab_ring_sketch2.*
+.. image:: /code/figure/ab_ring_sketch2.*
 
 After the lead has been attached, the system should look like this:
 
-.. image:: ../images/ab_ring_syst.*
+.. image:: /code/figure/ab_ring_syst.*
 
 The computation of the conductance goes in the same fashion as before.
 Finally you should get the following result:
 
-.. image:: ../images/ab_ring_result.*
+.. image:: /code/figure/ab_ring_result.*
 
 where one can observe the conductance oscillations with the
 period of one flux quantum.
@@ -328,7 +328,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/ab_ring_note1.*
+    .. image:: /code/figure/ab_ring_note1.*
 
   - Per default, `~kwant.builder.Builder.attach_lead` attaches
     the lead to the "outside" of the structure, by tracing the
@@ -343,7 +343,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/ab_ring_note2.*
+    .. image:: /code/figure/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/superconductors.rst b/doc/source/tutorial/superconductors.rst
index 67738b2f91e30e112e073267253e62dc3650b0ff..6113cb78a7be06bf57b7bf7f39b65095505b1f29 100644
--- a/doc/source/tutorial/superconductors.rst
+++ b/doc/source/tutorial/superconductors.rst
@@ -3,7 +3,7 @@ Superconductors: orbital degrees of freedom, conservation laws and symmetries
 
 .. seealso::
     The complete source code of this example can be found in
-    :download:`tutorial/superconductor.py <../../../tutorial/superconductor.py>`
+    :download:`superconductor.py </code/download/superconductor.py>`
 
 This example deals with superconductivity on the level of the
 Bogoliubov-de Gennes (BdG) equation. In this framework, the Hamiltonian
@@ -51,14 +51,14 @@ of freedom.
 Let us consider a system that consists of a normal lead on the left,
 a superconductor on the right, and a tunnel barrier in between:
 
-.. image:: ../images/superconductor_transport_sketch.*
+.. image:: /code/figure/superconductor_transport_sketch.*
 
 We implement the BdG Hamiltonian in Kwant using a 2x2 matrix structure
 for all Hamiltonian matrix elements, as we did
 previously in the :ref:`spin example <tutorial_spinorbit>`. We declare
 the square lattice and construct the scattering region with the following:
 
-.. literalinclude:: superconductor.py
+.. literalinclude:: /code/include/superconductor.py
     :start-after: #HIDDEN_BEGIN_nbvn
     :end-before: #HIDDEN_END_nbvn
 
@@ -80,7 +80,7 @@ of it). The next step towards computing conductance is to attach leads.
 Let's attach two leads: a normal one to the left end, and a superconducting
 one to the right end. Starting with the left lead, we have:
 
-.. literalinclude:: superconductor.py
+.. literalinclude:: /code/include/superconductor.py
     :start-after: #HIDDEN_BEGIN_ttth
     :end-before: #HIDDEN_END_ttth
 
@@ -120,7 +120,7 @@ of ascending eigenvalues of the conservation law.
 In order to move on with the conductance calculation, let's attach the second
 lead to the right side of the scattering region:
 
-.. literalinclude:: superconductor.py
+.. literalinclude:: /code/include/superconductor.py
     :start-after: #HIDDEN_BEGIN_zuuw
     :end-before: #HIDDEN_END_zuuw
 
@@ -140,7 +140,7 @@ confused by the fact that it says ``transmission`` -- transmission
 into the same lead is reflection), and reflection from electrons to holes
 is ``smatrix.transmission((0, 1), (0, 0))``:
 
-.. literalinclude:: superconductor.py
+.. literalinclude:: /code/include/superconductor.py
     :start-after: #HIDDEN_BEGIN_jbjt
     :end-before: #HIDDEN_END_jbjt
 
@@ -150,7 +150,7 @@ reflection of electrons to electrons, and from its size we can extract the numbe
 
 For the default parameters, we obtain the following conductance:
 
-.. image:: ../images/superconductor_transport_result.*
+.. image:: /code/figure/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
@@ -176,13 +176,13 @@ the electron and hole blocks. The exception is of course at zero energy, in whic
 case particle-hole symmetry transforms between the electron and hole blocks, resulting
 in a symmetric scattering matrix. We can check the symmetry explicitly with
 
-.. literalinclude:: superconductor.py
+.. literalinclude:: /code/include/superconductor.py
     :start-after: #HIDDEN_BEGIN_pqmp
     :end-before: #HIDDEN_END_pqmp
 
 which yields the output
 
-.. literalinclude:: ../images/check_PHS_out.txt
+.. literalinclude:: /code/figure/check_PHS_out.txt
 
 Note that :math:`\mathcal{P}` flips the sign of momentum, and for the parameters
 we consider here, there are two electron and two hole modes active at zero energy.
diff --git a/setup.py b/setup.py
index 70a9f19ccddf0e001d0d76af6d7e22c21632d34a..197505db81f8a782506519d5a4f70687f713aa84 100755
--- a/setup.py
+++ b/setup.py
@@ -227,40 +227,7 @@ Build configuration was:
 """
 
 
-class build_tut(Command):
-    description = "build the tutorial scripts"
-    user_options = []
-
-    def initialize_options(self):
-        pass
-
-    def finalize_options(self):
-        pass
-
-    def run(self):
-        tut_dir = 'tutorial'
-        if not os.path.exists(tut_dir):
-            os.mkdir(tut_dir)
-        for in_fname in glob.glob('doc/source/tutorial/*.py'):
-            basename = os.path.basename(in_fname)
-            if basename == 'faq.py':
-                # The FAQ script is not meant as an example.
-                continue
-            out_fname = os.path.join(tut_dir, basename)
-            with open(in_fname, 'rb') as in_file:
-                with open(out_fname, 'wb') as out_file:
-                    for line in in_file:
-                        if not line.startswith(b'#HIDDEN'):
-                            out_file.write(line)
-
-
-# Our version of the "build" command also makes sure the tutorial is made.
-# Even though the tutorial is not necessary for installation, and "build" is
-# supposed to make everything needed to install, this is a robust way to ensure
-# that the tutorial is present.
 class build(build_orig):
-    sub_commands = [('build_tut', None)] + build_orig.sub_commands
-
     def run(self):
         super().run()
         write_version(os.path.join(self.build_lib, *STATIC_VERSION_PATH))
@@ -638,7 +605,6 @@ def main():
           cmdclass={'build': build,
                     'sdist': sdist,
                     'build_ext': build_ext,
-                    'build_tut': build_tut,
                     'test': test},
           ext_modules=exts,
           install_requires=['numpy >= 1.8.1', 'scipy >= 0.14',