Mean-value formula for inhomogeneous harmonic functions I am working on Evans' PDE textbook problems, but I am stuck with the following problem about modification of the proof of the mean-value formula for harmonic functions. I cannot really see how to derive the second term of RHS of the formula below. I would appreciate it if someone could help me derive this formula. 

Modify the proof of the mean-value formulas to show for $n\ge 3$ that 
  $$
u(0)=\frac{1}{V(\partial B(0,r))}\int_{\partial B(0,r)}g \, dS
+\frac{1}{n(n-2)\alpha(n)}\int_{B(0,r)}\left[\frac{1}{|x|^{n-2}}-\frac{1}{r^{n-2}}\right]f \,dx,
$$
  provided 
  $$
\begin{cases}-\Delta u=f & \text{in } B^0(0,r) \\ \quad \, \, \, u=g & \text{on } \partial B(0,r).\end{cases} 
$$

 A: $\def\Vol{\operatorname{Vol}}$As Evans does in his proof of the mean value formula, define $\phi\colon [0,r] \to \mathbb R$ by 
\[
  \phi(s) = \frac 1{\Vol(\partial B_s)} \int_{\partial B_s} u(x)\, dS(x)
   = \frac 1{\Vol(\partial B_1)} \int_{\partial B_1} u(sx)\, dS(x)
\]
We have
\begin{align*}
  \phi'(s) &= \frac 1{\Vol(\partial B_1)} \int_{\partial B_1} Du(sx)x\, dS(x)\\
    &= \frac 1{n\alpha(n)s^{n-1}} \int_{\partial B_s} Du(x)\, \frac xs\, dS(x)\\
    &= \frac 1{n\alpha(n)s^{n-1}} \int_{\partial B_s} Du(x)\, \nu_{B_s}(x)\, dS(x)\\
    &= \frac 1{n\alpha(n)s^{n-1}} \int_{\partial B_s} \frac{\partial u}{\partial \nu}(x)\,dS(x)\\
    &= \frac{1}{n\alpha(n)s^{n-1}} \int_{B_s} \Delta u(x)\, dx\\
    &= -\frac 1{n\alpha(n)s^{n-1}}\int_{B_s}f(x)\, dx
\end{align*}
As $\phi$ is differentiable on $(0,r)$ and continuous on $[0,r]$, we have
\begin{align*}
  \phi(r) -\phi(0) &= \int_0^r \phi'(s)\, ds\\
       &= -\int_0^r \frac 1{n\alpha(n)s^{n-1}} \int_{B_s} f(x)\, dx\; ds\\
       &= -\frac 1{n\alpha(n)} \int_0^r s^{1-n}\int_0^s \int_{\partial B_\sigma} f(x)\, dS(x)\; d\sigma\; ds\\
       &= \frac 1{n\alpha(n)} \int_0^r \int_\sigma^r s^{1-n}\int_{\partial B_\sigma} f(x)\, dS(x)\; ds\; d\sigma\\
       &= -\frac 1{n(n-2)\alpha(n)} \int_0^r \int_{\partial B_\sigma} f(x)\, dS(x)\; d\sigma\\
       &= -\frac 1{n(n-2)\alpha(n)}\int_0^r \int_{\partial B_\sigma}\left(\frac 1{\sigma^{n-2}} - \frac 1{r^{n-2}}\right) f(x)\, dS(x)\;d\sigma\\
       &= -\frac 1{n(n-2)\alpha(n)}\int_0^r \int_{\partial B_\sigma}\left(\frac 1{|x|^{n-2}} - \frac 1{r^{n-2}}\right) f(x)\, dS(x)\;d\sigma\\
       &= -\frac 1{n(n-2)\alpha(n)}\int_{B_r} \left(\frac 1{|x|^{n-2}} - \frac 1{r^{n-2}}\right) f(x)\, dx
\end{align*}
Now, as $\phi(0) = u(0)$ and
\[ \phi(r) = \frac 1{\Vol(\partial B_r)}\int_{\partial B_r} g(x)\, dS(x) \]
the desired result follows.
A: Alternatively, look at the second term on the right side. It's the fundamental solution, shifted down to be $0$ on the boundary of $B_r$, integrated against $\Delta u$. Heuristically, integrate by parts and use that the Laplacian of the fundamental solution is $\delta_0$ to get the desired formula.
