In this paper, the Galerkin method is used to obtain numerical solutions to two-dimensional steady-state reaction-diffusion problems. Uncertainties in reaction and diffusion coefficients are modeled using parameterized stochastic processes. A stochastic version of the Lax–Milgram lemma is used in order to guarantee existence and uniqueness of the theoretical solutions. The space of approximate solutions is constructed by tensor product between finite dimensional deterministic functional spaces and spaces generated by chaos polynomials, derived from the Askey–Wiener scheme. Performance of the developed Galerkin scheme is evaluated by comparing first and second order moments and probability histograms obtained from approximate solutions with the corresponding estimates obtained via Monte Carlo simulation. Results for three example problems show very fast convergence of the approximate Galerkin solutions. Results also show that complete probability densities (histograms) of the responses are correctly approximated by the developed Galerkin basis.