In this article, we propose a 3D split-step algorithm (SSA) for seismic-wave simulation. We first transform the wave equations in 3D anisotropic media into a system of first-order partial differential equations with respect to time t. Then we use the multidimensional high-order interpolation method to approximate the high-order spatial derivatives, so that we obtain a system of semidiscrete ordinary differential equations (ODEs). Finally, the third-order implicit Adams method and truncated differentiator series method are applied to solve the semidiscrete ODEs. We provide the theoretical study on the properties of the 3DSSA, such as stability criteria, theoretical error, numerical error, numerical dispersion, and computational efficiency. We also compare some seismic modeling results of this method against those of some high-order finite-difference schemes. Theoretical analysis and numerical tests show that the 3DSSA is third-order accurate in time and fourth-order accurate in space. However, its computational costs and memory requirements are much less than those of the fourth-order Lax–Wendroff correction method and the fourth-order staggered-grid method. Using a multilayer elastic model with large velocity contrasts and free surface, we compare the result of the 3DSSA with that of the discrete-wavenumber method. We also present the synthetic seismograms in the 3D three-layer isotropic medium, the wave-field snapshots in the 3D two-layer medium, and the 3D transversely isotropic medium with a vertical symmetry axis. All these promising numerical results illustrate that the 3DSSA can suppress effectively the numerical dispersion caused by discretizing the wave equations when too few sampling points per minimum wavelength are used or when models have large velocity contrasts between adjacent layers, further resulting in both increasing the computational efficiency and saving the storage space when big grids are used.