Seismic wave propagation in a viscoelastic media can be described by a system of integro-differential equations. The solution of such equations requires special methods when using finite-difference techniques in the time domain. In the frequency domain, the integral terms are represented by complex elastic parameters. This paper presents an efficient algorithm for viscoelastic modelling based on the integral Laguerre transform for the approximation of temporal derivatives and for the calculation of convolution integrals. For the calculation of spatial derivatives, it is possible to use various methods: finite-difference and finite-element techniques, spectral and pseudo-spectral methods. We then obtain a system of algebraic equations with a matrix independent of the parameter m , i.e. the degree of the Laguerre polynomials. In this case, only the right-hand side of the system has recurrent dependence on the parameter m, which is an analogue of the temporal frequency in the frequency domain. The obtained system with a large number of right-hand sides can be solved using fast methods, where the matrix is transformed only once, as opposed to the frequency-domain approach, when the matrix is transformed for each temporal frequency.