A large number of studies have been published on the topic of acoustic wavefield modeling in anisotropic media. All of them are based on the choice of the suitable wave equation for numerical implementation. However, these wave equations are usually cumbersome, have an unclear physical nature, are computationally demanding, and generate artificial pseudo shear modes, which are considered as artifacts in the seismic imaging process. Duveneck and Bakker (2011) derived a system of coupled differential wave equations based on Hooke's law and equation of motion only. Despite all the advantages, these equations are unstable for a certain configuration of anisotropic parameters and generate S-wave artifacts. Liu et al. (2009), on the other hand, derived an unconditionally stable single wave equation that turned out to be difficult to model. Moreover, it is responsible only for the P-wave mode. Nikonenko and Charara (2020) have shown that this single wave equation is just one mode for the Duveneсk coupled equations and proposed a possible fully explicit scheme for its solution. We continue this approach making the solution optimal and extending it to other cases of anisotropy. Numerical examples illustrate the absence of artifacts and the accuracy of the proposed method.