Solution (source code)

= Solution

The continuous <harmonic function> has $m\le\phi(\mathbf x)\le M$ on the compact sphere. Integrate these inequalities over its area $4\pi r^2$ and use the <mean value property for harmonic functions>:
$$
4\pi r^2m\le\int_{\partial V_r}\phi\,dS\le4\pi r^2M,
\qquad \boxed{m\le\phi(0)\le M.}
$$
The centre value is an average with positive weights, so it lies between the smallest and largest boundary-sphere values.