Fractional differential equations are becoming increasingly popular as a modelling tool to describe a wide range of non-classical phenomena with spatial heterogeneities throughout the applied sciences and engineering. However, the non-local nature of the fractional operators causes essential difficulties and challenges for numerical approximations. We here investigate the numerical solution of fractional-in-space phase-field models such as Allen-Cahn and Cahn-Hilliard equations via the contour integral method (CIM) for computing the fractional power of a matrix times a vector. Time discretization is performed by the first-and second-order implicit-explicit schemes with an adaptive time-step size approach, whereas spatial discretization is performed by a symmetric interior penalty Galerkin (SIPG) method. Several numerical examples are presented to illustrate the effect of the fractional power.