In this work, we consider a convexity splitting scheme for a coupled phase field and energy equation, a modification of Stefan problem. The Stefan problem is a free boundary value problem that models the temperature in a homogeneous multiphase medium. Each phase is modeled using a heat diffusion parabolic partial differential equation and incorporates latent heat at the phase interface. The existence of the jump in the heat flux due to the latent heat leads to a nonlinear free boundary problem. The Stefan problem has a sharp interface, thus coupling the Stefan problem with a phase field model allows for a more realistic diffuse region of phase transition. Solving this phase field and temperature system implicitly can by computationally expensive and using time lagging methods can be unstable, thus we look into the viability of using an efficient and stable convexity splitting scheme for the coupled system, in several variants.