Slopes in nature usually present layered characteristics, and its stability is susceptible to rainfall events. Considering that current analytical solutions are only suited to simulate the rainfall infiltration of double‐layered infinite unsaturated slopes, an analytical procedure is hence proposed in this study to tackle the consideration of multiple layers. The variable separation method and transfer matrix method are combined to derive the analytical solution of pore water pressure (PWP) for simulating rainfall infiltration in layered infinite unsaturated slopes. After having validated the proposed model and analytical solutions by comparing with existing literature and numerical simulation, the closed‐form solution of PWP is incorporated into the limit equilibrium for assessing slope stability. A three‐layer slope is selected as an example for further discussion. PWP distribution and factor of safety are calculated, considering the effects of saturated hydraulic conductivity and thickness of the upper layer, intensity of antecedent and subsequent rainfall, and varied soil unit weight along the depth. The slope stability subjected to rainfall effects is consistent with the variation of PWP. The proposed analytical solutions provide a simple and practical avenue for computing PWP distribution and evaluating the stability of multi‐layered slopes under rainfall conditions, which can also serve as a benchmark for numerical solutions.