For elliptic boundary value problems of the form −ΔU+F(x, U, U x )=0 on Ω,B[U]=0 on ϖΩ, with a nonlinearityF growing at most quadratically with respect to the gradientU x and with a mixed-type linear boundary opeatorB, a numerical method is presented which can be used to prove the existence of a solution within a “close”H 1,4(Ω)-neighborhood of some approximate solution ω∈H 2(Ω) satisfying the boundary condition, provided that the defect-norm ∥−Δω +F(·, ω, ω x )∥2 is sufficiently small and, moreover, the linearization of the given problem at ω leads to an invertible operatorL. The main tools are explicit Sobolev imbeddings and eigenvalue bounds forL or forL*L. All kinds of monotonicity or inverse-positivity assumptions are avoided.