Abstract Phase field methods have become of great interest for the simulation of droplet and bubble dynamics, moving free interfaces and more recently phase change phenomena. One example is the Navier–Stokes–Korteweg (NSK) equations. With a particular focus on the van der Waals fluid, we present the numerical solution of the NSK system with the energy equation included. The least-squares spectral element formulation with a time-stepping procedure, a high-order continuity approximation and an element-by-element technique is implemented to provide a general and robust solver for the thermal NSK equations. A convergence analysis is conducted to verify our solver. Two numerical examples regarding phase transitions of a droplet and thermocapillary convection are provided to validate our solver.