The Cahn-Hilliard phase field method for two-phase flow has gained particular attention due to its unique features including its flexibility for complex morphological and topological changes, the intrinsic property of conserving mass, and the natural approach to account for the surface tension. The essential idea of the method is to use a phase field function to describe the two-phase system, while a thin smooth transition layer (interfacial area) connects the two immiscible fluids, where the value of phase field function varies continuously. The application of phase field method to two-phase flows has become more widespread recently, but to the best of our knowledge, very little progress has been made for the method being applied to the two-phase flows with phase change. This includes evaporation, condensation and boiling, which plays an important role in enhanced heat transfer in power electronics, energy, and aerospace engineering. In previous work, in order to face the challenge of large density contrast and high Reynolds number of practical engineering problems, we developed a stabilized phase field method that can handle two-phase flow with density ratio over 1000, at high Reynolds number over 10,000, and applied the method to simulate slug initiation in a long circular pipe. In this paper, inspired by the boiling model widely used in the level-set method, we propose a new boiling model that assumes that boiling takes place in the whole interfacial layer. The method is then used to solve the non-solenoidal Navier-Stokes equations. The boiling model is validated by simulating a vapor bubble growing in super-heated liquid. For this case, the growth rate of the bubble has an analytical solution, and it is used as a benchmark case in volume of fluid (VOF) and level-set methods extensively. For three different refrigerants, namely water, R134a and HFE7100, our phase field method with the boiling model can obtain accurate simulation results. Moreover, the method and model are applied to predict the three-dimensional boiling heat transfer in a rectangular micro-channel that contains a water vapor bubble with various inlet super-heat conditions. We found that the predicted bubble shape is very similar to that visualized in existing experiment. From our simulation of boiling flow using the phase field method, We have found that the required mesh resolution for the phase field method is comparable with that of VOF and level-set methods.