This work presents a new scheme for wave propagation simulation in three-dimensional (3-D) elastic-anisotropic media. The algorithm is based on the rapid expansion method (REM) as a time integration algorithm, and the Fourier pseudospectral method for computation of the spatial derivatives. The REM expands the evolution operator of the second-order wave equation in terms of Chebychev polynomials, constituting an optimal series expansion with exponential convergence. The modeling allows arbitrary elastic coefficients and density in lateral and vertical directions.Numerical methods which are based on finite-difference techniques (in time and space) are not efficient when applied to realistic 3-D models, since they require considerable computer memory and time to obtain accurate results. On the other hand, the Fourier method permits a significant reduction of the working space, and the REM algorithm gives machine accuracy with half the computational effort as the usual second-order temporal differencing scheme. The new algorithm provides spectral accuracy for band limited wave propagation with no temporal or spatial dispersion. Hence, the combination REM/Fourier method could be considered at present the fastest and the most accurate of all the algorithms based on spectral methods in terms of efficiency of computer time. The technique is successfully tested with the analytical solution in the symmetry axis of a 3-D homogeneous transversely isotropic solid. Computed radiation patterns in clay shale and sandstone show the characteristics predicted by the theory. A final example computes synthetic seismograms showing the effects of shear-wave splitting of a model where an azimuthally anisotropic cracked shale layer is inside a transversely isotropic sandstone.