We solve the non-local convection equations. The solutions for four model problems are compared with results of generalized smoothed particle hydrodynamics (GSPH) simulations. In each case we test two closure schemes: (1) where third moments are defined by the diffusion approximation; and (2) where the full third moment equations are used and fourth moments are defined by a modified form of the quasinormal approximation. In overshooting models, the convective flux becomes negative shortly after the stability boundary. The negative amplitude remains small, and the temperature gradient in the overshooting zone has nearly the radiative value. Turbulent velocities decay by a factor of e after (0.5-1.5)l, depending on the model, where l is the mixing length. Turbulent viscosity is more important than negative buoyancy in decelerating overshooting fluid blobs. These predictions are consistent with helioseismology. The equations and the GSPH code use the same physical approximations, so it was anticipated that, if the closures for high-order moments are accurate enough, solutions for the low-order moments will automatically agree with the GSPH results. Such internal consistency holds approximately. Unexpectedly, however, the best second moments are found with the first closure scheme, and the best third moments are obtained with the second. The relationship among moments from the solution of the equations is not the same as the relationship found by GSPH simulations. In particular, if we use the best fourth moment closure model suggested by Paper II, we cannot obtain steady-state solutions, but adding a sort of diffusion for stability helps.