Gradient estimate and divergence theorem in an open ball (PDE)

I need to prove the estimate

$ |\partial_{x_i}u(x_0)| \leq \frac{N}{R}u(x_0)$

where $ u \in C(B(x_0,R))$ and $ u$ is nonnegative and harmonic in $ B(x_0,R)$ , and $ i=1, \ldots,N$

I found this proof (see answer), but it’s the case where $ u$ is continuous up to the boundary $ \partial B$ of the ball of radius $ R$ . This hypothesis is crucial when I apply divergence theorem.

My attempt

Since $ u$ is harmonic in the ball, then also $ \partial x_i u(x)$ is harmonic for every $ i$ , as showed in the linked post.

I know that I have to consider smaller balls, so I choose $ \varepsilon>0$ and set $ r=(1-\varepsilon)R$ . Of course, $ u$ is continuous (and also harmonic) also at the boundary of the ball of radius $ r$ .

So I write \begin{align} \frac{\partial u}{\partial x_i}(x_0) = \frac{1}{|(B(x_0,r))|}\int_{B(x_0,r)}\frac{\partial u}{\partial x_i}(x)\, dx = \frac{1}{|(B(x_0,r))|}\int_{\partial B(x_0,r)}u(x)v_i\, dS(x) \end{align}

where $ v_i= \frac{x-x_{0i}}{r}$ is the i-th component of the outer normal.

Then I bound the absolute value of the previous term and by using the fact that $ u$ is nonnegative I get

\begin{align} |\frac{\partial u}{\partial x_i}(x_0) | \leq \frac{1}{|(B(x_0,r))|}\int_{\partial B(x_0,r)}u(x)\, dS(x)=\frac{N}{r}\frac{1}{|\partial B(x_0,r)|} \int_{\partial B(x_0,r)}u(x)\, dS(x)=\frac{N}{r} u(x_0) \end{align}

where the last equality comes from the MVP. So, up to now, I have

\begin{align} |\partial x_i u(x_0)| \leq \frac{N}{(1-\varepsilon)R} u(x_0) \end{align}

and by letting $ \varepsilon \rightarrow 0^{+}$ I have the thesis

Could it be okay?