This paper deals with the nonlinear two point boundary value problem y″ = f( x, y, y′, R 1,…, R n ), x 0 < x < x f S 1 y( x 0) + S 2 y′( x 0) = S 3, S 4 y( x f ) + S 5 y′( x f ) = S 6 where R 1,…, R n , S 1,…, S 6 are bounded continuous random variables. An approximate probability distribution function for y( x) is constructed by numerical integration of a set of related deterministic problems. Two distinct methods are described, and in each case convergence of the approximate distribution function to the actual distribution function is established. Primary attention is placed on problems with two random variables, but various generalizations are noted. As an example, a nonlinear one-dimensional heat conduction problem containing one or two random variables is studied in some detail.