diff --git a/doc/source/images/tutorial5a.py b/doc/source/images/tutorial5a.py
index 8d54f0efcdb2681e6a700a66e7985d709c0fb96c..359c3a3530cdc91c50207d76662a6a9fd9dcbe38 100644
--- a/doc/source/images/tutorial5a.py
+++ b/doc/source/images/tutorial5a.py
@@ -40,8 +40,7 @@ def make_lead(a=1, t=1.0, mu=0.7, Delta=0.1, W=10):
 
         lead[(1, j), (0, j)] = - t * tau_z
 
-    # return a finalized lead
-    return lead.finalized()
+    return lead
 
 
 def plot_bandstructure(flead, momenta):
@@ -66,7 +65,8 @@ def plot_bandstructure(flead, momenta):
 
 
 def main():
-    flead = make_lead()
+    # Make system and finalize it right away.
+    flead = make_lead().finalized()
 
     # list of momenta at which the bands should be computed
     momenta = np.arange(-1.5, 1.5 + .0001, 0.002 * pi)
diff --git a/doc/source/images/tutorial5b.py b/doc/source/images/tutorial5b.py
index 322b120130181bdab19cd82b24d88c8eecfa97f6..1ed8213f8797b1d911938bb064a9948f585b436d 100644
--- a/doc/source/images/tutorial5b.py
+++ b/doc/source/images/tutorial5b.py
@@ -80,12 +80,12 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
     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
 
-    #### Attach the leads and return the finalized system. ####
+    #### Attach the leads and return the system. ####
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
     sys.attach_lead(lead2)
 
-    return sys.finalized()
+    return sys
 
 def plot_conductance(fsys, energies):
     # Compute conductance
@@ -114,7 +114,7 @@ def plot_conductance(fsys, energies):
 
 
 def main():
-    fsys = make_system()
+    fsys = make_system().finalized()
 
     plot_conductance(fsys, energies=[0.002 * i for i in xrange(100)])
 
diff --git a/doc/source/tutorial/tutorial5.rst b/doc/source/tutorial/tutorial5.rst
index bd572b00476842405fc052d35730dada2782707f..c7f5a9a35465dabc3ba00cf84b3188d811786a36 100644
--- a/doc/source/tutorial/tutorial5.rst
+++ b/doc/source/tutorial/tutorial5.rst
@@ -31,7 +31,7 @@ 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
-    :lines: 21-45
+    :lines: 21-42
 
 As you see, the example is syntactically equivalent to our
 :ref:`spin example <tutorial_spinorbit>`, the only difference
diff --git a/examples/tutorial5a.py b/examples/tutorial5a.py
index 995c5b960364c9e7511620bdc4266f3a579810f2..60bf12ef9b413df2d6f312c8dba9b3e6fda300ea 100644
--- a/examples/tutorial5a.py
+++ b/examples/tutorial5a.py
@@ -39,8 +39,7 @@ def make_lead(a=1, t=1.0, mu=0.7, Delta=0.1, W=10):
 
         lead[(1, j), (0, j)] = - t * tau_z
 
-    # return a finalized lead
-    return lead.finalized()
+    return lead
 
 
 def plot_bandstructure(flead, momenta):
@@ -56,7 +55,8 @@ def plot_bandstructure(flead, momenta):
 
 
 def main():
-    flead = make_lead()
+    # Make system and finalize it right away.
+    flead = make_lead().finalized()
 
     # list of momenta at which the bands should be computed
     momenta = np.arange(-1.5, 1.5 + .0001, 0.002 * pi)
diff --git a/examples/tutorial5b.py b/examples/tutorial5b.py
index 19bd84662ecd9f4f678a01e7f4a754361945d787..2771599c90a68007e4e5d1008f19583cb8de755b 100644
--- a/examples/tutorial5b.py
+++ b/examples/tutorial5b.py
@@ -79,12 +79,12 @@ def make_system(a=1, W=10, L=10, barrier=1.5, barrierpos=(3, 4),
     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
 
-    #### Attach the leads and return the finalized system. ####
+    #### Attach the leads and return the system. ####
     sys.attach_lead(lead0)
     sys.attach_lead(lead1)
     sys.attach_lead(lead2)
 
-    return sys.finalized()
+    return sys
 
 def plot_conductance(fsys, energies):
     # Compute conductance
@@ -103,10 +103,13 @@ def plot_conductance(fsys, energies):
 
 
 def main():
-    fsys = make_system()
+    sys = make_system()
 
     # Check that the system looks as intended.
-    kwant.plot(fsys)
+    kwant.plot(sys)
+
+    # Finalize the system.
+    fsys = sys.finalized()
 
     plot_conductance(fsys, energies=[0.002 * i for i in xrange(100)])