$C^0$ estimate for solutions of elliptic PDE with Neumann BC I am interested in a reference for (or counterexample to) a particular
$C^0$ estimate for solutions of the Laplace equation with Neumann
boundary conditions.  More precisely, let $(M,g)$ be a $C^\infty$
compact Riemannian manifold with non-empty smooth boundary $\partial M$. Let
$\Delta_g=-div_g\nabla$ denote the Laplace operator and $\nu$ the
outward pointing unit normal on $\partial M$. 
Consider (strong)
solutions $u$ to the following system
\begin{cases}
\Delta_gu=f&\text{ in }M\\
\partial_\nu u=h&\text{ on }\partial M
\end{cases}
where $f\in C^\infty(M)$ and $h\in C^\infty(\partial M)$ satisfy
$\int_M fdvol_g=\int_{\partial M}hd\sigma_g$ (the smoothness 
condition is for simplicity but the later condition is necessary
for the problem to be well-posed).  Here $dvol_g$ denotes the Riemannian
volume form and $d\sigma_g$ denotes its pullback under the inclusion
$\partial M\to M$. To specify $u$ uniquely, let
us also suppose that it has zero average value i.e. $\int_Mudvol_g=0$.
My question is this: does there exist $C>0$, independent of $f$ and $h$,
such that $\sup_M|u|\leq C(\sup_M|f|+\sup_{\partial M}|h|)$?
Most classical sources (Evans, Gilbarg+Trudinger, and Aubin) either only consider Dirichlet conditions (in which case the analogous estimate is essentially just the maximum principle) or only provide Schauder estimates for the Neumann problem -- this is not quite what I want.
Does anyone know of a resource for this kind of estimate or of a counterexample? 
 A: Yes, I think the estimate you propose is true. As a simple case take $f = 0$
and $M = B_1 \subset \mathbb{R}^n$. By adding a constant assume that $\inf_{B_1} u = 0$ and $\sup_{B_1}u = K,$ and note that these extrema are on the boundary. Note that $u(0)$ is closer to $0$ or $K$, say $u(0) > K/2$ without loss of generality. Then by the Harnack inequality we have $u > cK$ in $B_{1/2}$, so the function
$$c(n)K(|x|^{2-n} - 1)$$
is a lower barrier in the annulus $B_1-B_{1/2}$. In particular, at the point where $u$ takes its minimum we have $|h| > c(n)K$ which gives the estimate.
In a more general domain we can argue similarly; take an interior ball touching the point on the boundary where $u = 0$ and argue that $u$ is no larger than $C\sup|h|$ at its center, and similarly take a ball touching the point where $u = K$ and argue that $u$ is larger than $K - C\sup|h|$ at its center. If $K$ is much larger than $\sup|h|$ we contradict the Harnack inequality (with a constant depending on the geometry of $M$).
Finally, if $f$ is nonzero similar arguments should work since we still have the Harnack inequality (now depending on $\sup|f|$) and for a barrier we can take $|x|^{-\alpha}$ for $\alpha$ large depending on $\sup|f|$. 
A: Suppose for simplicity that $M$ is a compact n-dimensional submanifold of $\mathbb{R}^n$ with boundary.
Extend $f$ to all $\mathbb{R}^n$ so that it's $0$ outside $M$ and define the function $u_1(x) := \displaystyle \int_{\mathbb{R}^n}\phi(y-x)f(y)dy$ where $\phi$ is the fundamental solution to the Laplace equation. Then $$-\Delta u_1 = f$$ in $M$, and we have the following estimates:
$$\rVert u_1 \rVert_{C^0(M)} + \rVert \nabla u_1 \rVert_{C^0(M)} \leq C_1 \lVert f \rVert_{C^0(M)}$$
Now we find an integral representation of a solution to the homogeneous system using single layer potentials. Define $S (g) := \int_{\partial M} \phi(x-y) g(y) dy$ for $g \in C^0(\partial M)$. It is known that $S$ is a compact bounded linear operator on $C^{0}(\partial M)$ into $C^0(M)$ (check "The Laplace Equation" by Medkova, Proposition 6.7.1). It is also known that there exists an isomorphism $T$ on $C^0(\partial M)$ such that if $v = S(\psi_1)$ and $\psi_2 = T(\psi_1)$ then $\Delta v= 0$ on $M$ and $\frac{\partial v}{\partial \nu} = \psi_2$ on $\partial M$. And so define $u_2:= S\circ T^{-1} (h-\frac{\partial u_1}{\partial \nu})$, which will then solve
$$\Delta u_2 = 0 \quad \text{in} \  M, \qquad \frac{\partial u_2}{\partial \nu} = h-\frac{\partial u_1}{\partial \nu} \quad \text{in} \ \partial M$$ And so we have
$$\lVert u_2 \rVert_{C^0(M)} \leq \lVert S \rVert \lVert T^{-1} \rVert \left(\lVert h\rVert_{C^0(\partial M)}+\lVert \frac{\partial u_1}{\partial \nu} \rVert_{C^0(\partial M)} \right) \leq C_2(\lVert h\rVert_{C^0(\partial M)} + \lVert f\rVert_{C^0(M)})$$
Finally, let $C_3 = -\int_M (u_1+u_2)$; then $u:= u_1 +u_2 + C_3$ solves
$$-\Delta u = f \quad \text{in} \ M, \qquad \frac{\partial u}{\partial \nu} = h \quad \text{in}  \ \partial M, \qquad \int_M u = 0$$
and we get
$$\lVert u\rVert_{C^0( M)} \leq \lVert u_1\rVert_{C^0(M)} + \lVert u_2\rVert_{C^0( M)} + |C_3| \leq C \left( \lVert h\rVert_{C^0(\partial M)} + \lVert f\rVert_{C^0( M)} \right) + |C_3| $$
and the desired estimate follows.
