$\Delta u=\text{constant}$ in a domain with von Neumann boundary condition

Solution 1:

That PDE with that boundary values has the weak formulation $$\int_{\Omega}\nabla u\cdot\nabla\varphi-k\varphi+\int_{\partial\Omega}\varphi=0$$ for every $\varphi\in C^{\infty}(\overline{\Omega})$. With this weak formulation, and with the trace inequality, it can be solved in $H^1(\Omega)$ in classical ways, for example with Lax-Milgram theorem.

If $\Omega$ is the unit ball, the solution is $\frac{1}{2}|x|^2$.