Solution (source code)

= Solution

Put $\mathbf x=r\boldsymbol\omega$, with $|\boldsymbol\omega|=1$. The <Jacobian determinant> in <spherical coordinates> gives $dS=r^2d\omega=r^2\sin\theta\,d\theta\,d\varphi$, so
$$
f(r)=\frac1{4\pi}\int_{S^2}\phi(r\boldsymbol\omega)\,d\omega.
$$
Differentiation under this fixed-domain integral yields
$$
f'(r)=\frac1{4\pi}\int_{S^2}\nabla\phi(r\boldsymbol\omega)\cdot\boldsymbol\omega\,d\omega
=\frac1{4\pi r^2}\int_{\partial V_r}\frac{\partial\phi}{\partial n}\,dS
=\frac1{4\pi r^2}\int_{V_r}\nabla^2\phi\,dV=0.
$$
The last step is the <divergence theorem> and the <Laplace equation>. Hence $f$ is constant on radii for which the closed ball lies inside $V$. Continuity at the origin gives $\lim_{r\to0}f(r)=\phi(0)$, proving the <mean value property for harmonic functions>:
$$
\boxed{f(r)=\phi(0).}
$$