A phase-field model is used to simulate the dynamics of step-flow during epitaxial growth. An efficient numerical method is developed using a second order accurate fully implicit time discretization together with a second order accurate finite difference spatial discretization. To test the algorithm, we focus on simulations of kinetic instabilities that arise due to the anisotropy of adsorption (attachment) rates at steps from the upper and lower terraces during growth. When the attachment rate from the lower terrace is larger than that from the upper terrace (Ehrlich-Schwoebel effect), step meandering occurs. When the opposite is true of the rates, step bunching occurs. Step-step interactions are seen to reduce the meandering instability.