In this paper, we propose a high‐order stereo‐modeling method (STEM) to approximate the high‐order spatial derivatives included in the wave equations using simultaneously wave‐field displacements and their gradients and propose a two‐step method for time marching, which is called the two‐step STEM in brief. The two‐step STEM has a higher‐order accuracy in space than conventional finite‐difference (FD) methods when the same number of spatial grid points in a wavelength is used. For example, the stereo‐modeling method that uses five points in one spatial direction can achieve an eighth‐order accuracy in space, whereas other FD methods such as conventional FD methods, Lax–Wendroff correction (LWC) methods, and other methods only have a fourth‐order accuracy. Theoretical properties of the two‐step STEM including stability and errors are analyzed for 1D and 2D cases. The numerical dispersion relationship provided by the two‐step STEM for 1D and 2D cases are also investigated in this study. Meanwhile, we present numerical results computed by the two‐step STEM and compare with the eighth‐order LWC method, the eighth‐order staggered‐grid FD method, and the fourth‐order staggered‐grid method. Numerical results show that the two‐step STEM can effectively suppress numerical dispersion caused by discretizing the wave equations when coarse spatial grids are used or models have strong velocity contrasts between adjacent layers. In contrast to other high‐order FD methods such as the eighth‐order LWC, the eighth‐order staggered‐grid FD, and the fourth‐order staggered‐grid method, the new method has substantially less computational time and requires less memory because large spatial and time increments can be used. Thus, the two‐step STEM can be potentially used to solve large‐scale wave‐propagation problems and seismic tomography.