Solution (source code)

= Solution

For the spherical average $M(r)$, rescaling the sphere and applying the <divergence theorem> gives
$$
M(r)=\frac1{n\omega_n}\int_{S^{n-1}}u(y+r\theta)\,dS_\theta,\qquad
M'(r)=\frac1{n\omega_n r^{n-1}}\int_{\partial B_r(y)}\partial_\nu u
=\frac1{n\omega_n r^{n-1}}\int_{B_r(y)}\Delta u=0.
$$
By <continuity>, $M(r)\to u(y)$ as $r\downarrow0$. Integrating these spherical averages in radius proves the <mean value property for harmonic functions>:
$$
\boxed{u(y)=\frac1{n\omega_n r^{n-1}}\int_{\partial B_r(y)}u
=\frac1{\omega_n r^n}\int_{B_r(y)}u.}
$$
For the <strong maximum principle for harmonic functions>, suppose an interior point attains the maximum $m$ on a connected domain. The continuous nonnegative function $m-u$ has zero average on each sufficiently small ball around that point, so it vanishes there. The maximum set is therefore open as well as closed, and <connectedness> makes it the entire domain. Apply the same argument to $-u$ for a minimum.