A numerical method to compute propagation constants and mode shapes of elastic waves in layered piezoelectric-piezomagnetic composites, potentially deposited on a substrate, is described. The basic feature of the method is the expansion of particle displacement, stress fields, electric and magnetic potentials in each layer on different polynomial bases: Legendre for a layer of finite thickness and Laguerre for the semi-infinite substrate. The exponential convergence rate of the method for propagation of Love waves is numerically verified. The main advantage of the method is to directly determine complex wave numbers for a given frequency via a matricial eigenvalue problem, in a way that no transcendental equation has to be solved. Results are presented and the method is discussed.