A numerical model is proposed to compute the eigenmodes and the forcedresponse of multilayered elastic spheres. The main idea is to describe analyti-cally the problem along the angular coordinates with spherical harmonics andto discretize the radial direction with one-dimensional finite elements. Theproper test function must be carefully chosen so that both vector and tensorspherical harmonics orthogonality relationships can be used. The proposedapproach yields a general one-dimensional formulation with a fully analyticaldescription of the angular behaviour, suitable for any interpolating technique.A linear eigenvalue problem, simple and fast to solve, is then obtained. Theeigensolutions are the spheroidal and torsional modes. They are favourablycompared with literature results for a homogeneous sphere. The eigensolu-tions are superposed to compute explicitly the forced response. The latter isused to reconstruct the propagation of surfaces waves. In particular, the col-limation of a Rayleigh wave (non-diffracted surface wave propagating with aquasi-constant width) excited by a line source in a homogeneous sphere is re-covered with the model. Based on the vibration eigenmodes, a modal analysisshows that such a wave is a superposition of fundamental spheroidal modeswith a displacement confined at the equator of the sphere. These modes arethe so-called Rayleigh modes, of sectoral type and high polar wavenumbers.When a thin viscoelastic coating is added to the sphere, the Rayleigh modebehaviour is recovered in a limited frequency range, allowing the generationof a collimating wave at the interface between the sphere and the coating