In this paper, we construct two linearly independent response processes to the random Legendre differential equation on (−1,1)∪(1,3), consisting of Lp(Ω) convergent random power series around the regular‐singular point 1. A theorem on the existence and uniqueness of Lp(Ω) solution to the random Legendre differential equation on the intervals (−1,1) and (1,3) is obtained. The hypotheses assumed are simple: initial conditions in Lp(Ω) and random input A in L∞(Ω) (this is equivalent to A having absolute moments that grow at most exponentially). Thus, this paper extends the deterministic theory to a random framework. Uncertainty quantification for the solution stochastic process is performed by truncating the random series and taking limits in Lp(Ω). In the numerical experiments, we approximate its expectation and variance for certain forms of the differential equation. The reliability of our approach is compared with Monte Carlo simulations and generalized polynomial chaos expansions.