Proving the Liouville Theorem I am trying to prove the Liouville Theorem using the Mean Value Property for Harmonic Functions. The question is:

Let $u$ be a Harmonic Function and $|u|$ is bounded. Prove that $u$ is constant. Interpret this when $u$ represents the temperature at a point $\textbf x$ in the plane. 

Here is my $\textit{attempt} $ at a solution.
We begin by showing:
$$|\partial_{x_i} u(\textbf x_0)|\leq \frac{3|n_i|M}{r}$$
Noting that $\partial_{x_i} u(\textbf x_0)$ is Harmonic.
Applying the MVP:
$$|\partial_{x_i} u(\textbf x_0)| = \bigg|\frac{3}{4\pi r^3}\iiint_{B(\textbf x_0,r)}\partial_{x_i} u(\vec x) d\vec x\bigg|$$
Using the Divergence Theorem:
$$\bigg|\frac{3}{4\pi r^3}\iiint_{B(\textbf x_0,r)}\partial_{x_i} u(\vec x) d\vec x\bigg| = \frac{3}{4\pi r^3}\bigg|\iint_{\partial B(\textbf x_0, r)} u(\vec x)\cdot n_i dS_{\vec x}\bigg|$$ 
Where $n_i$ denotes the $i$-th component of the unit normal vector. Now using the Cauchy Inequality:
$$\frac{3}{4\pi r^3}\bigg|\iint_{\partial B(\textbf x_0, r)} u(\vec x)\cdot n_i dS_{\vec x}\bigg| \leq \frac{3}{4\pi r^3}\iint_{\partial B(\textbf x_0,r)}| u(\vec x)|\cdot|n_i|dS_{\vec x}$$
The above yields:
$$\frac{3}{4\pi r^3}\iint_{\partial B(\textbf x_0,r)}| u(\vec x)|\cdot|n_i|dS_{\vec x}\leq \frac{3}{4\pi r^3}\cdot 4\pi r^2|n_i|\cdot M = \frac{3|n_i|\cdot M}{r}$$
Now, applying this to:
$$\bigg|\sum_{i=1}^3\partial x_i u(\textbf x_0)\bigg|\leq \frac{C M}{r}$$
Where $C$ is some constant.
Now, since $u$ and $\partial_{x_i} u$ are both Harmonic on $\mathbb R^3$, taking $r\rightarrow\infty$ yields that:
$$|\nabla u(\textbf x_0)| =0$$
Thus we conclude that $u$ is a constant.
Am I missing anything? Also if someone could explain how to interpret this as a temperature it would be appreciated. Thanks!
 A: Since the argument seems all right to me, I'll just add a few words on two generalizations (as partially requested in the comments). 

First, the argument actually shows that if $u$ has polynomial growth, i.e. $|u(x)| \lesssim (1+|x|)^k$ for some $k$, then $|\nabla u| \lesssim (1+|x|)^{k-1}$. Iterating this, we get 
$$
|u(x)| \lesssim (1+|x|)^k
\quad \Rightarrow \quad 
D^{k+1} u \equiv 0 
\quad \Rightarrow \quad
u \text{ is a polynomial of degree } \le k.
$$

Second, it is enough to assume that $u$ is bounded from below (or from above). 
Without loss of generality, assume that $u \ge 0$ in $\mathbb R^n$ (otherwise consider $\pm u + c$). Take any two points $x,y \in \mathbb R^n$ and denote $d = |x-y|$. For any $r>0$ we have $B(x,r) \subseteq B(y,r+d)$ and so
$$ \omega_n r^n u(x) = \int_{B(x,r)} u \le \int_{B(y,r+d)}u = \omega_n(r+d)^n u(y) $$
by mean value property. Thus, we obtained the inequality 
$$ u(x) \le \left( 1 + \frac d r \right)^n u(y) $$
valid for any $r>0$, in consequence $u(x) \le u(y)$. Since $x,y$ are arbitrary (and in particular the order doesn't matter), $u$ is constant. 
