We consider $u:U\to\mathbb{R}$ be a harmonic function. Now, I need to prove that for every $x_{0}\in U$, the function $\left( {0,{\rm{dist}}\left( {{x_0},\partial U} \right)} \right)\ni\rho \mapsto \dfrac{1}{{{\rho ^{n - 2}}}}\int_{{B_\rho }\left( {{x_0}} \right)} {{{\left| {Du} \right|}^2}dx} $
is non-decreasing.
We know that if $u$ is a harmonic function then $u$ be also a subharmonic function. But I am referring Example 2 in Section 8.6.2 in Partial Differential Equations Lawrence C.Evan with Monotonicity formula for harmonic functions . But I do not have an idea approaching the problem to solve because I confuse with them.
Thank you advance for your supports
As is mentioned in the comments, the proof you're looking is precisely the example in Evans that you've been referred to except with $p=2$. From what I understand, your issue is in the first line when Evans computes $ \frac d {dr} \big ( \frac 1 {r^{n-2}}\int_{B(0,r)} \vert Du \vert^2 d x\big ) $ and possibly integrating (10) over $B(0,r)$. The product rule for derivatives and polar coordinates (see Section C.3 Thm 4 (ii) in Evans) yield \begin{align} \frac d {dr} \bigg ( {r^{2-n}}\int_{B(0,r)} \vert Du \vert^2 d x\bigg ) &= (2-n)r^{1-n}\int_{B(0,r)} \vert Du \vert^2 d x + r^{2-n}\frac d {dr} \int_{B(0,r)} \vert Du \vert^2 d x \\ &= \frac{2-n}{r^{n-1}}\int_{B(0,r)} \vert Du \vert^2 d x + \frac{1}{r^{n-2}} \int_{\partial B(0,r)} \vert Du \vert^2 d S \tag{1} \end{align} which is the second line in Evans' computation. Setting $p=2$ in (10) in the same example gives \begin{align*} \sum_{i=1}^n \big ( \vert Du\vert^2 x_i\big)_{x_i} &=2 \sum_{i=1}^n \bigg ( \bigg ( Du\cdot x + \frac{n-2}{2} u\bigg) u_{x_i} \bigg)_{x_i} \end{align*} which can also be written as $$\mathrm{div} (\vert Du\vert^2x) = 2 \mathrm{div}\bigg ( \bigg ( \vert x \vert u_r + \frac{n-2}{2} u\bigg) Du \bigg ) $$ since $u_r = Du \cdot \frac x {\vert x \vert }$. Hence, the divergence theorem implies \begin{align} r \int_{\partial B(0,r)}\vert Du\vert^2 d S &= \int_{\partial B(0,r)} \vert Du\vert^2(x \cdot \nu )d S \\ &= \int_{B(0,r)} \mathrm{div} (\vert Du\vert^2x) dx \\ &= 2 \int_{B(0,r)}\mathrm{div}\bigg ( \bigg ( \vert x \vert u_r + \frac{n-2}{2} u\bigg) Du \bigg ) dx\\ &= 2 \int_{\partial B(0,r)} \vert x \vert u_r ( Du \cdot \nu ) d S+ (n-2)\int_{\partial B(0,r)} u( Du \cdot \nu )dS \\ &= 2r\int_{\partial B(0,r)} u_r^2dS + (n-2)\int_{\partial B(0,r)} u\frac{\partial u}{\partial \nu}dS \\ &= 2r\int_{\partial B(0,r)} u_r^2dS + (n-2)\int_{ B(0,r)} \vert Du \vert^2 dx \tag{2} \end{align} where the last line follow from the fact $u$ is harmonic. The above equality is again given in Evans. Substituting (2) into (1) gives \begin{align*} \frac d {dr} \bigg ( {r^{2-n}}\int_{B(0,r)} \vert Du \vert^2 d x\bigg ) &= \frac{1}{r^{n-1}} \bigg( 2r\int_{\partial B(0,r)} u_r^2dS - r \int_{\partial B(0,r)}\vert Du\vert^2 d S \bigg ) + \frac{1}{r^{n-2}} \int_{\partial B(0,r)} \vert Du \vert^2 d S \\ &= \frac 2 {r^{n-2}}\int_{\partial B(0,r)} u_r^2dS \geqslant 0 \end{align*} as required.