A systematic and efficient algorithm, the generalized reflection and transmission matrix method, has been developed for wave propagation in stratified fluid-saturated poroelastic half-space. The proposed method has the advantage of computational efficiency and numerical stability for high frequencies and large layer thickness. A wide class of seismic sources, ranging from a single-body force to double couples, is introduced by utilizing the moment tensor concept. In order to validate the proposed algorithm, we applied our formulation to calculate wave fields in a homogeneous poroelastic half-space. It is shown that the numerical results computed with the present approach agree well with those computed with the analytical solution. Numerical examples for a two-layer model subjected to various sources such as double couple, dipole, and explosive sources are provided. From the waveforms of surface displacements, the arrivals of transmitted and converted PS and SP waves at the interface of the two-layer model can be clearly observed. As expected, it is impossible to observe the arrivals of transmitted and transmitted and converted SP waves from the waveforms induced by fluid withdrawal.