Numerical simulation of wave propagation phenomena is a relevant subject in many different fields of sciences and engineering, like geophysics, ocean acoustics, material science and biomedical engineering. As a matter of fact, an effective modelling of wave propagation phenomena in real media requires to take into account several non-trivial but fundamental physical aspects, such as energy dissipation due to viscous effects and heterogeneous spatial distribution of material properties, which pose additional challenges from both a mathematical modelling and numerical point of view. This talk is concerned with the aforementioned framework and we in particular focus on the numerical aspects of the problem, thereby presenting a new numerical approach aimed at simulating the propagation of waves in viscoelastic and heterogeneous materials. Specifically, we consider a viscoelastic material behaviour described by means of the well-known Kelvin-Voigt rheology, whereby the stress tensor appearing in the Cauchy momentum equations is composed by an elastic term depending on the displacement as well as a viscous term depending on the velocity. Moreover, the material heterogeneity has been taken into account by considering a single domain (body) partitioned in several layers, each of them associated to a different material that can be either solid or fluid and characterised by a specific set of coefficients defining its physical properties. The system of Partial Differential Equations has been reduced from second order to first order in time by doubling the solution variables, that are the displacement and velocity fields. The model has been solved numerically in space by employing a Galerkin Spectral approach which approximate the displacement and velocity using Chebyshev polynomials defined on each subdomain. The continuity of the solution across the material interfaces is obtained by imposing weak continuity of the basis functions. The discretisation in space yields a system of Ordinary Differential Equations that are solved by means of standard numerical time schemes.