A sub-supersolution approach for Neumann boundary value problems with gradient dependence