Equivalent staggered-grid (ESG) as a new family of schemes has been utilized in seismic modeling, imaging, and inversion. Traditionally, the Taylor series expansion is often applied to calculate finite-difference (FD) coefficients on spatial derivatives, but the simulation results suffer serious numerical dispersion on a large frequency zone. We develop an optimized equivalent staggered-grid (OESG) FD method that can simultaneously suppress temporal and spatial dispersion for solving the second-order system of the 3D elastic wave equation. On the one hand, we consider the coupling relations between wave speeds and spatial derivatives in the elastic wave equation and give three sets of FD coefficients with respect to the P-wave, S-wave, and converted-wave (C-wave) terms. On the other hand, a novel plane wave solution for the 3D elastic wave equation is derived from the matrix decomposition method to construct the time–space dispersion relations. FD coefficients of the OESG method can be acquired by solving the new dispersion equations based on the Newton iteration method. Finally, we construct a new objective function to analyze P-wave, S-wave, and C-wave dispersion concerning frequencies. The dispersion analyses show that the presented method produces less modeling errors than the traditional ESG method. The synthetic examples demonstrate the effectiveness and superiority of the presented method.