We have presented a recursive integral time extrapolation method for modeling elastic wave propagation in general heterogeneous anisotropic media. The one-step formulation involves analytical wavefields that contain either positive or negative frequencies, which provides crucial information about the direction of wave propagation. The two-step formulation, on the other hand, involves only real-valued wavefields. The proposed method employs a low-rank approximation to efficiently apply the Fourier integral operators defined by the mixed-domain components of the modified Christoffel matrix. In practice, this reduces the computational cost to a small number of spatial fast Fourier transforms per time step. The low-rank decomposition only needs to be computed once prior to wave extrapolation. Numerical examples show that the proposed method has superior accuracy and stability compared with conventional finite-difference and pseudospectral methods.