We present efficient and accurate modelling of seismic wavefields in anelastic media. We use a first‐order viscoelastic wave equation based on the generalized Robertsson's formulation. In our work, viscoacoustic and viscoelastic wave equations use the standard linear solid mechanism. To numerically solve the first‐order wave equation, we employed a scheme derived from the Lie product formula, where the time evolution operator of the analytic solution is written as a product of exponential matrices, and each exponential matrix term is approximated by the Taylor series expansion. The accuracy of the proposed scheme is evaluated by comparison with the analytical solution for a homogeneous medium. We also present simulations of some geological models with different structural complexities, whose results confirm the accuracy of the proposed scheme and illustrate the attenuation effect on the seismic energy during its propagation in the medium. Our results demonstrate that the numerical scheme employed can be used to extrapolate wavefields stably for even larger time steps than those usually used by traditional finite‐difference schemes.