diff --git a/src/3_drude_model.md b/src/3_drude_model.md
index b797c2685be92a9729a585c33807d160474e9bde..02cd6a08163e126f6d800b9a884e0afbc7225cd2 100644
--- a/src/3_drude_model.md
+++ b/src/3_drude_model.md
@@ -25,6 +25,45 @@ Ohm's law states that $V=IR=I\rho\frac{l}{A}$. In this lecture we will investiga
 - At each scattering event an electron returns to momentum ${\bf p}=0$.
 - In-between scattering events electrons respond to the Lorentz force ${\bf F}_{\rm L}=-e\left({\bf E}+{\bf v}\times{\bf B}\right)$.
 
+```python
+import numpy as np
+import matplotlib.pyplot as plt
+
+walker_number = 20 # number of particles
+tau = 1 # relaxation time
+gamma = .3 # dissipation strength
+a = 1 # acceleration
+dt = .1 # infinitesimal
+T = 20 # simulation time 
+
+v = np.zeros((2, int(T // dt), walker_number), dtype=float) 
+
+scattering_events = np.random.binomial(1, dt/tau, size=v.shape[1:])
+angles = np.random.uniform(high=2*np.pi, size=scattering_events.shape) * scattering_events
+rotations = np.array(
+    [[np.cos(angles), np.sin(angles)],
+     [-np.sin(angles), np.cos(angles)]]
+)
+
+for step in range(1, v.shape[1]):
+    v[:, step] = v[:, step-1]
+    v[0, step] += a * dt
+    v[:, step] = np.einsum(
+        'ijk,jk->ik',
+        rotations[:, :, step-1, :],
+        v[:, step, :]
+    ) * (1 - gamma * scattering_events[step-1])
+
+r = np.cumsum(v * dt, axis=1)
+
+scattering_positions = np.copy(r)
+scattering_positions[:, ~scattering_events.astype(bool)] = np.nan
+
+plt.plot(*r[:, :100], alpha=.5, c='#1f77b4');
+plt.scatter(*scattering_positions[:, :100], s=10);
+plt.axis('off');
+```
+
 We start by considering only an electric field (_i.e._ ${\bf B}=0$). What velocity do electrons acquire in-between collisions?
 
 $$
@@ -124,4 +163,4 @@ $$\mathbf{E} = \begin{pmatrix} \rho_{xx} & \rho_{xy} \\ \rho_{yx} & \rho_{yy} \e
   2. Invert the resistivity matrix to obtain the conductivity matrix $$\begin{pmatrix} \sigma_{xx} & \sigma_{xy} \\ \sigma_{yx} & \sigma_{yy} \end{pmatrix} $$, allowing you to express $\mathbf{J}$ as a function of $\mathbf{E}$. 
   3. Sketch $\sigma_{xx}$ and $\sigma_{xy}$ as a function of the magnetic field $\bf B$. 
   4. Give the definition of the Hall coefficient. What does the sign of the Hall coefficient indicate?
-  
\ No newline at end of file
+