You cannot handle the dispersive term (third order derivative) explicitly as you do, it leads definitely to instability.
To make it implicit you can linearize rho around the value at the old time step. This gives
vdw1.edp (1.8 KB)
In particular you solve in a hidden way a bilaplace type equation in rho.
Your problem is stiff because of the dispersive term, and overall because of the temperature T that you take less than 1. T less than 1 implies that the underlying pressure/velocity wave system is not hyperbolic, hence it is very much unclear if the whole system is well-posed. Then numerical results are not guaranteed at all!
For T greater than 1 I am more confident.
Anyway you get numerical oscillations because of the wave system (to avoid them you would need upwinding, which is not easily done in FreeFem).