<p>This paper proposes a general computational method for predicting the plane-strain consolidation behavior of multi-layered saturated soils with distributed drainage boundary. The method first transforms the consolidation problem into the Laplace domain. Subsequently, employing the Fourier cosine series expansion and its orthogonality, the distributed drainage boundary condition is directly converted into a linear system involving the solution coefficients, obtaining the corresponding analytical expressions in the Laplace domain. Finally, the numerical inverse Laplace transform technique yields the final solution. Convergence analysis is conducted to determine the minimum number of summation terms required for an accurate solution under typical conditions. The analysis reveals that increasing the pave rate, thickness-width ratio, anisotropy coefficient, and relative permeability, as well as reducing the relative compressibility, reduces the number of series terms required for convergence, thereby improving computational efficiency. Parametric analysis further revealed that compared with the full drainage boundary, the excess pore-water pressure under the distributed drainage boundary is significantly lower in the early stage of consolidation, particularly near the drainage channel. This difference gradually diminishes in the middle and late stages. Additionally, variations in the topsoil parameters have a more pronounced impact on consolidation characteristics than those in the lower soil layer. Therefore, in engineering practice, the design of distributed drainage boundaries should prioritize the characteristics of the topsoil.</p>