One answer was given in the comments, but here is another one which might be closer to the "consider the average" suggestion. Namely, one can consider the average $M_u(0,r)$ as a function of $r$ and write a second-order linear differential equation for it (basically, Laplace's equation averaged over spheres). The solutions form a two-dimensional space, which is $\{c_1\log r+c_2:c_1,c_2\in\mathbb R\}$ in two dimensions and $\{c_1 r^{2-n}+c_2:c_1,c_2\in\mathbb R\}$ in higher dimensions. No element of this space can have distinct finite limits at $0$ and at $1$.

I'm afraid that what I'm about to say is not quite you what to hear: your Neumann PDE does not have a solution for arbitrary choices of $g$.

To see why, let's compute the integral of $g$ over the surface $\Gamma$:
$$ \int_\Gamma g \ dS = \int_\Gamma \frac{\partial u}{\partial n} \ dS=\int_\Omega \nabla^2 u \ dV =\int_\Omega 0\ dV = 0.$$
Therefore, *the average value of $g$ on the boundary surface $S$ must be zero*. If $g$ fails to satisfy this condition, your PDE cannot be solved.

For a similar reason, your definition of the Green's function must be modified. Indeed, if we integrate $\frac{\partial G}{\partial n}$ over the surface $\Gamma$, we find that
$$ \int_\Gamma \frac{\partial G(\vec r, \vec r')}{\partial n_{\vec r}}dS(\vec r) = \int_\Omega \nabla_{\vec r}^2 G(\vec r,\vec r') \ dV(\vec r) = \int_V \delta(\vec r-\vec r') \ dV(\vec r) = 1.$$
So it is inconsistent to set $\frac{\partial G}{\partial n}$ equal to zero on the boundary $\Gamma$.

Fortunately, all is not lost! We can redefine the Green's function $G$ so that it satisfies
$$ \begin{cases} \nabla_{\vec r}^2 G(\vec r,\vec r') = \delta(\vec r - \vec r') &{\rm \ \ on \ \ } \Omega ,
\\
\frac{\partial G(\vec r, \vec r')}{\partial n_{\vec r}}=\frac 1 A & {\rm \ \ on \ \ } \Gamma, \end{cases}
$$
where $A = \int_\Gamma dS$ is the area of the boundary surface $\Gamma$.

Now, Green's identity states that
\begin{multline}\int_\Omega \left( u(\vec r) \nabla_{\vec r}^2 G(\vec r , \vec r') - G(\vec r,\vec r') \nabla_{\vec r} u(\vec r) \right) \ dV(\vec r) \\ = \int_\Gamma \left( u(\vec r) \frac{\partial G(\vec r,\vec r')}{\partial n_{\vec r}} - G(\vec r,\vec r') \frac{\partial u(\vec r)}{\partial n_{\vec r}}\right) dS(\vec r).\end{multline}
Plugging in my new definition of $G(\vec r, \vec r')$, we see that any $u(\vec r)$ satisfying your PDE must satisfy
we can immediately see that any solution to the PDE must satisfy
$$ u(\vec r') = - \int_\Gamma G(\vec r,\vec r')g(\vec r) \ dS(\vec r) + c, \ \ \ \ \ \ (\ast )$$
where the constant $c$ is equal to $\frac 1 A \int_\Gamma u(\vec r) dS(\vec r)$, the average value of $u$ on $\Gamma$.

Having redefined the Green's function, I'll give you an explicit expression in the case where $\Omega$ is a two-dimensional circular disk of radius $1$. Here it is:
$$ G(\vec r, \vec r') = \tfrac 1 {2\pi} \ln | \vec r -\vec r' | + \tfrac 1 {2\pi} \ln |\vec r - \vec r''|,$$
where $\vec r'' = \vec r'/|\vec r'|^2$ is the image of $r'$ under an inversion about the unit circle. It is similar to the Dirichlet Green's function, expect that we have a plus sign in front of the "image" term instead of a minus sign.

I believe that if you plug this Green's function into $(\ast )$ (with the constant $c$ chosen arbitrarily), you do indeed get a solution to your PDE, provided that $g$ obeys the consistency condition $\int_\Gamma g \ dS = 0$. This is stated in these lecture notes from my university, and also in Riley, Hobson and Bence, though I would love to see a rigorous proof! I wonder if anyone knows a good reference?

## Best Answer

Observe first that the problem is invariant under rotations, therefore solutions $u$ of the given problem are of the form $u(x)=f(|x|)$. Since \begin{eqnarray} \Delta u(x)&=&\sum_{i=1}^n\frac{\partial^2u}{\partial x_i^2}(x)=\sum_{i=1}^n\frac{\partial}{\partial x_i}\left(\frac{x_i}{|x|}f'(|x|)\right)=\sum_{i=1}^n\left[\frac{x_i^2}{|x|^2}f''(|x|)+\left(\frac{1}{|x|}-\frac{x_i^2}{|x|^3}\right)f'(|x|)\right]\cr &=&f''(|x|)+\frac{n-1}{|x|}f'(|x|), \end{eqnarray} we have to solve the ODE $$\tag{1} f''(t)+\frac{n-1}{t}f'(t)=0. $$ Integrating (1) we get $$ f(t)= \begin{cases} at^{2-n}+b & \text{ if } n \ge 3\\ a\ln t+b & \text{ if } n=2 \end{cases}, $$ where $a,b$ are real constants. Because of the conditions $f(1)=1$ and $\lim_{t \to \infty}f(t)=0$, there is no such $f$ for $n=2$. For $n \ge 3$, the solution is given by $$ f(t)=\frac{1}{t^{n-2}} \quad \forall t \ge 1. $$ Hence $$ u(x)=\frac{1}{|x|^{n-2}}\quad x \in \overline{\Omega}. $$