A function satisfying the mean value property is harmonic Here is the problem. 

I know that if $u$ is harmonic the equation holds, but I don't know how to prove it from the other direction.
 A: First of all, we note that by a change of variables
$$\begin{align*} u(x) &= \frac{1}{n \omega_n r^{n-1}} \int_{\partial B(x,r)} u(z) \, dS(z) \\ &= \frac{1}{\sigma_{n-1}} \int_{S_{n-1}} u(x+r \xi) \, dS(\xi)  \tag{1} \end{align*}$$
where $\sigma_{n-1}$ denotes the volume of the ($(n-1)$-dimensional) volume of the sphere $S_{n-1} = \partial B_{\mathbb{R}^n}(0,1)$. Now suppose that $\Delta u(x)<0$ for some $x$. Then it follows from the continuity of $\Delta u$ that we can choose $r_0>0$ such that $\Delta u(y)<0$ for all $y \in B(x,r_0)$. From the Lemma below, we conclude $$\frac{d}{dr} \int_{S_{n-1}} u(x+r \xi) \, dS(\xi) < 0$$ and therefore $$\int_{S_{n-1}} u(x+r \xi) \, dS(\xi) > \int_{S_{n-1}} u(x+r_0 \xi) \, dS(\xi) \quad \text{for all} \,  0<r<r_0.$$ Consequently, we obtain $$\begin{align*} \sigma_{n-1} u(x) &= \lim_{r \to 0} \int_{S_{n-1}} u(x+r \xi) \, dS(\xi) \\ &> \int_{S_{n-1}} u(x+r_0 \xi) \, dS(\xi). \end{align*}$$ Obviously, this contradicts $(1)$. A similar argumentation yields a contradiction if $\Delta u(x)>0$ for some $x$.

Lemma  $$\frac{d}{dr} \int_{S_{n-1}} u(x+r \xi) \, dS(\xi) = \frac{1}{r^{n-1}} \int_{B(x,r)} \Delta u(z) \, dz.$$

Proof: This equality follows from Gauß divergence theorem: $$\begin{align*} \frac{d}{dr} \int_{S_{n-1}} u(x+r \xi) \, dS(\xi) &= \int \xi \cdot \nabla u(x+r \xi) \, dS(\xi) \\ &= \frac{1}{r^{n-1}} \int_{\partial B(x,r)} (\nabla u(x+y) \mid \nu(y)) \, dS(y) \\ &= \frac{1}{r^{n-1}} \int_{B(x,r)} \underbrace{\text{div} \nabla}_{\Delta} u(z) \, dz.\end{align*}$$ Here $\nu(y)$ denotes the outer pointing normal.
