Thermal phase change problems are widespread in mathematics, nature, and science. They are particularly useful in simulating the phenomena of melting and solidification in materials science. In this paper we propose a novel class of arbitrarily high-order and unconditionally energy stable schemes for a thermal phase change model, which is the coupling of a heat transfer equation and a phase field equation. The unconditional energy stability and consistency error estimates are rigorously proved for the proposed schemes. A detailed implementation demonstrates that the proposed method requires only the solution of a system of linear elliptic equations at each time step, with an efficient scheme of sufficient accuracy to calculate the solution at the first step. It is observed from the comparison with the classical explicit Runge-Kutta method that the new schemes allow to use larger time steps. Adaptive time step size strategies can be applied to further benefit from this unconditional stability. Numerical experiments are presented to verify the theoretical claims and to illustrate the accuracy and effectiveness of our method.