diff --git a/PS10_Jung/PS10.py b/PS10_Jung/PS10.py
index 25c76158d23905f1cd46dde4a24a29151adcb845..2e8dbb3ab3e00206131055f4c0b704d04f20ccc0 100644
--- a/PS10_Jung/PS10.py
+++ b/PS10_Jung/PS10.py
@@ -7,6 +7,8 @@ plt.rcParams.update({
     "font.sans-serif": "Helvetica",
 })
 
+from tqdm import tqdm
+
 # a - Implementation based on Leimkuhler & Matthews
 gamma = 1
 omega_0 = 1/3 * gamma
@@ -15,8 +17,8 @@ d = 3
 k_BT = 1
 l = 1 # gamma**-1 * np.sqrt(k_BT / m)
 N = 100
-t_max = 1000 
-dt = 0.1 
+t_max = 1000
+dt = 0.1
 t = np.arange(0,t_max,dt)
 
 def force(r: np.ndarray) -> np.ndarray:
@@ -35,7 +37,7 @@ def BAOAB(force, x0, p0, m, dt, N):
         r = np.random.normal(0,1,(N,d))
         p[i] = p[i-1] + 1/2 * force(x[i-1]) * dt # same
         x[i] = x[i-1] + p[i] / (2*m) * dt # diff
-        p[i] = np.exp(-gamma*dt) * p[i] + xi_2 * r/np.sqrt(m)
+        p[i] = np.exp(-gamma*dt) * p[i] + xi_2 * r * np.sqrt(m)
         x[i] = x[i] + (dt/2) * p[i]/m
         p[i] = p[i] + 1/2 * force(x[i]) * dt
 
@@ -94,4 +96,23 @@ E_pot = (m/(2*N)) * omega_0**2 * np.linalg.norm(r, axis=(1,2))**2
 #     plt.plot(hist_p, label=f'p_{coords[i]}')
 #     plt.legend()
 #     plt.tight_layout()
-#     plt.savefig(f'p_{coords[i]}_dist.png')
\ No newline at end of file
+#     plt.savefig(f'p_{coords[i]}_dist.png')
+
+# c)
+r_f, p_f = BAOAB(lambda r: 0, r0, p0, m, dt, N)
+
+dr2 = np.zeros(t.shape)
+dr2_f = np.zeros(t.shape)
+for i in tqdm(range(1, len(dr2))):
+    dr2[i] = np.average(np.linalg.norm(r[i:] - r[:-i], axis=2)**2)
+    dr2_f[i] = np.average(np.linalg.norm(r_f[i:] - r_f[:-i], axis=2)**2)
+
+plt.loglog(t, dr2, label="harmonic")
+plt.loglog(t, dr2_f, label="free")
+plt.loglog(t, 3 * k_BT / m * t**2, ls="--", label=r"$3 k_B T t^2 / m$")
+plt.loglog(t, 6 * k_BT * t / (m * gamma) - 6 * k_BT / (m * gamma**2) * (1 - np.exp(-gamma * t)), ls="-.", label="free, analytic")
+plt.legend()
+plt.xlabel(r"$\log t$")
+plt.ylabel(r"$\log \delta r^2$")
+plt.savefig("problem10.1c.png", dpi=300)
+plt.show()