Solution (source code)

= Solution

Fix $x\in\Omega$ and $R<\operatorname{dist}(x,\partial\Omega)$. The <mean value property for harmonic functions> is
$$
u(x)=\frac1{|\partial B_r|}\int_{\partial B_r(x)}u\,dS
=\frac1{|B_r|}\int_{B_r(x)}u\,dy
\qquad(0<r\leq R).
$$
To prove the spherical identity, translate $x$ to the origin and put $F(r)=|S^{n-1}|^{-1}\int_{S^{n-1}}u(r\theta)\,dS_\theta$. The <divergence theorem> gives
$$
F'(r)=\frac1{|S^{n-1}|r^{n-1}}\int_{B_r}\Delta u\,dy=0.
$$
Hence $F(r)=\lim_{s\downarrow0}F(s)=u(0)$. Integrating the spherical identity in polar coordinates gives the ball identity.