diff --git a/content/kwant_workflow.py b/content/kwant_workflow.py
new file mode 100644
index 0000000000000000000000000000000000000000..62df7be38864057b2b17db9fd1709fedc67ecb0f
--- /dev/null
+++ b/content/kwant_workflow.py
@@ -0,0 +1,42 @@
+import kwant
+import matplotlib.pyplot as pyplot
+from scipy import misc
+import numpy as np
+
+# shape.png is an image containing the system shape colored in black.
+im = misc.imread('shape.png')[:, :, 0]
+w, h = im.T.shape
+
+sys = kwant.Builder()
+lat = kwant.lattice.square()
+
+for pos in np.argwhere(im != 255):
+    sys[lat(pos[1], -pos[0])] = 4
+
+def make_wire_shape(period, center, r, ymin, ymax):
+    direction = np.array(period / np.linalg.norm(period))
+    def wire_shape(pos):
+        rel_pos = pos - np.array(center)
+        projection = direction * np.dot(direction, rel_pos) - rel_pos
+        return sum(projection**2) <= r**2 and ymin < pos[1] < ymax
+
+    return wire_shape
+
+wire_dir = [(-1., -1.), (3., -4.), (2., 4.)]
+wire_center = [(h/2 - 50, -w/2 + 50), (h/2, -w/2 + 180),
+               (h/2 + 30, -w/2 + 200)]
+ylims = [(-400, -200), (-440, -100), (-100, 60)]
+
+for direction, center, ylim in zip(wire_dir, wire_center, ylims):
+    sys[lat.shape(make_wire_shape(direction, center, 50, *ylim), center)] = 4
+    lead = kwant.Builder(kwant.TranslationalSymmetry(direction))
+    lead[lat.wire(center, 50)] = 4
+    lead[lat.neighbors()] = -1
+    sys.attach_lead(lead)
+
+sys[lat.neighbors()] = -1
+
+sys = sys.finalized()
+
+wfs = kwant.wave_function(sys, energy=0.14)(0)[0]
+kwant.plotter.map(sys, abs(wfs), oversampling=1, cmap="PuBu")
diff --git a/content/shape.png b/content/shape.png
new file mode 100644
index 0000000000000000000000000000000000000000..250612f5607f96ef197e4ca93e74ad3f71b27daf
Binary files /dev/null and b/content/shape.png differ