Your approach is correct, I would use the mean value property to try to derive an analogue of the Cauchy estimates. Let $\Omega$ be the domain, and let $u$ be harmonic. See if you can show that if $B(z,r) \subset \Omega$, then
$$ |\partial_i u(z)| \le Cr^{-1} \|u \|_{L^{\infty}(\partial B(z,r))}$$
You will be using both the mean value property (since $\partial_i u$ is itself harmonic) and the divergence theorem. If you have further questions I can elaborate more.
Elaboration: The $L^{\infty}$, in the context of continuous functions (such as the harmonic functions in this problem), is just the norm corresponding to uniform convergence. So the sequence of harmonic functions $u_n$ converging uniformly on a compact set $K$ can be rewritten as $\|u - u_n\|_{L^{\infty}(K)} \rightarrow 0$ as $n \rightarrow \infty$.
To obtain the estimate in question, first use the mean value property to obtain:
$$ \partial_i u(z) = \frac{1}{\pi r^2} \int_{B(z,r)} \partial_i u(x,y) \, dx \, dy, $$
which follows from $\partial_i u$ itself being harmonic. By the divergence theorem, then this is equal to
$$ \frac{1}{\pi r^2} \int_{\partial B(z,r)} u \nu^i \, dS, $$
where $\nu^i$ is the $i$-th component of the unit vector normal to the surface $\partial B(z,r)$. Thus
$$ |\partial_i u(z)| \le \frac{1}{\pi r^2} \int_{\partial B(z,r)} |u| \, dS = \frac{2}{r} \|u\|_{L^{\infty}(\partial B(z,r))}$$
This allows you to control the pointwise convergence of the partial derivatives of $u_n$ by the uniform convergence of $u_n$. However, this in turn allows you to control the uniform convergence of the derivatives on a compact set, since then you can just use larger balls.
The general idea is to express $u$ in terms of $z$ and $\bar z$, guess which analytic function $f$ satisfies $2u=f+\bar f$ (which is $u=\operatorname{Re}f$), and verify the guess. The preliminary computations can be a bit sloppy (not paying attention to branching), since at the end the result is checked.
Examples:
- $u(x,y)=x^2-y^2 = \frac14 ((z+\bar z)^2+(z-\bar z)^2) = \frac12 z^2+\frac12 \bar z^2$. So we try $f(z) = z^2$.
- $u(x,y)=\ln(x^2+y^2) = \ln (z\bar z) = \ln z+\ln \bar z$. So we try $f(z) = 2\ln z$.
Best Answer
The derivative of a holomorphic function $f(x+iy) = u(x,y) + i v(x,y)$ is given by the partial derivatives of its real and imaginary parts: $f'(x+iy) = u_x(x,y) + i v_x(x,y)$. Now you can use the Cauchy-Riemann equations to get what you are after.