In heteroepitaxy, the mismatch of lattice constants between the crystal film and the substrate causes misfit strain and stress in the bulk of the film, driving the surface of the film to self-organize into various nanostructures. Below the roughening transition temperature, an epitaxial surface consists of facets and steps and changes its morphology by lateral motion of steps. In this paper, we present a $2+1$-dimensional continuum model for the long-range elastic interaction on stepped surface of a strained film. The continuum model is derived rigorously from the discrete model for the interaction between steps; thus it incorporates the discrete features of the stepped surfaces. Examples show that our continuum model is much more accurate as an approximation to the discrete model than the traditional continuum approximation. Moreover, in the linear instability of a planar surface, our continuum model gives the transition from step bunching instability to step undulation instability as the distance between adjacent steps increases, which agrees with the experimental observations and the results of discrete models and is missing using the traditional continuum approximation. Numerical simulations of the surface evolution using our model in the nonlinear regime show several different surface morphologies, including the morphology of step bunching which cannot be obtained using the traditional continuum approximation.