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. (C) 2017 Elsevier Ltd. All rights reserved.