A computationally efficient numerical scheme is presented for the phase-field model of two-phase systems for anisotropic interfacial energy. The scheme is solved by using a nonlinear multigrid method. When the coefficient for the anisotropic interfacial energy is sufficiently high, the interface of the system shows corners or missing crystallographic orientations. Numerical simulations with high and low anisotropic coefficients show excellent agreement with exact equilibrium shapes. We also present spinodal decomposition, which shows the robustness of the pro-posed scheme.