Derive mean value property of harmonic function using this method We can view harmonic function $u(x,y,z)$ as a wave function that does not depend on time $t$. Let $\overline{u}(r)$ be its average value on the circle with radius $r$.
Assume we have proved that it satisfies $0=(\overline{u})_{rr}+2\frac{1}{r}(\overline{u})_r$.
I think it remains to show that $\overline{u}$ is constant so that $\overline{u}$ equals $u(0)$, which is equivalent to show $(\overline{u})_r=0$. We can see $0$ is one solution of the ODE, but I cannot disprove that others are not.
Any hints would be appreciated!
 A: Let $v=\overline{u}$.
$v_{rr}+\frac{2}{r} v_r=0$, so $\frac{d}{dr} \left ( e^{2 \log(r)} v_r \right ) =0$, so $r^2 v_r=C$, so $v_r=\frac{C}{r^2}$, so $v=C/r+A$. 
This is the general solution to the equation on any interval not containing $r=0$. This means that if you delete the origin, then there is a nonconstant radial function which satisfies the ODE. However, this function is singular at the origin. If you assume the function is not singular at the origin, then you must have $C=0$ in which case $v=A$ as desired.
Note that this singular solution does not furnish any contradiction to the mean value theorem as it is correctly formulated (the domain of the harmonic function must contain the entire ball whose boundary is the sphere in question).
A: METHODOLOGY $1$:
I thought it might be instructive to show that if $\bar u(a)$ is the average value, taken over a sphere of radius $a$, of a harmonic function $u$, then $\bar u(a)$ is constant, independent of the radius of the sphere.
Let $V$ be a sphere of radius $a$ centered at $\vec r$ and let $S$ be its boundary described by $|\vec r-\vec r'|=a$.  The outer unit normal on $S$ at $\vec r'$ is given by $\hat n'=-\frac1a (\vec r-\vec r')$
Let $G(\vec r|\vec r')=\frac{1}{4\pi |\vec r-\vec r'|}$ be the Green Function (or alternatively the Green's Function) for Laplace' Equation so that $G$ satisfies
$$\nabla^2 G(\vec r|\vec r')=\delta(\vec r-\vec r') \tag 1$$  
Note that on the sphere, we have 
$$\left.\left(G(\vec r|\vec r')\right)\right|_{S}=\frac{1}{4\pi a} \tag 2$$ 
and 
$$\left.\left(\nabla' G(\vec r|\vec r')\right)\right|_{S}=- \frac{\hat n'}{4\pi a^2} \tag 3$$
Using Green's Identity, we can write 
$$\int_{V}\left(u(\vec r')\nabla'^2 G(\vec r|\vec r')-G(\vec r|\vec r')\nabla'^2 u(\vec r')\right)dV'=\oint_{S}\left(u(\vec r')\nabla' G(\vec r|\vec r')-G(\vec r|\vec r')\nabla' u(\vec r')\right)\cdot \hat n'dS' \tag 4$$
Since $u$ is harmonic, and since $G$ satisfies $(1)$, the left-hand side of $(4)$ becomes 
$\int_{V}\left(u(\vec r')\nabla'^2 G(\vec r|\vec r')-G(\vec r|\vec r')\nabla'^2 u(\vec r')\right)dV'=-u(\vec r) \tag 5$.  
We next exploit the fact that $G$ and its normal derivative are constant on the sphere as given by $(2)$ and $(3)$.  Then, the first term on the right-hand side of $(4)$ can be simplified as 
$$\begin{align}
\oint_{S}u(\vec r')\nabla' G(\vec r|\vec r')\cdot \hat n'dS' &=-\frac{1}{4\pi a^2}\oint_S u(\vec r')\,dS' \tag 6
\end{align}$$
while the second term on the right-hand side of $(1)$ becomes
$$\begin{align}
\oint_{S}G(\vec r|\vec r')\nabla' u(\vec r')\cdot \hat n'dS'&=\frac{1}{4\pi a}\oint_S \nabla u(\vec r')\cdot \hat n'\,dS'\\\\
&=\frac{1}{4\pi a}\int_V \nabla'^2 u(\vec r')\,dV'\\\\
&=0 \tag 7
\end{align}$$
Therefore, putting $(4)-(7)$ together, we find that 
$$u(\vec r)=\frac{1}{4\pi a^2}\oint_S u(\vec r')\,dS'$$
which states that the average of the harmonic function, taken over a sphere of radius $a$ that is centered at $\vec r$, is equal to $u(\vec r)$, the value of $u$ at the center of the sphere.  
Inasmuch as the value of $u(\vec r)$ for any $\vec r$ depends in no way on the radius of a sphere that we constructed, then the average cannot depend on this radius.  It is therefore independent of both $a$ and coordinate variables rendering it constant, as was to be shown.

METHODOLOGY $2$:
Here, we will show that if $\bar u(a)$ is the average value, taken over a sphere of radius $a$, of a harmonic function $u$, and $\bar u(a)$ satisfies the ODE
$$\frac{1}{a^2}\frac{d}{da}\left(a^2\frac{d\bar u(a)}{da}\right)=0 \tag{8}$$
then $\bar u$ is a constant, independent of the radius $a$.
We assume that $(8)$ is true.  If so, then we have
$$\frac{d}{da}\left(a^2\frac{d\bar u(a)}{da}\right)=0\implies \frac{d\bar u(a)}{da}=-B/a^2 \implies \bar u(a)=A+Ba^{-1}$$
where $A$ and $B$ are integration constants.  
Next, we note that $\bar u(a)= \frac{1}{4\pi}\int_{\Omega} u(\vec r') \,d\Omega'$, where $\Omega$ is the solid angle of the sphere.  Since $u$ is bounded inside the sphere, then 
$$\lim_{a\to 0}\bar u(a)=\frac{1}{4\pi}\int_{\Omega}\,\lim_{a\to 0} u(\vec r') \,d\Omega'$$
must be finite.  We conclude, therefore, that we must have $B=0$ and that $\bar u$ is independent on the radius $a$, as was to be shown.

