Why the real and imaginary parts of a complex analytic function are not independent?
It's really just a question of the definition of the derivative. If $z=x+yi,$ $f(z)=u(x,y)+iv(x,y)$ can be any pair of functions $u,v.$
But if $f$ is differentiable, then:
$$f'(z)=\lim_{h\to 0}\frac{f(z+h)-f(z)}{h}\tag{1}$$
then $h$ can approach $0$ in many different ways, since $h$ is complex.
For example, you can have $h\to 0$ on the real line. Then: $$f'(z)=\frac{\partial u}{\partial x} +i\frac{\partial v}{\partial x}$$
But if $h\to 0$ along the imaginary part, then:
$$\begin{align}f'(z)&=\frac{1}{i}\left(\frac{\partial u}{\partial y}+i\frac{\partial v}{\partial y}\right)\\ &=\frac{\partial v}{\partial y}-i\frac{\partial u}{\partial y} \end{align}$$
So for the limit to be independent of any path you take $h\to 0$ you must have at minimum that $$\frac{\partial u}{\partial x}=\frac{\partial v}{\partial y},\\\frac{\partial v}{\partial x}=-\frac{\partial u}{\partial y}\tag{2}$$
So for (1) to be true, we need $u,v$ to satisfy the differential equations in (2).
It turns out that $(2)$ is enough to ensure that $(1)$ converges to a single value, but that is not 100% obvious.
The equations in (2) are called the Cauchy-Riemann equations.
Another way of looking at it is, given a function $f:\mathbb R^2\to\mathbb R^2$ mapping $\begin{pmatrix}x\\y\end{pmatrix}\mapsto \begin{pmatrix}u(x,y)\\v(x,y)\end{pmatrix}$ there is a matrix derivative standard from multi-variable calculus:
$$Df\begin{pmatrix}x\\y\end{pmatrix}=\begin{pmatrix}\frac{\partial u}{\partial x}&\frac{\partial u}{\partial y}\\\frac{\partial v}{\partial x}&\frac{\partial v}{\partial y}\end{pmatrix}\tag{3}$$
For small vectors $$\mathbf h=\begin{pmatrix}h_1\\h_2\end{pmatrix}$$ you get $f\left(\begin{pmatrix}x\\y\end{pmatrix}+\mathbf h\right)\approx f\begin{pmatrix}x\\y\end{pmatrix}+Df\begin{pmatrix}x\\y\end{pmatrix}\mathbf h.$
In particular, $Df$ is in some sense the "best" matrix, $\mathbf A,$ for estimating $f(\mathbf v+\mathbf h)\approx f(\mathbf v)+\mathbf A\mathbf h.$
Now, these matrices are not complex numbers. But an interesting fact is that the set of matrices of the form:
$$\begin{pmatrix}a&-b\\b&a\end{pmatrix}\tag{4}$$
are a ring isomorphic to the ring of complex numbers. Specifically, the above matrix corresponds to $a+bi.$
We also have that:
$$\begin{pmatrix}a&-b\\b&a\end{pmatrix}\begin{pmatrix}x\\y\end{pmatrix}= \begin{pmatrix}ax-by\\bx+ay\end{pmatrix}$$
compare that with:
$$(a+bi)(x+yi)=(ax-by)+(ay+bx)i.$$
So these matrices (4) act on $(x,y)^T$ the same way that $a+bi$ acts on $x+yi$ by multiplication.
The Cauchy-Riemann equations (2) just mean that $Df\begin{pmatrix}x\\y\end{pmatrix}$ is an example of (4) - that is, when the Cauchy-Riemann equations are true for $u,v$ then the multi-variate derivative (3) can be thought of as a complex number.
So we see that when we satisfy Cauchy-Riemann, $Df\begin{pmatrix}x\\y\end{pmatrix}\cdot\mathbf h$ can be seen as multiplication of complex numbers, $f'(z)$ and $h=h_1+h_2i.$ Then you have:
$$f(z+h)\approx f(z)+f'(z)h.$$
where $f'(z)$ is not just the best estimating complex number for this approximation, but also $f'(z)$ is the best linear operation on $h$ for this estimation.
So complex analysis is taking the vector function and asking, $f$ "when does it make sense to think of the derivative of the $\mathbb R^2\to\mathbb R^2$ as a complex number?" That is exactly when Cauchy-Riemann is true.
In the general case $f:\mathbb R^2\to\mathbb R^2,$ we can't really take the second derivative and get an estimate $f(z+h)\approx f(z)+Df(z)\cdot h +\frac{1}{2}D^2f(z)\cdot h^2+\cdots.$ We can't get easy equivalents to power series approximations of $f.$
But when $Df$ satisfies Cauchy-Riemann, we can think if $Df$ as a complex-valued function.
So complex analysis is a subset of the real analysis of functions $\mathbb R^2\to\mathbb R^2$ such that the derivative matrix $Df$ can be thought of as a complex number. This set of functions turns out to have a lot of seemingly magical properties.
This complex differentiability turns out to be fairly strong property on the functions we study. The niceness of the Cauchy-Riemann equations gives up some truly lovely results.
The most sensible way, I believe, to understand this is that differentiation characterizes the micro-local behavior of a function at a point. In particular, one of several interpretations of the derivative of any function $f$ at some input point $x_0$ is that in a suitably tiny region around that point, $f$ "acts" like (up to some shifts to get things suitably centered) a multiplication by $f'(x_0)$.
When $f$ is a complex function, then complex differentiability means that it must act like multiplication by a complex number, namely the complex derivative $f'(z_0)$ for a now complex test point $z_0$. And multiplication by a general complex number scrambles together the real and imaginary parts of the number so multiplied.
Since the independence was nicely addressed in the accepted answer, I address the relationship between the Cauchy integral formula and the maximum principle first. (Below the line there's a very similar explanation for why real and imaginary parts are not independent)
The Cauchy integral formula intuitively states that the value at a point is the average of the values in a circle around a point, weighted somehow according to their distance (and angle). This is not easy to see directly from the Cauchy-Riemann equations (i.e., linearity of the differential). It is quite a deep theorem that $f$ is holomorphic if and only if $f$ satisfies the Cauchy integral formula (the proof usually involves this cycle: holomorphic $\Rightarrow$ integral formula $\Rightarrow$ power series $\Rightarrow$ holomorphic). The nicest proof of the Cauchy integral formula I have seen so far just uses homotopy invariance and the fact that integrals along contractible curves are zero, but this is a story for another question I guess. Let me just remark that the Cauchy integral formula is true more generally, e.g.,
- The Mean Value Equality for harmonic functions
- The Stokes Equation for differential forms on smooth manifolds with boundary
Both of these results may offer geometric insight when the time comes. Perhaps it is best for now to think about "functions satisfying the Cauchy integral formula" just as "some class of functions that has this mean value property", knowing in the back of your head that you will eventually understand that these functions are actually the same as "functions that have $\mathbb C$-linear differential".
Believing in the Cauchy integral formula has the advantage of making the maximum principle seem intuitive: If the function is the weighted average of itself at all the circles around it, it cannot get bigger inside these circles than it is on the circles. If you look closely, this is probably basically the argument in your proof of the maximum principle.
Recall that the differential of a (real-)differentiable function $f: \mathbb R^2 \rightarrow \mathbb R^2$ at a point $x \in \mathbb R^2$ is the $\mathbb R$-linear map $Df_x: \mathbb R^2 \rightarrow \mathbb R^2$ that approximates $f$ best at the point $x$. You might have seen this as $$\lim_{h \rightarrow 0}\frac{f(x + h) - f(x) - Df_x(h)}{\Vert h \Vert} = 0.$$ This could be visualised by saying that every (real-)differentiable function $f: \mathbb R^2 \rightarrow \mathbb R^2$, after zooming in enough, looks like a linear transformation. Maybe check out some visualizations by 3blue1brown if you do not have a concrete picture in mind.
Now consider a holomorphic function $f: \mathbb C \rightarrow \mathbb C$. In this case, the differential of $f$ at $z$ is the $\mathbb C$-linear map $Df_z: \mathbb C \rightarrow \mathbb C$ that approximates $f$ best at $z$. Emphasis lies on the fact that $Df_z$ is $\mathbb C$-linear this time, i.e., $Df_z$ is just multiplication by some complex number denoted $f'(z)$.
Now comes the important observation: Let $a \in \mathbb C$ be a complex number. Then the map $z \mapsto az$ is given by stretching and rotating, but not shearing. As a map $\mathbb R^2 \rightarrow \mathbb R^2$, it has the form $$\begin{pmatrix} x & -y \\ y & x \end{pmatrix},$$ where $a = x + iy$ if that helps you to visualize it (otherwise, please ignore the matrix).
Why is this important? This shows, that when zooming in, $f$ looks like stretching and rotating. But, as you may gather from the above matrix, from intuition or just by blind faith, if you know what stretching and rotating does to one vector, then you also know what it does to all other vectors. This is why real and imaginary parts are not independent.