A planarly-stratified multilayer is considered with material properties depending on the Cartesian coordinate normal to the layers. Upon the assumption that the time dependence be harmonic, the equations of motion and the constitutive equations (of a viscoelastic solid) are given the form of a first-order system of ordinary differential equations. The propagator of the whole multilayer is determined and hence the reflection and transmission matrices are obtained for different boundary conditions. Next a new algorithm, which avoids some drawbacks of other procedures, is outlined. The reflection and transmission matrices of the multilayer are determined by recursive relations via the matrices, of reflection and transmission, associated with the single layers.