<p>This paper proposes the application of the λ-Neumann Monte Carlo method to estimate the statistical moments of the solution to the bending problem of a Timoshenko beam resting on a stochastic Winkler foundation. Uncertainty is modeled in the beam and foundation stiffness coefficients, via parameterized stochastic processes. The methodology proposed in this work differs from the usual approach, as it is developed based on theoretical results regarding the existence and uniqueness of the realizations. The consistency of the numerical approximations is ensured through rigorous analysis of the well-posedness of the underlying stochastic boundary value problem. This theoretical result guides the choice of approximation spaces via the Galerkin method, leading to realizations consistent with this result. Consequently, statistically consistent estimators of the expected value and variance of the transverse and angular displacement stochastic processes are obtained. The Monte Carlo simulation method is used to evaluate the performance of the proposed methodology, in two practical examples.</p>