diff --git a/src/1_einstein_model.md b/src/1_einstein_model.md
index ae615fb5d2fe65acfa8fb9e01e372011246e8434..0e825e63eb8726bd38f1ab546bbacb8fb550264e 100644
--- a/src/1_einstein_model.md
+++ b/src/1_einstein_model.md
@@ -108,8 +108,6 @@ V = 0.5*(omega**2)*(x**2)
 
 fig, ax = pyplot.subplots(figsize=(10, 7))
 
-ax.plot(x, V) #plot harmonic potential
-
 for i in range(no_states):
     
     ax.hlines(h0_ener(i), x[0], x[len(x)-1], linestyles='dotted', colors='k')
@@ -139,8 +137,10 @@ for i in range(no_states):
                     xytext=(x[0]+1/2, h0_ener(i-1)), 
                     arrowprops=dict(arrowstyle="<->"))
         
-    ax.fill_between(x, h0_ener(i), ho_evec(x, i, no_states) + h0_ener(i))
+    ax.fill_between(x, h0_ener(i), ho_evec(x, i, no_states) + h0_ener(i), alpha=0.5)
         
+ax.plot(x, V, 'k', linewidth=1) #plot harmonic potential
+
 # Move left y-axis and bottim x-axis to centre, passing through (0,0)
 ax.spines['left'].set_position('center')
 ax.spines['bottom'].set_position(('data', 0.0))