The paper presents an efficient algorithm, based on the application of the finite integral Fourier–Bessel transforms with respect to spatial coordinate and the spectral Laguerre method for approximation of temporal derivatives as applied to the problem of seismic wave propagation in the heterogeneous viscoelastic medium. Such an approach allows reducing the solution of the original problem to a series of one-dimensional problems, which are then numerically solved. As a elaxation function the generalized model for the standard linear solid rheology is considered. Some examples of the calculation of synthetic seismograms for different models of viscoelastic media are presented.