We present a new numerical method for elastic wave modeling in 3D isotropic and anisotropic media, which is called the 3D optimal nearly analytic discrete method (onadm) in this article. This work is an extension of the 2D onadm (Yang et al., 2004) that models acoustic and elastic waves propagating in 2D media. The formulation is derived by using a multivariable truncated Taylor series expansion and high-order interpolation approximations. Our 3D onadm enables wave propagation to be simulated in three dimensions through isotropic and anisotropic models. Promising numerical results show that the error of the onadm for the 3D case is less than those of the conventional finite-difference (fd) method and the so-called Lax–Wendroff correction (lwc) schemes, measured quantitatively by the root- mean-square deviation from analytical solution. The seismic wave fields in the 3D isotropic and anisotropic media are simulated and compared with those obtained by using the fourth-order lwc method and exact solutions for the acoustic wave case. We show that, compared with the conventional fd method and the lwc schemes, the onadm for the 3D case can reduce significantly the storage space and computational costs. Numerical experiments illustrate that the onadm provides a useful tool for the 3D large-scale isotropic and anisotropic problems and it can suppress effectively numerical dispersions caused by discretizing the 3D wave equations when too-coarse grids are used, which is the same as the 2D onadm. Numerical modeling also implies that simultaneously using both the wave displacement and its gradients to approximate the high-order derivatives is important for both decreasing the numerical dispersion caused by the discretization of wave equations and compensating the important wave-field information included in wave-displacement gradients.