We have that for a function $f$ defined on some open subset $U \subset \mathbb{C}$ then the following if true:
Suppose $u=\mathrm{Re}(f), v=\mathrm{Im}(f)$ and that all partial derivatives $u_x,u_y,v_x,v_y$ exists and are continuous on $U$. Suppose further that they satisfy the Cauchy-Riemann equations. Then $f$ is holomorphic on $U$.
The proof for this is readily available, though there is a subtlety that I can't understand.
We essentially want to compute $\lim_{h \rightarrow 0} \dfrac{f(z+h)-f(z)}{h}$ where $h=p+qi \in \mathbb{C}$.
We need a relationship like $u(x+p,y+q)-u(x,y)=pu_x(x,y)+qu_y(x,y)+o(|p|+|q|)$. Why is this relationship true?
One way to look at this is by interpreting $f$ as a function from $\mathbb{R}^2 \to \mathbb{R}^2$. This allows you to apply whatever you know about differentiability in $\mathbb{R}^2$ to $\mathbb{C}$. In other words, you interpret $f$ as function $$ f \,:\, (u,v) \to \left(\textrm{Re}(f(u+iv)), \textrm{Im}(f(u+iv)\right) =: (f_1(u,v),f_2(u,v)) \text{.} $$
If all the partial derivatives $\frac{\partial f_1}{\partial u}$, $\frac{\partial f_1}{\partial v}$, $\frac{\partial f_2}{\partial u}$, $\frac{\partial f_2}{\partial v}$ exist and are continuous, $f$ is differentiable in the $\mathbb{R}^2$-sense, i.e. for every $\mathbf{x}=(u,v)$ there's a linear function $df$ which approximates $f$ around $\mathbf{x}$. In other words, $$\begin{aligned} f(\mathbf{x} + \Delta) &= f(\mathbf{x}) + df_{\mathbf{x}}\cdot \Delta + R_{\mathbf{x}}(\Delta) \quad\text{where}\quad \lim_{\Delta\to 0} \frac{\left|\left| R_{\mathbf{x}}(\Delta) \right|\right|}{||\Delta||} = 0 \quad\text{and}\quad \\ df_{\mathbf{x}} &= \left(\begin{matrix} \frac{\partial f_1}{\partial u} &\frac{\partial f_1}{\partial v} \\ \frac{\partial f_2}{\partial u} &\frac{\partial f_2}{\partial v} \end{matrix}\right) \text{.} \end{aligned}$$
If the cauchy-riemann differential equations hold, $df_\mathbf{x}$ is then a scaled rotation, i.e. there are $A_\mathbf{x}$ and $B_\mathbf{x}$ such that $$ df_\mathbf{x} = \left(\begin{matrix} A_\mathbf{x}& B_\mathbf{x} \\ -B_\mathbf{x}& A_\mathbf{x} \end{matrix}\right) \text{.} $$ This form of $df_\mathbf{x}$ allows you to translate back into $\mathbb{C}$, because its exactly these matrices which correspond to complex multiplication. Simply set $$ f'(z) = f'(u+iv) = A_{(u,v)} + iB_{(u,v)} \text{.} $$ and translate the approximation back to $\mathbb{C}$, i.e. $$ f(z+\Delta) = f(z) + f'(z)\cdot\Delta + R_z(\Delta) \quad\text{where}\quad \lim_{\Delta\to 0} \frac{\left| R_{\mathbf{x}}(\Delta) \right|}{|\Delta|} = 0 \text{,} $$ which is exactly what you want to show.