diff --git a/doc/source/images/closed_system.py.diff b/doc/source/images/closed_system.py.diff
index 7c858ecf565831815ef5a0ee7ca7a8492f20d91e..4c2816d01d5f9c80c3cb3979a4da3c8bbe6e6e5b 100644
--- a/doc/source/images/closed_system.py.diff
+++ b/doc/source/images/closed_system.py.diff
@@ -1,6 +1,6 @@
 --- original
 +++ modified
-@@ -17,6 +17,7 @@
+@@ -18,6 +18,7 @@
  
  # For plotting
  from matplotlib import pyplot
@@ -8,9 +8,9 @@
  
  
  def make_system(a=1, t=1.0, r=10):
-@@ -69,19 +70,24 @@
-         # we only plot the 15 lowest eigenvalues
-         energies.append(ev[:15])
+@@ -70,19 +71,24 @@
+ 
+         energies.append(ev)
  
 -    pyplot.figure()
 +    fig = pyplot.figure()
diff --git a/doc/source/tutorial/closed_system.py b/doc/source/tutorial/closed_system.py
index 5a84896d212956f6bd2709080e04740e4e2d1621..c474b8c5060bc688c7d2fb15bfd0c54b0db2ab51 100644
--- a/doc/source/tutorial/closed_system.py
+++ b/doc/source/tutorial/closed_system.py
@@ -10,11 +10,12 @@
 
 
 from cmath import exp
+import numpy as np
 import kwant
 
 # For eigenvalue computation
 #HIDDEN_BEGIN_tibv
-import scipy.linalg as la
+import scipy.sparse.linalg as sla
 #HIDDEN_END_tibv
 
 # For plotting
@@ -67,12 +68,12 @@ def plot_spectrum(sys, Bfields):
         B = Bfield
 
         # Obtain the Hamiltonian as a dense matrix
-        ham_mat = sys.hamiltonian_submatrix()
+        ham_mat = sys.hamiltonian_submatrix(sparse=True)
 
-        ev = la.eigh(ham_mat, eigvals_only=True)
+        # we only calculate the 15 lowest eigenvalues
+        ev = sla.eigsh(ham_mat, k=15, which='SM', return_eigenvectors=False)
 
-        # we only plot the 15 lowest eigenvalues
-        energies.append(ev[:15])
+        energies.append(ev)
 
     pyplot.figure()
     pyplot.plot(Bfields, energies)
diff --git a/doc/source/tutorial/tutorial3.rst b/doc/source/tutorial/tutorial3.rst
index 06fc4df2f0265ca7f431379e627fed582328dade..f907fba80a37e35516e3d210e62cf705d08fb9f1 100644
--- a/doc/source/tutorial/tutorial3.rst
+++ b/doc/source/tutorial/tutorial3.rst
@@ -56,16 +56,17 @@ Hamiltonian is approximating.
 Closed systems
 ..............
 
-Although kwant is (currently) mainly aimed towards transport problem, it
+Although kwant is (currently) mainly aimed towards transport problema, it
 can also easily be used to compute properties of closed systems -- after
 all, a closed system is nothing more than a scattering region without leads!
 
-In this example, we compute the spectrum of a closed, (approximately)
-circular quantum dot as a function of magnetic field
-(Fock-Darwin spectrum).
+In this example, we compute the wave functions of a closed, (approximately)
+circular quantum dot and its spectrum as a function
+of magnetic field (Fock-Darwin spectrum).
 
-To compute the eigenenergies, we will make use of the linear algebra
-functionality of `scipy <www.scipy.org>`_:
+To compute the eigenenergies and eigenstates, we will make use of the sparse
+linear algebra functionality of `scipy <www.scipy.org>`_, which interfaces
+the ARPACK package:
 
 .. literalinclude:: closed_system.py
     :start-after: #HIDDEN_BEGIN_tibv
@@ -90,9 +91,8 @@ system using `~kwant.system.System.hamiltonian_submatrix`:
     :start-after: #HIDDEN_BEGIN_yvri
     :end-before: #HIDDEN_END_yvri
 
-In this toy model we use dense matrices and dense matrix algebra since
-the system is very small. (In a real application one would probably
-want to use sparse matrix methods.) Finally, we obtain the result:
+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.*