We present a procedure for simulating epitaxial growth based on the phase-field method. We consider a basic model in which growth is initiated by a flux of atoms onto a heated surface. The deposited atoms diffuse in the presence of this flux and eventually collide to form islands which grow and decay by the attachment and detachment of migrating atoms at their edges. Our implementation of the phase-field method for this model includes uniform deposition, isotropic surface diffusion, and stochastic nucleation (in both space and time), which creates islands whose boundaries evolve as the surface atoms "condense" into and "evaporate" from the islands. Computations using this model in the submonolayer regime, prior to any appreciable coalescence of islands, agree with the results of kinetic Monte Carlo (KMC) simulations for the coverage-dependence of adatom and island densities and island-size distributions, for both reversible and irreversible growth. The scaling of the island density, as obtained from homogeneous rate equations, agrees with KMC simulations for irreversible growth and for reversible growth for varying deposition flux at constant temperature. For reversible growth with varying temperature but constant flux, agreement relies on an estimate of the formation energy of the critical cluster. Taken together, our results provide a comprehensive analysis of the phase-field method in the submonolayer regime of epitaxial growth, including the verification of the main scaling laws for adatoms and island densities and the scaling functions for island-size distributions, and point to the areas where the method can be extended and improved.