Let $\Omega$ be a bounded smooth domain with $$-\Delta u + ku = f$$ $$\partial_\nu u|_{\partial\Omega} = 0$$ where $k > 0$ is a constant and $f \in L^\infty(\Omega)$. It follows that $u \in H^2(\Omega)\cap C^{1,\alpha}(\bar \Omega)$.

I want an estimate on $\lVert u \rVert_{C^{1,\alpha}(\bar \Omega)}$ (or at least an estimate on the $L^\infty$ norm of the gradient of $u$) in terms of the data but I cannot find any literature for the Neumann problem giving an a priori estimate.

Does anyone know this estimate? I'd rather keep dimension as general as possible so I don't want to use Sobolev embeddings to obtain this.