We propose a new algorithm to compute the teleseismic P and S receiver function synthetics for a multilayered Cartesian structure with anisotropic flat layers. The algorithm is based on the first-order perturbation theory in which the layered background structure is assumed one-dimensional with isotropic variations in vertical direction. Anisotropic velocity perturbations acting as secondary sources constitute the heterogeneities in the medium. The total wavefield is solved using a convolutional type integral equation along with the Green's function of the one-dimensional reference medium extracted using the reflectivity method. The integral equation involves a five-fold integration in space and wavenumber domains. Four of these integrals are achieved analytically and the fifth integral, which is spatial integral in the vertical direction, is performed numerically for which the Born single scattering approximation greatly suffices.The proposed algorithm is demonstrated on some selected numerical examples adapted from published work in the literature. © 2013 Elsevier B.V.