We consider steady state unsaturated flow in bounded, randomly heterogeneous soils under the influence of random boundary and source terms. Our aim is to predict pressure heads and fluxes without resorting to Monte Carlo simulation, upscaling, or linearization of the constitutive relationship between unsaturated hydraulic conductivity and pressure head. We represent this relationship through Gardner's exponential model while treating its exponent α as a random constant and saturated hydraulic conductivity K_s as a spatially correlated random field. We linearize the steady state unsaturated flow equations by means of the Kirchhoff transformation and integrate them in probability space to obtain exact integro-differential equations for the conditional mean and variance-covariance of transformed pressure head and flux. After approximating these equations recursively to second order in the standard deviation σσ_Y of Y = ln K_s, we solve them by finite elements for superimposed mean uniform and divergent flows in the vertical plane, with and without conditioning on measured Y values. Comparison with Monte Carlo solutions demonstrates that whereas our nonlocal solution is nominally restricted to mildly nonuniform media with σ_Y² « 1, it yields remarkably accurate results for strongly nonuniform media with σ_Y² at least as large as 2. This accords well with a previous theoretical analysis, which shows that the solution may remain asymptotic for values of σ_Y² as large as 2.