A method is presented to study the three-dimensional quasi-static response of a multi-layered poroelastic half-space with compressible constituents. The system under consideration consists of N layers of different thickness and material properties overlying a homogeneous half-space. Fourier expansion, Laplace transforms and Hankel transforms with respect to the circumferential, time and radial coordinates, respectively, are used in the formulation. Laplace-Hankel transforms of displacements and pore pressure at layer interfaces are considered as the basic unknowns. Exact stiffness matrices describing the relationship between generalized displacement and force vectors of a finite layer and a half-space are derived explicitly in the transform space. The global stiffness matrix of a layered system is assembled by considering the continuity of tractions and fluid flow at layer interfaces. The time histories of displacements, stresses and pore pressure are obtained by solving the stiffness equation system for discrete values of Laplace and Hankel transform parameters, and using numerical quadrature schemes for Laplace and Hankel transform inversions. Selected numerical results for different layered systems are presented to portray the influence of layering and poroelastic material properties. The advantage of the present method is that for an N-layered system, it yields a numerically stable symmetric stiffness matrix of order 4N × 4N when compared to the unsymmetric and numerically unstable coefficient matrix of order 8N × 8N associated with the conventional method based on the determination of layer arbitrary coefficients.