diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index b2b491e26378ab460960a5c28f2c2f7ec91ea535..13db2d37e5f7356bc1d7ef6602840c3fa3bc5975 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -24,11 +24,8 @@ build and upload the contents:
 
   script:
     # Compile lectures
-    - python code/band_structures.py
     - python code/bands_2d.py
-    - python code/misc.py
-    - python code/semiconductors.py
-    - python code/heat_capacity.py
+    - python execute.py
     - mkdocs build
     - "rsync -rv site/* solidstate@tnw-tn1.tudelft.net:"
   only:
diff --git a/code/band_structures.py b/code/band_structures.py
deleted file mode 100644
index 7c4e380b883d17e240259705d8759d99fc0972fd..0000000000000000000000000000000000000000
--- a/code/band_structures.py
+++ /dev/null
@@ -1,227 +0,0 @@
-import os
-
-import numpy as np
-import matplotlib
-import mpl_toolkits
-matplotlib.use('agg')
-from matplotlib import pyplot
-
-import common
-from common import draw_classic_axes
-
-pi = np.pi
-
-
-def fermi_dirac():
-    fig = pyplot.figure()
-    ax = fig.add_subplot(1,1,1)
-    xvals = np.linspace(0, 2, 200)
-    mu = .75
-    beta = 20
-    ax.plot(xvals, xvals < mu, ls='dashed', label='$T=0$')
-    ax.plot(xvals, 1/(np.exp(beta * (xvals-mu)) + 1),
-            ls='solid', label='$T>0$')
-    ax.set_xlabel(r'$\varepsilon$')
-    ax.set_ylabel(r'$f(\varepsilon, T)$')
-    ax.set_yticks([0, 1])
-    ax.set_yticklabels(['$0$', '$1$'])
-    ax.set_xticks([mu])
-    ax.set_xticklabels([r'$\mu$'])
-    ax.set_ylim(-.1, 1.1)
-    ax.legend()
-    draw_classic_axes(ax)
-    pyplot.tight_layout()
-    fig.savefig('fermi_dirac.svg')
-
-
-def phonons_1d():
-    k = np.linspace(-2*pi, 6*pi, 500)
-    fig, ax = pyplot.subplots()
-
-    pyplot.plot(k, np.abs(np.sin(k/2)))
-
-    ax.set_ylim(bottom=0, top=1.2)
-    ax.set_xlabel('$ka$')
-    ax.set_ylabel(r'$\omega$')
-    pyplot.xticks(list(2*pi*np.arange(-1, 4)) + [-pi, pi],
-                  [r'$-2\pi$', '$0$', r'$2\pi$', r'$4\pi$', r'$6\pi$',
-                   r'$-\pi$', r'$\pi$'])
-    pyplot.yticks([1], [r'$2\sqrt{\frac{\kappa}{m}}$'])
-    pyplot.vlines([-pi, pi], 0, 1.1, linestyles='dashed')
-    pyplot.hlines([1], .1*pi, 1.3*pi, linestyles='dashed')
-    draw_classic_axes(ax)
-    ax.annotate(s='', xy=(-pi, -.15), xytext=(pi, -.15),
-                arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
-    ax.text(0, -.25, '1st Brillouin zone', ha='center')
-    ax.set_ylim(bottom=-.7)
-    pyplot.savefig('phonons3.svg')
-
-
-def aliasing():
-    x = np.linspace(-.2, 2.8, 500)
-    fig, ax = pyplot.subplots()
-    ax.plot(x, np.cos(pi*(x - .3)), label=r'$k=\pi/a$')
-    ax.plot(x, np.cos(3*pi*(x - .3)), label=r'$k=3\pi/a$')
-    ax.plot(x, np.cos(5*pi*(x - .3)), label=r'$k=5\pi/a$')
-    sites = np.arange(3) + .3
-    ax.scatter(sites, np.cos(pi*(sites - .3)), c='k', s=64, zorder=5)
-    ax.set_xlabel('$x$')
-    ax.set_ylabel('$u_n$')
-    ax.set_xlim((-.1, 3.2))
-    ax.set_ylim((-1.3, 1.3))
-    ax.legend(loc='lower right')
-    draw_classic_axes(ax)
-    ax.annotate(s='', xy=(.3, -1.1), xytext=(1.3, -1.1),
-                arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
-    ax.text(.3 + .5, -1.25, '$a$', ha='center')
-    pyplot.savefig('phonons4.svg')
-
-
-def tight_binding_1d():
-    pyplot.figure()
-    k = np.linspace(-pi, pi, 300)
-    pyplot.plot(k, -np.cos(k))
-    pyplot.xlabel('$ka$'); pyplot.ylabel('$E$')
-    pyplot.xticks([-pi, 0, pi], [r'$-\pi$', 0, r'$\pi$'])
-    pyplot.yticks([-1, 0, 1], ['$E_0+2t$', '$E_0$', '$E_0-2t$'])
-    pyplot.savefig('tight_binding_1d.svg')
-
-
-def meff_1d_tb():
-    pyplot.figure(figsize=(8, 5))
-    k = np.linspace(-pi, pi, 300)
-    meff = 1/np.cos(k)
-    color = list(matplotlib.rcParams['axes.prop_cycle'])[0]['color']
-    pyplot.plot(k[meff > 0], meff[meff > 0], c=color)
-    pyplot.plot(k[meff < 0], meff[meff < 0], c=color)
-    pyplot.ylim(-5, 5)
-    pyplot.xlabel('$ka$'); pyplot.ylabel('$m_{eff}$')
-    pyplot.xticks([-pi, 0, pi], [r'$-\pi$', 0, r'$\pi$']);
-    pyplot.savefig('meff_1d_tb.svg')
-
-
-def dispersion_2m(k, kappa=1, M=1.4, m=1, acoustic=True):
-    Mm = M*m
-    m_harm = (M + m) / Mm
-    root = kappa * np.sqrt(m_harm**2 - 4*np.sin(k/2)**2 / Mm)
-    if acoustic:
-        root *= -1
-    return np.sqrt(kappa*m_harm + root)
-
-
-def phonons_1d_2masses():
-    # TODO: Add panels with eigenvectors
-    k = np.linspace(-2*pi, 6*pi, 300)
-    fig, ax = pyplot.subplots()
-    ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
-    ax.plot(k, dispersion_2m(k), label='acoustic')
-    ax.set_xlabel('$ka$')
-    ax.set_ylabel(r'$\omega$')
-    pyplot.xticks([-pi, 0, pi], [r'$-\pi$', '$0$', r'$\pi$'])
-    pyplot.yticks([], [])
-    pyplot.vlines([-pi, pi], 0, 2.2, linestyles='dashed')
-    pyplot.legend()
-    pyplot.xlim(-1.75*pi, 3.5*pi)
-    pyplot.ylim(bottom=0)
-    draw_classic_axes(ax, xlabeloffset=.2)
-    pyplot.savefig('phonons6.svg')
-
-
-def DOS_2masses():
-    matplotlib.rcParams['font.size'] = 24
-    k = np.linspace(-.25*pi, 1.5*pi, 300)
-    k_dos = np.linspace(0, pi, 20)
-    fig, (ax, ax2) = pyplot.subplots(ncols=2, sharey=True, figsize=(10, 5))
-    ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
-    ax.plot(k, dispersion_2m(k), label='acoustic')
-    ax.vlines(k_dos, 0, dispersion_2m(k_dos, acoustic=False),
-              colors=(0.5, 0.5, 0.5, 0.5))
-    ax.hlines(
-        np.hstack((dispersion_2m(k_dos, acoustic=False), dispersion_2m(k_dos))), 
-        np.hstack((k_dos, k_dos)), 
-        1.8*pi,
-        colors=(0.5, 0.5, 0.5, 0.5)
-    )
-    ax.set_xlabel('$ka$')
-    ax.set_ylabel(r'$\omega$')
-    ax.set_xticks([0, pi])
-    ax.set_xticklabels(['$0$', r'$\pi$'])
-    ax.set_yticks([], [])
-    ax.set_xlim(-pi/4, 2*pi)
-    ax.set_ylim((0, dispersion_2m(0, acoustic=False) + .2))
-    draw_classic_axes(ax, xlabeloffset=.2)
-    
-    k = np.linspace(0, pi, 1000)
-    omegas = np.hstack((
-        dispersion_2m(k, acoustic=False), dispersion_2m(k)
-    ))
-    ax2.hist(omegas, orientation='horizontal', bins=75)
-    ax2.set_xlabel(r'$g(\omega)$')
-    ax2.set_ylabel(r'$\omega$')
-    
-    # Truncate the singularity in the DOS
-    max_x = ax2.get_xlim()[1]
-    ax2.set_xlim((0, max_x/2))
-    draw_classic_axes(ax2, xlabeloffset=.1)
-    matplotlib.rcParams['font.size'] = 16
-    pyplot.savefig('phonons8.svg')
-
-
-def consistency_1_2_masses():
-    k = np.linspace(0, 2*pi, 300)
-    k_dos = np.linspace(0, pi, 20)
-    fig, ax = pyplot.subplots()
-    ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
-    ax.plot(k, dispersion_2m(k), label='acoustic')
-    omega_max = dispersion_2m(0, acoustic=False)
-    ax.plot(k, omega_max * np.sin(k/4), label='equal masses')
-    ax.set_xlabel('$ka$')
-    ax.set_ylabel(r'$\omega$')
-    ax.set_xticks([0, pi, 2*pi])
-    ax.set_xticklabels(['$0$', r'$\pi/2a$', r'$\pi/a$'])
-    ax.set_yticks([], [])
-    ax.set_xlim(-pi/8, 2*pi+.4)
-    ax.set_ylim((0, dispersion_2m(0, acoustic=False) + .2))
-    ax.legend(loc='lower right')
-    pyplot.vlines([pi, 2*pi], 0, 2.2, linestyles='dashed')
-    draw_classic_axes(ax, xlabeloffset=.2)
-    pyplot.savefig('phonons7.svg')
-
-
-def DOS_finite_phonon_chain(n, output_name):
-    rhs = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
-    rhs[0, 0] -= 1
-    rhs[-1, -1] -= 1
-    pyplot.figure()
-    pyplot.hist(np.sqrt(np.abs(np.linalg.eigvalsh(rhs))), bins=30)
-    pyplot.xlabel("$\omega$")
-    pyplot.ylabel("Number of levels")
-    pyplot.savefig(output_name)
-
-
-def DOS_finite_electron_chain(n, output_name):
-    rhs = - np.eye(n, k=1) - np.eye(n, k=-1)
-    pyplot.figure()
-    pyplot.hist(np.linalg.eigvalsh(rhs), bins=30)
-    pyplot.xlabel("$E$")
-    pyplot.ylabel("Number of levels")
-    pyplot.savefig(output_name)
-
-
-def main():
-    os.chdir('docs/figures')
-    fermi_dirac()
-    phonons_1d()
-    phonons_1d_2masses()
-    consistency_1_2_masses()
-    DOS_2masses()
-    aliasing()
-    meff_1d_tb()
-    tight_binding_1d()
-    DOS_finite_phonon_chain(3, '3_phonons.svg')
-    DOS_finite_phonon_chain(300, '300_phonons.svg')
-    DOS_finite_electron_chain(3, '3_electrons.svg')
-    DOS_finite_electron_chain(300, '300_electrons.svg')
-
-if __name__ == "__main__":
-    main()
diff --git a/code/common.py b/code/common.py
index e51a0987838dbf035ce642d660a580a0e4d92824..df2d1fb38b6cfb89106cd66930c8a89eef5df0e6 100644
--- a/code/common.py
+++ b/code/common.py
@@ -1,22 +1,33 @@
-import os
-
-import numpy as np
+from IPython import get_ipython
 import matplotlib
-import mpl_toolkits
 
-matplotlib.rcParams['text.usetex'] = True
-matplotlib.rcParams['figure.figsize'] = (8, 5)
-matplotlib.rcParams['font.size'] = 16
-matplotlib.rcParams['savefig.transparent'] = True
+
+def configure_plotting():
+    # Configure matplotlib
+    ip = get_ipython()
+    ip.config['InlineBackend']['figure_format'] = 'svg'
+
+    # Fix for https://stackoverflow.com/a/36021869/2217463
+    ip.config['InlineBackend']['rc'] = {}
+    ip.enable_matplotlib(gui='inline')
+
+    matplotlib.rcParams['text.usetex'] = False
+    matplotlib.rcParams['figure.figsize'] = (10, 7)
+    matplotlib.rcParams['font.size'] = 16
+
 
 def draw_classic_axes(ax, x=0, y=0, xlabeloffset=.1, ylabeloffset=.07):
     ax.set_axis_off()
     x0, x1 = ax.get_xlim()
     y0, y1 = ax.get_ylim()
-    ax.annotate(ax.get_xlabel(), xytext=(x1, y), xy=(x0, y),
-            arrowprops=dict(arrowstyle="<-"), va='center')
-    ax.annotate(ax.get_ylabel(), xytext=(x, y1), xy=(x, y0),
-            arrowprops=dict(arrowstyle="<-"), ha='center')
+    ax.annotate(
+        ax.get_xlabel(), xytext=(x1, y), xy=(x0, y),
+        arrowprops=dict(arrowstyle="<-"), va='center'
+    )
+    ax.annotate(
+        ax.get_ylabel(), xytext=(x, y1), xy=(x, y0),
+        arrowprops=dict(arrowstyle="<-"), ha='center'
+    )
     for pos, label in zip(ax.get_xticks(), ax.get_xticklabels()):
         ax.text(pos, y - xlabeloffset, label.get_text(),
                 ha='center', va='bottom')
diff --git a/code/heat_capacity.py b/code/heat_capacity.py
deleted file mode 100644
index 41282ed3616f68b6f617ad077efb4c216d2c65ee..0000000000000000000000000000000000000000
--- a/code/heat_capacity.py
+++ /dev/null
@@ -1,106 +0,0 @@
-import os
-
-import matplotlib
-matplotlib.use('agg')
-from matplotlib import pyplot
-
-import numpy as np
-from scipy.optimize import curve_fit
-from scipy.integrate import quad
-
-import common
-from common import draw_classic_axes
-
-
-# Data from Einstein's paper
-T = [222.4, 262.4, 283.7, 306.4, 331.3, 358.5, 413.0, 479.2, 520.0, 879.7, 1079.7, 1258.0],
-c = [0.384, 0.578, 0.683, 0.798, 0.928, 1.069, 1.343, 1.656, 1.833, 2.671, 2.720, 2.781]
-
-def plot_n_e():
-    fig, (ax, ax2) = pyplot.subplots(ncols=2, figsize=(10, 5))
-    omega = np.linspace(0.1, 2)
-    ax.plot(omega, 1/(np.exp(omega) - 1))
-    ax.set_ylim(0, top=3)
-    ax.set_xlim(left=0)
-    ax.set_xlabel(r'$\hbar \omega$')
-    ax.set_xticks([0, 1])
-    ax.set_xticklabels(['$0$', '$k_B T$'])
-    ax.set_ylabel('$n$')
-    ax.set_yticks([1, 2])
-    ax.set_yticklabels(['$1$', '$2$'])
-    draw_classic_axes(ax, xlabeloffset=.2)
-    T = np.linspace(0.01, 2)
-    ax2.plot(T, 1/2 + 1/(np.exp(1/T)-1))
-    ax2.set_ylim(bottom=0)
-    ax2.set_xlabel('$k_B T$')
-    ax2.set_xticks([0, 1])
-    ax2.set_xticklabels(['$0$', r'$\hbar \omega$'])
-    ax2.set_ylabel(r"$\bar\varepsilon$")
-    ax2.set_yticks([1/2])
-    ax2.set_yticklabels([r'$\hbar\omega/2$'])
-    draw_classic_axes(ax2, xlabeloffset=.15)
-    fig.savefig('bose_einstein.svg')
-
-
-def c_einstein(T, T_E=1):
-    x = T_E / T
-    return 3 * x**2 * np.exp(x) / (np.exp(x) - 1)**2
-
-
-def integrand(y):
-    return y**4 * np.exp(y) / (np.exp(y) - 1)**2
-
-
-@np.vectorize
-def c_debye(T, T_D=1):
-    x = T / T_D
-    return 9 * x**3 * quad(integrand, 0, 1/x)[0]
-
-
-def plot_einstein():
-    T = np.linspace(0.01, 1.5, 500)
-    fig, ax = pyplot.subplots()
-
-    ax.plot(T, c_einstein(T)/3)
-    ax.fill_between(T, c_einstein(T)/3, 1, alpha=0.5)
-
-    ax.set_ylim(bottom=0, top=1.2)
-    ax.set_xlabel('$T$')
-    ax.set_ylabel(r'$\omega$')
-    ax.set_xticks([1])
-    ax.set_xticklabels([r'$\hbar \omega/k_B$'])
-    ax.set_yticks([1])
-    ax.set_yticklabels(['$k_B$'])
-    pyplot.hlines([1], 0, 1.5, linestyles='dashed')
-    draw_classic_axes(ax)
-    pyplot.savefig('zeropoint.svg')
-
-
-def plot_debye():
-    T = np.linspace(0.01, 1.5, 500)
-    fig, ax = pyplot.subplots()
-
-    ax.plot(T, c_einstein(T), label="Einstein model")
-    ax.plot(T, c_debye(T), label="Debye model")
-
-    ax.set_ylim(bottom=0, top=3.4)
-    ax.set_xlabel('$T$')
-    ax.set_ylabel(r'$\omega$')
-    ax.set_xticks([1])
-    ax.set_xticklabels([r'$\Theta_D$'])
-    ax.set_yticks([3])
-    ax.set_yticklabels(['$3Nk_B$'])
-    ax.legend(loc='lower right')
-    pyplot.hlines([3], 0, 1.5, linestyles='dashed')
-    draw_classic_axes(ax, xlabeloffset=0.3)
-    pyplot.savefig('debye.svg')
-
-
-def main():
-    os.chdir('docs/figures')
-    plot_n_e()
-    plot_einstein()
-    plot_debye()
-
-if __name__ == "__main__":
-    main()
diff --git a/code/misc.py b/code/misc.py
deleted file mode 100644
index 28345843e0ccb9d37c1dbc107f96e774cd6d8cd5..0000000000000000000000000000000000000000
--- a/code/misc.py
+++ /dev/null
@@ -1,75 +0,0 @@
-import os
-
-import matplotlib
-matplotlib.use('agg')
-from matplotlib import pyplot
-
-import numpy as np
-from scipy.optimize import curve_fit
-from scipy.integrate import quad
-
-import common
-from common import draw_classic_axes
-
-def plot_curie():    
-    B = np.linspace(-2, 2, 1000)
-    fig, ax = pyplot.subplots()
-
-    sqrt_plus = lambda x: np.sqrt(x * (x >= 0))
-    ax.plot(B, np.tanh(B * 3), label="low $T$")
-    ax.plot(B, np.tanh(B), label="high $T$")
-
-    ax.legend()
-    ax.set_ylabel("$M$")
-    ax.set_xlabel("$B$")
-    draw_classic_axes(ax, xlabeloffset=.2)
-    fig.savefig('curie.svg')
-
-
-def plot_free_electron_DOS():    
-    E = np.linspace(0, 2, 500)
-    fig, ax = pyplot.subplots()
-
-    ax.plot(E, np.sqrt(E))
-
-    ax.set_ylabel(r"$g(\varepsilon)$")
-    ax.set_xlabel(r"$\varepsilon$")
-    draw_classic_axes(ax, xlabeloffset=.2)
-    fig.savefig('sqrt.svg')
-
-
-def plot_finite_T_FEM():
-    E = np.linspace(0, 2, 500)
-    fig, ax = pyplot.subplots()
-    ax.plot(E, np.sqrt(E), linestyle='dashed')
-    ax.text(1.7, 1.4, r'$g(\varepsilon)\propto \sqrt{\varepsilon}$', ha='center')
-    ax.fill_between(E, np.sqrt(E) * (E < 1), alpha=.3)
-    
-    n = np.sqrt(E) / (1 + np.exp(20*(E-1)))
-    ax.plot(E, n)
-    ax.fill_between(E, n, alpha=.5)
-    w = .17
-    ax.annotate(s='', xy=(1, 1), xytext=(1-w, 1),
-                arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
-    ax.text(1-w/2, 1.1, r'$\sim k_BT$', ha='center')
-    ax.plot([1-w, 1+w], [1, 0], c='k', linestyle='dashed')
-    ax.annotate(s='', xy=(1, 0), xytext=(1, 1),
-                arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
-    ax.text(1.2, .7, r'$g(\varepsilon_F)$', ha='center')
-    ax.set_xticks([1])
-    ax.set_xticklabels([r'$\varepsilon_F$'])
-
-    ax.set_ylabel(r"$g(\varepsilon)$")
-    ax.set_xlabel(r"$\varepsilon$")
-    draw_classic_axes(ax, xlabeloffset=.2)
-    
-    fig.savefig('FD_DOS.svg')
-
-def main():
-    os.chdir('docs/figures')
-    plot_curie()
-    plot_free_electron_DOS()
-    plot_finite_T_FEM()
-
-if __name__ == "__main__":
-    main()
diff --git a/code/semiconductors.py b/code/semiconductors.py
deleted file mode 100644
index 3d578281e29824a83db77f56093b8f28a28bd733..0000000000000000000000000000000000000000
--- a/code/semiconductors.py
+++ /dev/null
@@ -1,70 +0,0 @@
-import os
-
-import matplotlib
-matplotlib.use('agg')
-from matplotlib import pyplot
-
-import numpy as np
-from scipy.optimize import curve_fit
-from scipy.integrate import quad
-
-import common
-from common import draw_classic_axes
-
-E_V, E_C, E_F = -1.2, 1.8, .4
-E_D, E_A = E_C - .7, E_V + .5
-m_h, m_e = 1, .5
-sqrt_plus = lambda x: np.sqrt(x * (x >= 0))
-
-
-def plot_dos():
-    E = np.linspace(-3, 3, 1000)
-    fig, ax = pyplot.subplots()
-
-    n_F = 1/(np.exp(2*(E - E_F)) + 1)
-    g_e = m_e * sqrt_plus(E - E_C)
-    g_h = m_h * sqrt_plus(E_V - E)
-    ax.plot(E, g_h, label="$g_e$")
-    ax.plot(E, g_e, label="$g_h$")
-    ax.plot(E, 10 * g_h * (1-n_F), label="$n_h$")
-    ax.plot(E, 10 * g_e * n_F, label="$n_e$")
-    ax.plot(E, n_F, label="$n_F$", linestyle='dashed')
-    ax.set_ylim(top=1.5)
-
-    ax.set_xlabel('$E$')
-    ax.set_ylabel('$g$')
-    ax.set_xticks([E_V, E_C, E_F])
-    ax.set_xticklabels(['$E_V$', '$E_C$', '$E_F$'])
-    ax.legend()
-    draw_classic_axes(ax, xlabeloffset=.2)
-    fig.savefig('intrinsic_DOS.svg')
-
-def plot_doping():    
-    E = np.linspace(-3, 3, 1000)
-    fig, ax = pyplot.subplots()
-
-    n_F = 1/(np.exp(2*(E - E_F)) + 1)
-    g_e = m_e * sqrt_plus(E - E_C)
-    g_h = m_h * sqrt_plus(E_V - E)
-    ax.plot(E, g_h, label="$g_e$")
-    ax.plot(E, g_e, label="$g_h$")
-
-    sigma = 0.01
-    g_D = np.exp(-(E_D - E)**2 / sigma**2)
-    g_A = .7 * np.exp(-(E_A - E)**2 / sigma**2)
-    ax.plot(E, g_D, label='$g_D$')
-    ax.plot(E, g_A, label='$g_A$')
-    ax.legend()
-    ax.set_xticks([E_V, E_C, E_A, E_D])
-    ax.set_xticklabels(['$E_V$', '$E_C$', '$E_A$', '$E_D$'])
-    draw_classic_axes(ax, xlabeloffset=.2)
-    fig.savefig('doping.svg')
-
-
-def main():
-    os.chdir('docs/figures')
-    plot_dos()
-    plot_doping()
-
-if __name__ == "__main__":
-    main()
diff --git a/docs/mathjaxhelper.js b/docs/mathjaxhelper.js
deleted file mode 100644
index bec7b177c7c1d2ca063f416db4af27fed5e6b824..0000000000000000000000000000000000000000
--- a/docs/mathjaxhelper.js
+++ /dev/null
@@ -1,5 +0,0 @@
-MathJax.Hub.Config({
-  config: ["MMLorHTML.js"],
-  jax: ["input/TeX", "output/HTML-CSS", "output/NativeMML"],
-  extensions: ["MathMenu.js", "MathZoom.js"]
-});
diff --git a/execute.py b/execute.py
new file mode 100644
index 0000000000000000000000000000000000000000..71c7b23ea2d2bffa62a0ec3efd43f7a936c2ccdc
--- /dev/null
+++ b/execute.py
@@ -0,0 +1,43 @@
+from pathlib import Path
+import shutil
+
+import nbconvert
+import notedown
+from traitlets.config import Config
+
+reader = notedown.MarkdownReader()
+
+src = Path('src')
+target = Path('docs')
+shutil.rmtree(target, ignore_errors=True)
+target.mkdir(exist_ok=True)
+shutil.copytree(src / 'figures', target / 'figures')
+
+exporter = nbconvert.MarkdownExporter(
+    config=Config(dict(
+        MarkdownExporter=dict(
+            preprocessors=[
+                nbconvert.preprocessors.ExecutePreprocessor,
+                nbconvert.preprocessors.ExtractOutputPreprocessor,
+            ],
+            exclude_input=True
+        )
+    ))
+)
+
+writer = nbconvert.writers.FilesWriter(build_directory=str(target))
+
+for source_file in src.glob('*.md'):
+    fname = source_file.name[:-len('.md')]
+    notebook = reader.reads(source_file.read_text())
+
+    output, resources = exporter.from_notebook_node(
+        notebook,
+        resources={
+            'unique_key': fname,
+            'output_files_dir': 'figures',
+            'metadata': {'path': str(Path('.') / 'code')}
+        }
+    )
+
+    writer.write(output, resources, fname)
diff --git a/mkdocs.yml b/mkdocs.yml
index e55b2c3abb03e4bc6239d4cd541d1b6cdc542a5e..a6b4cc1fdb778ab82464b653b62670c8052f81b1 100644
--- a/mkdocs.yml
+++ b/mkdocs.yml
@@ -13,7 +13,8 @@ pages:
   - Crystal structure and diffraction: 'lecture_5.md'
   - Tight binding and nearly free electrons: 'lecture_6.md'
   - Semiconductors: 'lecture_7.md'
-  - Magnetism: 'lecture_8.md'
+  - Attic:
+    - Magnetism: 'lecture_8.md'
 
 theme:
   name: material
@@ -22,7 +23,6 @@ theme:
     primary: 'white'
     accent: 'indigo'
 
-
 markdown_extensions:
   - mdx_math:
       enable_dollar_delimiter: True
@@ -33,6 +33,5 @@ markdown_extensions:
 
 extra_javascript:
   - 'https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.3/MathJax.js?config=TeX-AMS_HTML'
-  - 'mathjaxhelper.js'
 
 copyright: "Copyright © 2017-2018 Delft University of Technology, CC-BY-SA-NC 4.0."
diff --git a/docs/figures/DOS.svg b/src/figures/DOS.svg
similarity index 100%
rename from docs/figures/DOS.svg
rename to src/figures/DOS.svg
diff --git a/docs/figures/DOS_periodic.svg b/src/figures/DOS_periodic.svg
similarity index 100%
rename from docs/figures/DOS_periodic.svg
rename to src/figures/DOS_periodic.svg
diff --git a/docs/figures/E_F_and_carrier_density.svg b/src/figures/E_F_and_carrier_density.svg
similarity index 100%
rename from docs/figures/E_F_and_carrier_density.svg
rename to src/figures/E_F_and_carrier_density.svg
diff --git a/docs/figures/RKKY.svg b/src/figures/RKKY.svg
similarity index 100%
rename from docs/figures/RKKY.svg
rename to src/figures/RKKY.svg
diff --git a/docs/figures/adsorption.svg b/src/figures/adsorption.svg
similarity index 100%
rename from docs/figures/adsorption.svg
rename to src/figures/adsorption.svg
diff --git a/docs/figures/band_diagram_question.svg b/src/figures/band_diagram_question.svg
similarity index 100%
rename from docs/figures/band_diagram_question.svg
rename to src/figures/band_diagram_question.svg
diff --git a/docs/figures/band_diagram_solution.svg b/src/figures/band_diagram_solution.svg
similarity index 100%
rename from docs/figures/band_diagram_solution.svg
rename to src/figures/band_diagram_solution.svg
diff --git a/docs/figures/band_structure_sketch.svg b/src/figures/band_structure_sketch.svg
similarity index 100%
rename from docs/figures/band_structure_sketch.svg
rename to src/figures/band_structure_sketch.svg
diff --git a/docs/figures/bonding.svg b/src/figures/bonding.svg
similarity index 100%
rename from docs/figures/bonding.svg
rename to src/figures/bonding.svg
diff --git a/docs/figures/bonding_with_repulsion.svg b/src/figures/bonding_with_repulsion.svg
similarity index 100%
rename from docs/figures/bonding_with_repulsion.svg
rename to src/figures/bonding_with_repulsion.svg
diff --git a/docs/figures/bragg.svg b/src/figures/bragg.svg
similarity index 100%
rename from docs/figures/bragg.svg
rename to src/figures/bragg.svg
diff --git a/docs/figures/bravais1.svg b/src/figures/bravais1.svg
similarity index 100%
rename from docs/figures/bravais1.svg
rename to src/figures/bravais1.svg
diff --git a/docs/figures/bravais2.svg b/src/figures/bravais2.svg
similarity index 100%
rename from docs/figures/bravais2.svg
rename to src/figures/bravais2.svg
diff --git a/docs/figures/brillouin.svg b/src/figures/brillouin.svg
similarity index 100%
rename from docs/figures/brillouin.svg
rename to src/figures/brillouin.svg
diff --git a/docs/figures/brillouin_mod.svg b/src/figures/brillouin_mod.svg
similarity index 100%
rename from docs/figures/brillouin_mod.svg
rename to src/figures/brillouin_mod.svg
diff --git a/docs/figures/crystalfield.svg b/src/figures/crystalfield.svg
similarity index 100%
rename from docs/figures/crystalfield.svg
rename to src/figures/crystalfield.svg
diff --git a/docs/figures/cubic.svg b/src/figures/cubic.svg
similarity index 100%
rename from docs/figures/cubic.svg
rename to src/figures/cubic.svg
diff --git a/docs/figures/cubic_mod.svg b/src/figures/cubic_mod.svg
similarity index 100%
rename from docs/figures/cubic_mod.svg
rename to src/figures/cubic_mod.svg
diff --git a/docs/figures/direct_indirect.svg b/src/figures/direct_indirect.svg
similarity index 100%
rename from docs/figures/direct_indirect.svg
rename to src/figures/direct_indirect.svg
diff --git a/docs/figures/dos_1d_tb.svg b/src/figures/dos_1d_tb.svg
similarity index 100%
rename from docs/figures/dos_1d_tb.svg
rename to src/figures/dos_1d_tb.svg
diff --git a/docs/figures/einstein.svg b/src/figures/einstein.svg
similarity index 100%
rename from docs/figures/einstein.svg
rename to src/figures/einstein.svg
diff --git a/docs/figures/electric_field.svg b/src/figures/electric_field.svg
similarity index 100%
rename from docs/figures/electric_field.svg
rename to src/figures/electric_field.svg
diff --git a/docs/figures/electrons.svg b/src/figures/electrons.svg
similarity index 100%
rename from docs/figures/electrons.svg
rename to src/figures/electrons.svg
diff --git a/docs/figures/ewald.svg b/src/figures/ewald.svg
similarity index 100%
rename from docs/figures/ewald.svg
rename to src/figures/ewald.svg
diff --git a/docs/figures/extended_nearly_free_electron_bands.svg b/src/figures/extended_nearly_free_electron_bands.svg
similarity index 100%
rename from docs/figures/extended_nearly_free_electron_bands.svg
rename to src/figures/extended_nearly_free_electron_bands.svg
diff --git a/docs/figures/fcc.svg b/src/figures/fcc.svg
similarity index 100%
rename from docs/figures/fcc.svg
rename to src/figures/fcc.svg
diff --git a/docs/figures/fcc_mod.svg b/src/figures/fcc_mod.svg
similarity index 100%
rename from docs/figures/fcc_mod.svg
rename to src/figures/fcc_mod.svg
diff --git a/docs/figures/fermi_circle.svg b/src/figures/fermi_circle.svg
similarity index 100%
rename from docs/figures/fermi_circle.svg
rename to src/figures/fermi_circle.svg
diff --git a/docs/figures/fermi_circle_periodic.svg b/src/figures/fermi_circle_periodic.svg
similarity index 100%
rename from docs/figures/fermi_circle_periodic.svg
rename to src/figures/fermi_circle_periodic.svg
diff --git a/docs/figures/fermi_level.svg b/src/figures/fermi_level.svg
similarity index 100%
rename from docs/figures/fermi_level.svg
rename to src/figures/fermi_level.svg
diff --git a/docs/figures/free_electron.svg b/src/figures/free_electron.svg
similarity index 100%
rename from docs/figures/free_electron.svg
rename to src/figures/free_electron.svg
diff --git a/docs/figures/graphite.svg b/src/figures/graphite.svg
similarity index 100%
rename from docs/figures/graphite.svg
rename to src/figures/graphite.svg
diff --git a/docs/figures/graphite_mod.svg b/src/figures/graphite_mod.svg
similarity index 100%
rename from docs/figures/graphite_mod.svg
rename to src/figures/graphite_mod.svg
diff --git a/docs/figures/hall_effect.svg b/src/figures/hall_effect.svg
similarity index 100%
rename from docs/figures/hall_effect.svg
rename to src/figures/hall_effect.svg
diff --git a/docs/figures/harmonic.svg b/src/figures/harmonic.svg
similarity index 100%
rename from docs/figures/harmonic.svg
rename to src/figures/harmonic.svg
diff --git a/docs/figures/holes.svg b/src/figures/holes.svg
similarity index 100%
rename from docs/figures/holes.svg
rename to src/figures/holes.svg
diff --git a/docs/figures/hund.svg b/src/figures/hund.svg
similarity index 100%
rename from docs/figures/hund.svg
rename to src/figures/hund.svg
diff --git a/docs/figures/interatomic_interaction.svg b/src/figures/interatomic_interaction.svg
similarity index 100%
rename from docs/figures/interatomic_interaction.svg
rename to src/figures/interatomic_interaction.svg
diff --git a/docs/figures/lattice_potential.svg b/src/figures/lattice_potential.svg
similarity index 100%
rename from docs/figures/lattice_potential.svg
rename to src/figures/lattice_potential.svg
diff --git a/docs/figures/laue.svg b/src/figures/laue.svg
similarity index 100%
rename from docs/figures/laue.svg
rename to src/figures/laue.svg
diff --git a/docs/figures/laue_mod.svg b/src/figures/laue_mod.svg
similarity index 100%
rename from docs/figures/laue_mod.svg
rename to src/figures/laue_mod.svg
diff --git a/docs/figures/matthiessen.svg b/src/figures/matthiessen.svg
similarity index 100%
rename from docs/figures/matthiessen.svg
rename to src/figures/matthiessen.svg
diff --git a/docs/figures/metal_semiconductor.svg b/src/figures/metal_semiconductor.svg
similarity index 100%
rename from docs/figures/metal_semiconductor.svg
rename to src/figures/metal_semiconductor.svg
diff --git a/docs/figures/miller.svg b/src/figures/miller.svg
similarity index 100%
rename from docs/figures/miller.svg
rename to src/figures/miller.svg
diff --git a/docs/figures/models.svg b/src/figures/models.svg
similarity index 100%
rename from docs/figures/models.svg
rename to src/figures/models.svg
diff --git a/docs/figures/nearly_free_electron_bands.svg b/src/figures/nearly_free_electron_bands.svg
similarity index 100%
rename from docs/figures/nearly_free_electron_bands.svg
rename to src/figures/nearly_free_electron_bands.svg
diff --git a/docs/figures/packing.svg b/src/figures/packing.svg
similarity index 100%
rename from docs/figures/packing.svg
rename to src/figures/packing.svg
diff --git a/docs/figures/phonons1.svg b/src/figures/phonons1.svg
similarity index 100%
rename from docs/figures/phonons1.svg
rename to src/figures/phonons1.svg
diff --git a/docs/figures/phonons2.svg b/src/figures/phonons2.svg
similarity index 100%
rename from docs/figures/phonons2.svg
rename to src/figures/phonons2.svg
diff --git a/docs/figures/phonons5.svg b/src/figures/phonons5.svg
similarity index 100%
rename from docs/figures/phonons5.svg
rename to src/figures/phonons5.svg
diff --git a/docs/figures/pn_charge_density.svg b/src/figures/pn_charge_density.svg
similarity index 100%
rename from docs/figures/pn_charge_density.svg
rename to src/figures/pn_charge_density.svg
diff --git a/docs/figures/pn_junction_bias.svg b/src/figures/pn_junction_bias.svg
similarity index 100%
rename from docs/figures/pn_junction_bias.svg
rename to src/figures/pn_junction_bias.svg
diff --git a/docs/figures/realistic.svg b/src/figures/realistic.svg
similarity index 100%
rename from docs/figures/realistic.svg
rename to src/figures/realistic.svg
diff --git a/docs/figures/reciprocal.svg b/src/figures/reciprocal.svg
similarity index 100%
rename from docs/figures/reciprocal.svg
rename to src/figures/reciprocal.svg
diff --git a/docs/figures/reciprocal_mod.svg b/src/figures/reciprocal_mod.svg
similarity index 100%
rename from docs/figures/reciprocal_mod.svg
rename to src/figures/reciprocal_mod.svg
diff --git a/docs/figures/reduced_nearly_free_electron_bands.svg b/src/figures/reduced_nearly_free_electron_bands.svg
similarity index 100%
rename from docs/figures/reduced_nearly_free_electron_bands.svg
rename to src/figures/reduced_nearly_free_electron_bands.svg
diff --git a/docs/figures/semiconductor.svg b/src/figures/semiconductor.svg
similarity index 100%
rename from docs/figures/semiconductor.svg
rename to src/figures/semiconductor.svg
diff --git a/docs/figures/single_atom.svg b/src/figures/single_atom.svg
similarity index 100%
rename from docs/figures/single_atom.svg
rename to src/figures/single_atom.svg
diff --git a/docs/figures/soft_xray.svg b/src/figures/soft_xray.svg
similarity index 100%
rename from docs/figures/soft_xray.svg
rename to src/figures/soft_xray.svg
diff --git a/docs/figures/solar_cell.svg b/src/figures/solar_cell.svg
similarity index 100%
rename from docs/figures/solar_cell.svg
rename to src/figures/solar_cell.svg
diff --git a/docs/figures/stacking.svg b/src/figures/stacking.svg
similarity index 100%
rename from docs/figures/stacking.svg
rename to src/figures/stacking.svg
diff --git a/docs/figures/superexchange.svg b/src/figures/superexchange.svg
similarity index 100%
rename from docs/figures/superexchange.svg
rename to src/figures/superexchange.svg
diff --git a/docs/figures/transport.svg b/src/figures/transport.svg
similarity index 100%
rename from docs/figures/transport.svg
rename to src/figures/transport.svg
diff --git a/docs/figures/two_atoms.svg b/src/figures/two_atoms.svg
similarity index 100%
rename from docs/figures/two_atoms.svg
rename to src/figures/two_atoms.svg
diff --git a/docs/figures/xray.svg b/src/figures/xray.svg
similarity index 100%
rename from docs/figures/xray.svg
rename to src/figures/xray.svg
diff --git a/docs/index.md b/src/index.md
similarity index 100%
rename from docs/index.md
rename to src/index.md
diff --git a/docs/lecture_1.md b/src/lecture_1.md
similarity index 74%
rename from docs/lecture_1.md
rename to src/lecture_1.md
index 8a43009edcabf1b7008a62a0e870c2f6e7c5c855..babc7992537dbd11bc1109270e1ec28eec12fd71 100644
--- a/docs/lecture_1.md
+++ b/src/lecture_1.md
@@ -1,3 +1,14 @@
+```python
+from matplotlib import pyplot
+
+import numpy as np
+from scipy.optimize import curve_fit
+from scipy.integrate import quad
+
+from common import draw_classic_axes, configure_plotting
+
+configure_plotting()
+```
 
 # Lecture 1 – Phonons and specific Heat
 _Based on chapter 2 of the book_
@@ -25,7 +36,31 @@ $$
 n(\omega,T)=\frac{1}{ {\rm e}^{\hbar\omega/k_{\rm B}T}-1}\Rightarrow\bar{\varepsilon}=\frac{1}{2}\hbar\omega+\frac{\hbar\omega}{ {\rm e}^{\hbar\omega/k_{\rm B}T}-1}
 $$
 
-![](figures/bose_einstein.svg)
+```python
+
+fig, (ax, ax2) = pyplot.subplots(ncols=2, figsize=(10, 5))
+omega = np.linspace(0.1, 2)
+ax.plot(omega, 1/(np.exp(omega) - 1))
+ax.set_ylim(0, top=3)
+ax.set_xlim(left=0)
+ax.set_xlabel(r'$\hbar \omega$')
+ax.set_xticks([0, 1])
+ax.set_xticklabels(['$0$', '$k_B T$'])
+ax.set_ylabel('$n$')
+ax.set_yticks([1, 2])
+ax.set_yticklabels(['$1$', '$2$'])
+draw_classic_axes(ax, xlabeloffset=.2)
+T = np.linspace(0.01, 2)
+ax2.plot(T, 1/2 + 1/(np.exp(1/T)-1))
+ax2.set_ylim(bottom=0)
+ax2.set_xlabel('$k_B T$')
+ax2.set_xticks([0, 1])
+ax2.set_xticklabels(['$0$', r'$\hbar \omega$'])
+ax2.set_ylabel(r"$\bar\varepsilon$")
+ax2.set_yticks([1/2])
+ax2.set_yticklabels([r'$\hbar\omega/2$'])
+draw_classic_axes(ax2, xlabeloffset=.15)
+```
 
 The term $\frac{1}{2}\hbar\omega$ is the _zero point energy_, which follows from the uncertainty principle.
 
@@ -40,7 +75,31 @@ C = \frac{\partial\bar{\varepsilon}}{\partial T}
 \end{multline}
 $$
 
-![](figures/zeropoint.svg)
+```python
+# Data from Einstein's paper
+T = [222.4, 262.4, 283.7, 306.4, 331.3, 358.5, 413.0, 479.2, 520.0, 879.7, 1079.7, 1258.0],
+c = [0.384, 0.578, 0.683, 0.798, 0.928, 1.069, 1.343, 1.656, 1.833, 2.671, 2.720, 2.781]
+
+def c_einstein(T, T_E=1):
+    x = T_E / T
+    return 3 * x**2 * np.exp(x) / (np.exp(x) - 1)**2
+
+T = np.linspace(0.01, 1.5, 500)
+fig, ax = pyplot.subplots()
+
+ax.plot(T, c_einstein(T)/3)
+ax.fill_between(T, c_einstein(T)/3, 1, alpha=0.5)
+
+ax.set_ylim(bottom=0, top=1.2)
+ax.set_xlabel('$T$')
+ax.set_ylabel(r'$\omega$')
+ax.set_xticks([1])
+ax.set_xticklabels([r'$\hbar \omega/k_B$'])
+ax.set_yticks([1])
+ax.set_yticklabels(['$k_B$'])
+pyplot.hlines([1], 0, 1.5, linestyles='dashed')
+draw_classic_axes(ax)
+```
 
 The dashed line signifies the classical value, $k_{\rm B}$. Shaded area $=\frac{1}{2}\hbar\omega$, the zero point energy that cannot be removed through cooling.
 
@@ -126,4 +185,29 @@ $$
 
 where $x=\frac{\hbar\omega}{k_{\rm B}T}$ and $\Theta_{\rm D}\equiv\frac{\hbar\omega_{\rm D}}{k_{\rm B}}$, the _Debye temperature_.
 
-![](figures/debye.svg)
+```python
+def integrand(y):
+    return y**4 * np.exp(y) / (np.exp(y) - 1)**2
+
+@np.vectorize
+def c_debye(T, T_D=1):
+    x = T / T_D
+    return 9 * x**3 * quad(integrand, 0, 1/x)[0]
+
+T = np.linspace(0.01, 1.5, 500)
+fig, ax = pyplot.subplots()
+
+ax.plot(T, c_einstein(T), label="Einstein model")
+ax.plot(T, c_debye(T), label="Debye model")
+
+ax.set_ylim(bottom=0, top=3.4)
+ax.set_xlabel('$T$')
+ax.set_ylabel(r'$\omega$')
+ax.set_xticks([1])
+ax.set_xticklabels([r'$\Theta_D$'])
+ax.set_yticks([3])
+ax.set_yticklabels(['$3Nk_B$'])
+ax.legend(loc='lower right')
+pyplot.hlines([3], 0, 1.5, linestyles='dashed')
+draw_classic_axes(ax, xlabeloffset=0.3)
+```
diff --git a/docs/lecture_2.md b/src/lecture_2.md
similarity index 84%
rename from docs/lecture_2.md
rename to src/lecture_2.md
index 8af060a5ff5b5d62ac8ec8d369c45388657eb521..e4d0bc1692cea27e8a821719f92505f6b1baf27b 100644
--- a/docs/lecture_2.md
+++ b/src/lecture_2.md
@@ -1,3 +1,12 @@
+```python
+from matplotlib import pyplot
+
+import numpy as np
+
+from common import draw_classic_axes, configure_plotting
+
+configure_plotting()
+```
 
 # Lecture 2 – Free electron model
 _Based on chapters 3 and 4 of the book_
@@ -51,7 +60,7 @@ How fast do electrons travel through a copper wire? Let's take $E$ = 1 volt/m, $
 $\rightarrow v=\mu E=\frac{e\tau}{m_{\rm e}}E=\frac{10^{-19}\times 2.5\times 10^{-14}}{10^{-30}}=2.5\times10^{-3}=2.5$ mm/s ! (= 50 $\mu$m @ 50 Hz AC)
 
 ### Hall effect
-Consider a conductive wire in a magnetic field ${\bf B} \rightarrow$ electrons are deflected in a direction perpendicular to ${\bf B}$ and ${\bf j}$. 
+Consider a conductive wire in a magnetic field ${\bf B} \rightarrow$ electrons are deflected in a direction perpendicular to ${\bf B}$ and ${\bf j}$.
 
 ![](figures/hall_effect.svg)
 
@@ -94,7 +103,7 @@ Comparable to phonons, but: electrons are _fermions_.
 
 - Only 2 (due to spin) allowed per $k$-value
 - Fill up from the lowest energy until you run out of electrons
-   
+
 $\rightarrow$ Calculate when you are out of electrons $\rightarrow$ _Fermi energy_.
 
 ![](figures/fermi_circle_periodic.svg)
@@ -117,7 +126,16 @@ $$
 g(\varepsilon)=\frac{ {\rm d}N}{ {\rm d}\varepsilon}=\frac{Vm^{3/2}\sqrt{2\varepsilon}}{\pi^2\hbar^3}\propto\sqrt{\varepsilon}
 $$
 
-![](figures/sqrt.svg)
+```python
+E = np.linspace(0, 2, 500)
+fig, ax = pyplot.subplots()
+
+ax.plot(E, np.sqrt(E))
+
+ax.set_ylabel(r"$g(\varepsilon)$")
+ax.set_xlabel(r"$\varepsilon$")
+draw_classic_axes(ax, xlabeloffset=.2)
+```
 
 Similarly,
 
@@ -133,7 +151,7 @@ $$
 with $\varepsilon_{\rm F}$ the _Fermi energy_ = highest filled energy at $T=0$.
 
 $$
-\varepsilon_{\rm F}=\frac{\hbar^2}{2m}\left( 3\pi^2\frac{N}{V} \right)^{2/3}\equiv \frac{\hbar^2 k_{\rm F}^2}{2m},\ 
+\varepsilon_{\rm F}=\frac{\hbar^2}{2m}\left( 3\pi^2\frac{N}{V} \right)^{2/3}\equiv \frac{\hbar^2 k_{\rm F}^2}{2m},\
 k_{\rm F}=\left( 3\pi^2\frac{N}{V} \right)^{1/3}
 $$
 
@@ -155,7 +173,27 @@ $$
 f(\varepsilon,T)=\frac{1}{ {\rm e}^{(\varepsilon-\mu)/k_{\rm B}T}+1}
 $$
 
-![](figures/fermi_dirac.svg)
+```python
+fig = pyplot.figure()
+ax = fig.add_subplot(1,1,1)
+xvals = np.linspace(0, 2, 200)
+mu = .75
+beta = 20
+ax.plot(xvals, xvals < mu, ls='dashed', label='$T=0$')
+ax.plot(xvals, 1/(np.exp(beta * (xvals-mu)) + 1),
+        ls='solid', label='$T>0$')
+ax.set_xlabel(r'$\varepsilon$')
+ax.set_ylabel(r'$f(\varepsilon, T)$')
+ax.set_yticks([0, 1])
+ax.set_yticklabels(['$0$', '$1$'])
+ax.set_xticks([mu])
+ax.set_xticklabels([r'$\mu$'])
+ax.set_ylim(-.1, 1.1)
+ax.legend()
+draw_classic_axes(ax)
+pyplot.tight_layout()
+```
+
 
 Chemical potential $\mu=\varepsilon_{\rm F}$ if $T=0$. Typically $\varepsilon_{\rm F}/k_{\rm B}$~70 000 K (~7 eV), whereas room temperature is only 300 K (~30 meV) $\rightarrow$ thermal smearing occurs only very close to Fermi surface.
 
@@ -167,7 +205,31 @@ $$
 
 We can use this to calculate the electronic contribution to the heat capacity.
 
-![](figures/FD_DOS.svg)
+```python
+E = np.linspace(0, 2, 500)
+fig, ax = pyplot.subplots()
+ax.plot(E, np.sqrt(E), linestyle='dashed')
+ax.text(1.7, 1.4, r'$g(\varepsilon)\propto \sqrt{\varepsilon}$', ha='center')
+ax.fill_between(E, np.sqrt(E) * (E < 1), alpha=.3)
+
+n = np.sqrt(E) / (1 + np.exp(20*(E-1)))
+ax.plot(E, n)
+ax.fill_between(E, n, alpha=.5)
+w = .17
+ax.annotate(s='', xy=(1, 1), xytext=(1-w, 1),
+            arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+ax.text(1-w/2, 1.1, r'$\sim k_BT$', ha='center')
+ax.plot([1-w, 1+w], [1, 0], c='k', linestyle='dashed')
+ax.annotate(s='', xy=(1, 0), xytext=(1, 1),
+            arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+ax.text(1.2, .7, r'$g(\varepsilon_F)$', ha='center')
+ax.set_xticks([1])
+ax.set_xticklabels([r'$\varepsilon_F$'])
+
+ax.set_ylabel(r"$g(\varepsilon)$")
+ax.set_xlabel(r"$\varepsilon$")
+draw_classic_axes(ax, xlabeloffset=.2)
+```
 
 Electrons in the top triangle are being excited to the bottom triangle due to temperature increase. Number of excited electrons $\approx\frac{1}{2}g(\varepsilon_{\rm F})k_{\rm B}T=n_{\rm exc}$. Total extra energy $E(T)-E(0)=n_{\rm exc}k_{\rm B}T=\frac{1}{2}g(\varepsilon_{\rm F})k_{\rm B}^2T^2$.
 
diff --git a/docs/lecture_3.md b/src/lecture_3.md
similarity index 92%
rename from docs/lecture_3.md
rename to src/lecture_3.md
index f0c7645d9d4893c85f805fd4ffb14eaa29dc66f6..c8cbb70c9b95b1fdfa846e3c2865aa8b572df49a 100644
--- a/docs/lecture_3.md
+++ b/src/lecture_3.md
@@ -1,3 +1,14 @@
+```python
+from matplotlib import pyplot
+
+import numpy as np
+
+from common import draw_classic_axes, configure_plotting
+
+configure_plotting()
+
+pi = np.pi
+```
 
 # Atoms and bonds
 
@@ -246,7 +257,7 @@ Same as a single covalent bond, only more atoms in a line. Considering 3 atoms:
 $$
 E \begin{pmatrix}
 c_1 \\ c_2 \\ c_3
-\end{pmatrix} = 
+\end{pmatrix} =
 \begin{pmatrix}
 E_0 & t & 0 \\
 t & E_0 & t \\
@@ -263,22 +274,45 @@ Diagonalizing large matrices is unwieldy, but let's try and check it numerically
 
 Eigenfrequencies of 3 atoms: `[0.0 1.0 1.732050]`
 
-![](figures/3_phonons.svg)
+```python
+def DOS_finite_phonon_chain(n):
+    rhs = 2 * np.eye(n) - np.eye(n, k=1) - np.eye(n, k=-1)
+    rhs[0, 0] -= 1
+    rhs[-1, -1] -= 1
+    pyplot.figure()
+    pyplot.hist(np.sqrt(np.abs(np.linalg.eigvalsh(rhs))), bins=30)
+    pyplot.xlabel("$\omega$")
+    pyplot.ylabel("Number of levels")
+
+DOS_finite_phonon_chain(3)
+```
 
 Energies of 3 orbitals: `[-1.41421356 0.0  1.41421356]`
 
-![](figures/3_electrons.svg)
+```python
+def DOS_finite_electron_chain(n):
+    rhs = - np.eye(n, k=1) - np.eye(n, k=-1)
+    pyplot.figure()
+    pyplot.hist(np.linalg.eigvalsh(rhs), bins=30)
+    pyplot.xlabel("$E$")
+    pyplot.ylabel("Number of levels")
+
+DOS_finite_electron_chain(3)
+```
 
 ### From 3 atoms to 300
 
 Frequencies of many phonons:
 
-![](figures/300_phonons.svg)
+```python
+DOS_finite_phonon_chain(300)
+```
 
 Energies of many electrons:
 
-![](figures/300_electrons.svg)
-
+```python
+DOS_finite_electron_chain(300)
+```
 ## Summary
 
 * Electrons in atoms occupy "shells" in the energetically favorable order
diff --git a/docs/lecture_4.md b/src/lecture_4.md
similarity index 68%
rename from docs/lecture_4.md
rename to src/lecture_4.md
index a93448b2f18324c377e875f6b405a9f965978aa9..55be75d7e1d4bd68b6c7827765f151fa49b24517 100644
--- a/docs/lecture_4.md
+++ b/src/lecture_4.md
@@ -1,3 +1,16 @@
+```python
+import matplotlib
+from matplotlib import pyplot
+
+import numpy as np
+
+from common import draw_classic_axes, configure_plotting
+
+configure_plotting()
+
+pi = np.pi
+```
+
 # Electrons and phonons in 1D
 
 Last lecture:
@@ -10,7 +23,7 @@ This lecture:
 * Infinitely many atoms
 * Main idea: use *periodicity* in space, similar to periodicity in time
 
-## Equations of motion 
+## Equations of motion
 
 ### Phonons
 
@@ -53,7 +66,24 @@ In that case $c_n/c_0 = \exp(2 \pi n p a/L) = \exp(2 \pi n p/N)$, and changing $
 
 We can see that the solutions with different $k$ (or different $p$) are identical by plotting them:
 
-![](figures/phonons4.svg)
+```python
+x = np.linspace(-.2, 2.8, 500)
+fig, ax = pyplot.subplots()
+ax.plot(x, np.cos(pi*(x - .3)), label=r'$k=\pi/a$')
+ax.plot(x, np.cos(3*pi*(x - .3)), label=r'$k=3\pi/a$')
+ax.plot(x, np.cos(5*pi*(x - .3)), label=r'$k=5\pi/a$')
+sites = np.arange(3) + .3
+ax.scatter(sites, np.cos(pi*(sites - .3)), c='k', s=64, zorder=5)
+ax.set_xlabel('$x$')
+ax.set_ylabel('$u_n$')
+ax.set_xlim((-.1, 3.2))
+ax.set_ylim((-1.3, 1.3))
+ax.legend(loc='lower right')
+draw_classic_axes(ax)
+ax.annotate(s='', xy=(.3, -1.1), xytext=(1.3, -1.1),
+            arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+ax.text(.3 + .5, -1.25, '$a$', ha='center');
+```
 
 Let's count how many different solutions we expect to find. We have a system with $N$ degrees of freedom (either $u_n$ or $c_n$), and therefore we need $N$ normal modes (or eigenstates).
 
@@ -73,7 +103,27 @@ $$\omega = \sqrt{\frac{2\kappa}{m}}\sqrt{1-\cos(ka)} = 2\sqrt{\frac{\kappa}{m}}|
 
 So we arrive to the dispersion relation
 
-![](figures/phonons3.svg)
+```python
+k = np.linspace(-2*pi, 6*pi, 500)
+fig, ax = pyplot.subplots()
+
+pyplot.plot(k, np.abs(np.sin(k/2)))
+
+ax.set_ylim(bottom=0, top=1.2)
+ax.set_xlabel('$ka$')
+ax.set_ylabel(r'$\omega$')
+pyplot.xticks(list(2*pi*np.arange(-1, 4)) + [-pi, pi],
+              [r'$-2\pi$', '$0$', r'$2\pi$', r'$4\pi$', r'$6\pi$',
+               r'$-\pi$', r'$\pi$'])
+pyplot.yticks([1], [r'$2\sqrt{\frac{\kappa}{m}}$'])
+pyplot.vlines([-pi, pi], 0, 1.1, linestyles='dashed')
+pyplot.hlines([1], .1*pi, 1.3*pi, linestyles='dashed')
+draw_classic_axes(ax)
+ax.annotate(s='', xy=(-pi, -.15), xytext=(pi, -.15),
+            arrowprops=dict(arrowstyle='<->', shrinkA=0, shrinkB=0))
+ax.text(0, -.25, '1st Brillouin zone', ha='center')
+ax.set_ylim(bottom=-.7);
+```
 
 Here $k_p = 2\pi p / L$ for $0 ≤ p < N$ are *all the normal modes* of the crystal, and therefore we can rederive a better Debye model!
 
@@ -96,7 +146,14 @@ E = E_0 + 2t\cos(ka),
 $$
 so we come to the dispersion relation:
 
-![png](figures/tight_binding_1d.svg)
+```python
+pyplot.figure()
+k = np.linspace(-pi, pi, 300)
+pyplot.plot(k, -np.cos(k))
+pyplot.xlabel('$ka$'); pyplot.ylabel('$E$')
+pyplot.xticks([-pi, 0, pi], [r'$-\pi$', 0, r'$\pi$'])
+pyplot.yticks([-1, 0, 1], ['$E_0+2t$', '$E_0$', '$E_0-2t$']);
+```
 
 
 Usually electron dispersion has multiple options for $E(k)$, each called a *band*. The complete dispersion relation is also called a *band structure*.
@@ -147,8 +204,17 @@ Substituting, we get
 
 $$m_{eff} = \hbar^2\left(\frac{d^2 E(k)}{dk^2}\right)^{-1}$$
 
-![](figures/meff_1d_tb.svg)
-
+```python
+pyplot.figure(figsize=(8, 5))
+k = np.linspace(-pi, pi, 300)
+meff = 1/np.cos(k)
+color = list(matplotlib.rcParams['axes.prop_cycle'])[0]['color']
+pyplot.plot(k[meff > 0], meff[meff > 0], c=color)
+pyplot.plot(k[meff < 0], meff[meff < 0], c=color)
+pyplot.ylim(-5, 5)
+pyplot.xlabel('$ka$'); pyplot.ylabel('$m_{eff}$')
+pyplot.xticks([-pi, 0, pi], [r'$-\pi$', 0, r'$\pi$']);
+```
 
 ### density of states
 
@@ -251,7 +317,30 @@ $$\omega^2=\frac{\kappa(M+m)}{Mm}\pm \kappa\left\{\left(\frac{M+m}{Mm}\right)^2-
 
 Looking at the eigenvectors we see that for every $k$ there are now two values of $\omega$: one corresponding to in-phase motion (–) and anti-phase (+).
 
-![](figures/phonons6.svg)
+```python
+def dispersion_2m(k, kappa=1, M=1.4, m=1, acoustic=True):
+    Mm = M*m
+    m_harm = (M + m) / Mm
+    root = kappa * np.sqrt(m_harm**2 - 4*np.sin(k/2)**2 / Mm)
+    if acoustic:
+        root *= -1
+    return np.sqrt(kappa*m_harm + root)
+
+# TODO: Add panels with eigenvectors
+k = np.linspace(-2*pi, 6*pi, 300)
+fig, ax = pyplot.subplots()
+ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
+ax.plot(k, dispersion_2m(k), label='acoustic')
+ax.set_xlabel('$ka$')
+ax.set_ylabel(r'$\omega$')
+pyplot.xticks([-pi, 0, pi], [r'$-\pi$', '$0$', r'$\pi$'])
+pyplot.yticks([], [])
+pyplot.vlines([-pi, pi], 0, 2.2, linestyles='dashed')
+pyplot.legend()
+pyplot.xlim(-1.75*pi, 3.5*pi)
+pyplot.ylim(bottom=0)
+draw_classic_axes(ax, xlabeloffset=.2)
+```
 
 'A' is called the _acoustic branch_, 'B' is called the _optical branch_.
 
@@ -261,9 +350,46 @@ The quantity $\frac{ {\rm d}k}{ {\rm d}\omega}$ can be calculated from the dispe
 
 When plotted, the DOS may look something like this:
 
-![](figures/phonons8.svg)
-
-Note that normally $g(\omega)$ is plotted along the vertical axis and $\omega$ along the horizontal axis – the plot above is just to demonstrate the relation between the dispersion and the DOS. The singularities in $g(\omega)$ at the bottom and top of each branch are called _van Hove singularities_. 
+```python
+matplotlib.rcParams['font.size'] = 24
+k = np.linspace(-.25*pi, 1.5*pi, 300)
+k_dos = np.linspace(0, pi, 20)
+fig, (ax, ax2) = pyplot.subplots(ncols=2, sharey=True, figsize=(10, 5))
+ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
+ax.plot(k, dispersion_2m(k), label='acoustic')
+ax.vlines(k_dos, 0, dispersion_2m(k_dos, acoustic=False),
+          colors=(0.5, 0.5, 0.5, 0.5))
+ax.hlines(
+    np.hstack((dispersion_2m(k_dos, acoustic=False), dispersion_2m(k_dos))),
+    np.hstack((k_dos, k_dos)),
+    1.8*pi,
+    colors=(0.5, 0.5, 0.5, 0.5)
+)
+ax.set_xlabel('$ka$')
+ax.set_ylabel(r'$\omega$')
+ax.set_xticks([0, pi])
+ax.set_xticklabels(['$0$', r'$\pi$'])
+ax.set_yticks([], [])
+ax.set_xlim(-pi/4, 2*pi)
+ax.set_ylim((0, dispersion_2m(0, acoustic=False) + .2))
+draw_classic_axes(ax, xlabeloffset=.2)
+
+k = np.linspace(0, pi, 1000)
+omegas = np.hstack((
+    dispersion_2m(k, acoustic=False), dispersion_2m(k)
+))
+ax2.hist(omegas, orientation='horizontal', bins=75)
+ax2.set_xlabel(r'$g(\omega)$')
+ax2.set_ylabel(r'$\omega$')
+
+# Truncate the singularity in the DOS
+max_x = ax2.get_xlim()[1]
+ax2.set_xlim((0, max_x/2))
+draw_classic_axes(ax2, xlabeloffset=.1)
+matplotlib.rcParams['font.size'] = 16
+```
+
+Note that normally $g(\omega)$ is plotted along the vertical axis and $\omega$ along the horizontal axis – the plot above is just to demonstrate the relation between the dispersion and the DOS. The singularities in $g(\omega)$ at the bottom and top of each branch are called _van Hove singularities_.
 
 ### Consistency check with 1 atom per cell
 
@@ -275,7 +401,25 @@ For $M\rightarrow m$, the dispersion diagram should come back to the diagram obt
 
 To reconsile the two pictures let's plot two unit cells in the reciprocal space. (Definition: each of these unit cells is called a **Brillouin zone**, a concept that will come up later.)
 
-![](figures/phonons7.svg)
+```python
+k = np.linspace(0, 2*pi, 300)
+k_dos = np.linspace(0, pi, 20)
+fig, ax = pyplot.subplots()
+ax.plot(k, dispersion_2m(k, acoustic=False), label='optical')
+ax.plot(k, dispersion_2m(k), label='acoustic')
+omega_max = dispersion_2m(0, acoustic=False)
+ax.plot(k, omega_max * np.sin(k/4), label='equal masses')
+ax.set_xlabel('$ka$')
+ax.set_ylabel(r'$\omega$')
+ax.set_xticks([0, pi, 2*pi])
+ax.set_xticklabels(['$0$', r'$\pi/2a$', r'$\pi/a$'])
+ax.set_yticks([], [])
+ax.set_xlim(-pi/8, 2*pi+.4)
+ax.set_ylim((0, dispersion_2m(0, acoustic=False) + .2))
+ax.legend(loc='lower right')
+pyplot.vlines([pi, 2*pi], 0, 2.2, linestyles='dashed')
+draw_classic_axes(ax, xlabeloffset=.2)
+```
 
 Doubling the lattice constant "folds" the band structure on itself, doubling all the bands.
 
@@ -289,7 +433,7 @@ $\rightarrow \rho_{\rm S}(k)=$ (no. of $k$-values / unit of $k$) $=\frac{L}{\pi}
 
 For open boundary conditions $\rightarrow$ running waves $\rightarrow k=$ $...\frac{-4\pi}{L},\frac{-2\pi}{L},0,\frac{2\pi}{L},\frac{4\pi}{L}...$
 
-$\rightarrow \rho_{\rm R}(k)=\frac{L}{2\pi}$, which is lower than for the case of closed boundary conditions, however, in this case the entire $k$-space is filled. 
+$\rightarrow \rho_{\rm R}(k)=\frac{L}{2\pi}$, which is lower than for the case of closed boundary conditions, however, in this case the entire $k$-space is filled.
 
 ## Summary
 
diff --git a/docs/lecture_5.md b/src/lecture_5.md
similarity index 100%
rename from docs/lecture_5.md
rename to src/lecture_5.md
diff --git a/docs/lecture_6.md b/src/lecture_6.md
similarity index 100%
rename from docs/lecture_6.md
rename to src/lecture_6.md
diff --git a/docs/lecture_7.md b/src/lecture_7.md
similarity index 91%
rename from docs/lecture_7.md
rename to src/lecture_7.md
index f4ab67aadb72d74e1b5e13cbabcea121e0b38d3f..d7949cfdc11880e8bbbfaa08e92ed40633d6faa2 100644
--- a/docs/lecture_7.md
+++ b/src/lecture_7.md
@@ -1,3 +1,23 @@
+```python
+from matplotlib import pyplot
+
+import numpy as np
+from scipy.optimize import curve_fit
+from scipy.integrate import quad
+
+from common import draw_classic_axes, configure_plotting
+
+configure_plotting()
+
+def sqrt_plus(x):
+    return np.sqrt(x * (x >= 0))
+
+# Band structure parameters.
+E_V, E_C, E_F = -1.2, 1.8, .4
+E_D, E_A = E_C - .7, E_V + .5
+m_h, m_e = 1, .5
+```
+
 # Semiconductor physics
 
 ## Review of band structure properties
@@ -102,7 +122,27 @@ Observe that because we are describing particles in the valence band as holes, $
 
 ### Part 1: pristine semiconductor
 
-![](figures/intrinsic_DOS.svg)
+```python
+E = np.linspace(-3, 3, 1000)
+fig, ax = pyplot.subplots()
+
+n_F = 1/(np.exp(2*(E - E_F)) + 1)
+g_e = m_e * sqrt_plus(E - E_C)
+g_h = m_h * sqrt_plus(E_V - E)
+ax.plot(E, g_h, label="$g_e$")
+ax.plot(E, g_e, label="$g_h$")
+ax.plot(E, 10 * g_h * (1-n_F), label="$n_h$")
+ax.plot(E, 10 * g_e * n_F, label="$n_e$")
+ax.plot(E, n_F, label="$n_F$", linestyle='dashed')
+ax.set_ylim(top=1.5)
+
+ax.set_xlabel('$E$')
+ax.set_ylabel('$g$')
+ax.set_xticks([E_V, E_C, E_F])
+ax.set_xticklabels(['$E_V$', '$E_C$', '$E_F$'])
+ax.legend()
+draw_classic_axes(ax, xlabeloffset=.2)
+```
 
 Electrons
 $$ E = E_G + {p^2}/{2m_e}$$
@@ -198,7 +238,26 @@ Likewise an acceptor adds an extra state at $E_A$ (close to the top of the valen
 
 ### Density of states with donors and acceptors
 
-![](figures/doping.svg)
+```python
+E = np.linspace(-3, 3, 1000)
+fig, ax = pyplot.subplots()
+
+n_F = 1/(np.exp(2*(E - E_F)) + 1)
+g_e = m_e * sqrt_plus(E - E_C)
+g_h = m_h * sqrt_plus(E_V - E)
+ax.plot(E, g_h, label="$g_e$")
+ax.plot(E, g_e, label="$g_h$")
+
+sigma = 0.01
+g_D = np.exp(-(E_D - E)**2 / sigma**2)
+g_A = .7 * np.exp(-(E_A - E)**2 / sigma**2)
+ax.plot(E, g_D, label='$g_D$')
+ax.plot(E, g_A, label='$g_A$')
+ax.legend()
+ax.set_xticks([E_V, E_C, E_A, E_D])
+ax.set_xticklabels(['$E_V$', '$E_C$', '$E_A$', '$E_D$'])
+draw_classic_axes(ax, xlabeloffset=.2)
+```
 
 All donor/acceptor states at the same energy:
 $$g_A(E) = N_A \delta(E-E_A),\quad g_D(E) = N_D \delta(E- (E_G - E_D))$$
@@ -336,7 +395,9 @@ See the book for details.
 
 Density of states in a doped semiconductor:
 
-![](figures/doping.svg)
+```python
+fig
+```
 
 Charge balance determins the number of electrons and holes as well as the position of the Fermi level.
 
diff --git a/docs/lecture_8.md b/src/lecture_8.md
similarity index 95%
rename from docs/lecture_8.md
rename to src/lecture_8.md
index 894df24b84e0620b222d5f4a28cc10600eb11864..9fda813c226665c4874bb4de9a1283ffb1df3823 100644
--- a/docs/lecture_8.md
+++ b/src/lecture_8.md
@@ -1,3 +1,15 @@
+```python
+from matplotlib import pyplot
+
+import numpy as np
+from scipy.optimize import curve_fit
+from scipy.integrate import quad
+
+from common import draw_classic_axes, configure_plotting
+
+configure_plotting()
+```
+
 # Lecture 8 – Magnetism
 
 _Based on chapters 19 and 20 of the book_
@@ -81,7 +93,18 @@ $$
 
 with $C$ a constant. This is known as the _Curie law_, and it shows that at high temperature a paramagnetic material becomes less susceptible to magnetisation.
 
-![](figures/curie.svg)
+```python
+B = np.linspace(-2, 2, 1000)
+fig, ax = pyplot.subplots()
+
+ax.plot(B, np.tanh(B * 3), label="low $T$")
+ax.plot(B, np.tanh(B), label="high $T$")
+
+ax.legend()
+ax.set_ylabel("$M$")
+ax.set_xlabel("$B$")
+draw_classic_axes(ax, xlabeloffset=.2)
+```
 
 ### Paramagnetism: free spin $J$
 For an atom with magnetic moments $L$ and $S$ the Hamiltonian will become:
@@ -193,4 +216,4 @@ $$
 {\mathcal H}=\sum_i \left( -J \sigma_i \sigma_{i+1}+ g\mu_{\rm B}B\sigma_i\right)
 $$
 
-where $\sigma_i=\pm S$.
\ No newline at end of file
+where $\sigma_i=\pm S$.