Phase field equations are used to model a wide range of multiphase problems such as separation of fluids, solidification, viscous fingering, fracture and fatigue. A wide variety of methods to numerically solve phase field equations can be found in the literature. In particular, high order methods are an effective option when accuracy improvement is desired. In the first part of this work, we analyze the accuracy and computational efficiency of the high order finite element method (FEM) and discontinuous Galerkin (DG) method applied to the second-order Allen–Cahn (AC) and fourth-order Cahn–Hilliard (CH) equations. Several schemes for time integration are used for these equations. The explicit schemes are the forward Euler, classical fourth-order Runge–Kutta (RK4) and the strong stability preserving ten stages fourth-order Runge–Kutta (RKSSP-10,4) described in Gottlieb et al. (2011). The backward Euler and trapezoidal implicit methods are adopted in the full and semi implicit schemes, as proposed in Eyre (unpublished). Manufactured solutions for one dimensional problems are used in order to evaluate the errors and to compare the different numerical methods. By choosing an adequate discretization for AC equations resulting from the previous analysis of the first part of the work, in the second part, we propose a numerical semi implicit scheme to solve the damage and fracture model described in Boldrini et al. (2016). This procedure employs the FEM for spatial discretization, the Newmark method for time integration of the kinematics equation and the backward Euler for the damage phase field evolution. Finally, results for 2D benchmark tests are presented for the fracture phase field model and the convergence to a sharp crack for a small width of the damage phase field layer γ is verified.