Since $u(x,y) = \dfrac{x+y}{x^2+y^2}$,
$$\frac{\partial u}{\partial x} = \frac{1}{x^2+y^2}-\frac{2x(x+y)}{(x^2+y^2)^2} = \frac{x^2+y^2-2x^2-2xy}{(x^2+y^2)^2} = \frac{-x^2-2xy+y^2}{(x^2+y^2)^2}.$$
Likewise, by symmetry,
$$\frac{\partial u}{\partial y} = \frac{-y^2-2xy+x^2}{(x^2+y^2)^2}.$$
Putting these together shows that
$$\begin{align} f'(z) &= \frac{-x^2-2xy+y^2-i(-y^2-2xy+x^2)}{(x^2+y^2)^2} \\ &= \frac{-x^2-2xy+y^2+iy^2+2ixy-ix^2}{(x^2+y^2)^2}.\end{align}$$
The numerator can be rewritten a bit:
$$\begin{align}(-x^2+2ixy+y^2)+i(-x^2+2ixy+y^2) &= (1+i)(-x^2+2ixy+y^2) \\ &= -(1+i)(x^2-2ixy-y^2) \\ &= -(1+i)(x-iy)^2.\end{align}$$
Thus we have that
$$f'(z) = \frac{-(1+i)(x-iy)^2}{((x+iy)(x-iy))^2} = -\frac{1+i}{(x+iy)^2} = -\frac{1+i}{z^2}.$$
From here, it is a simple application of integration to get the final answer.
I'll take $A\subseteq\mathbb{C}$ to be an open set. A function $f:A\rightarrow
\mathbb{C}$ is holomorphic in $A$ if for every $z_{0}\in A$ there exists the
limit
$$
\lim_{z\rightarrow z_{0}}\frac{f(z)-f(z_{0})}{z-z_{0}}=\ell\in\mathbb{C},
$$
which just says that the derivative $f^{\prime}(z_{0})$ exists in $\mathbb{C}%
$. Some books unfortunately require $f^{\prime}$ to be continuous, which just
muddles the waters. It is not needed, you get it for free.
Once you know that $f$ is holomorphic, you can prove that $f^{\prime}$ is
continuous and it is itself holomorphic, and then you can prove that there
exist derivatives of any order and are all holomorphic. To prove this you fix
an open ball $B$ in $A$ and prove that for every rectifiable close curve
$\gamma$ contained in $B$ you have that $\int_{\gamma}f(z)\,dz=0$. Note that
to make sense of the integral you only need $f$ to be continuous and not even
differentiability. Cauchy proved this theorem assuming that $f^{\prime}$ was
continuous (and this is why some books add this to the definition) but later
Goursat proved that the result continues to hold assuming that $f$ is only differentiable.
Once you have $\int_{\gamma}f(z)\,dz=0$ for every rectifiable close curve
$\gamma$ contained in $B$ you prove Cauchy's formula, that is, that for every
$z_{0}\in B$
\begin{equation}
f(z_{0})=\frac{1}{2\pi i}\int_{\partial B(z_{0},r)}\frac{f(z)}{z-z_{0}%
}dz\label{cauchy formula}%
\end{equation}
where the closed ball $\overline{B(z_{0},r)}$ is contained in $B$. The point
now is that the right-hand side is a function $h(z_{0})$ which is given by an
integral depending on a parameter and as a function of $z_{0}$ you have that
$\frac{1}{z-z_{0}}$ is as regular as you want. So you can use theorems on
differentiation under the integral sign to conclude that that the right-hand
side has derivatives of any order and they are all continuous. In turn the
same is true for the left-hand side. So now you know that $f$ has derivatives
of any order and they are all continuous. At this point, if you define
$u(x,y):=$ real part of $f(x+iy)$ and $v(x,y):=$immaginary part of $f(x+iy)$,
then $u$ and $v$ are $C^{\infty}$ because they are given by the composition of
the function $(x,y)\mapsto x+iy$ which is $C^{\infty}$ with the function $f$
which is also $C^{\infty}$. In particular you can use the chain rule and take
second derivatives of $u$ and $v$ and prove that the Cauchy-Riemann equations
are satisfied and that $u$ and $v$ is harmonic (the mixed derivatives cancel
out because $u$ and $v$ are $C^{\infty}$). Hope this answers your second questions.
Now the converse implication is not perfect. There is a nice review article of
Gray and Morris
. One good result is the following. Take $U\subseteq
\mathbb{R}^{2}$ be an open set and let $u:U\rightarrow\mathbb{R}$ and
$v:U\rightarrow\mathbb{R}$ be such that there exist the partial derivatives
$\partial_{x}u$ and $\partial_{y}u$ and $\partial_{x}v$ and $\partial_{y}v$.
Note that the existence of partial derivatives DOES not imply that $u$ and $v$
are differentiable in $U$. You also assume that $u$ and $v$ satisfy the
Cauchy-Riemann equations.
These hypotheses are not enough, since they do not even imply that $u$ and $v$
are continuous. The function
$$
f(z):=\left\{
\begin{array}
[c]{ll}%
\exp(-z^{-4}) & \text{if }z\neq0\\
0 & \text{if }z=0,
\end{array}
\right.
$$
does not have a derivative at $z=0$ but the corresponding functions $u$ and
$v$ satisfy the Cauchy-Riemann equations in $\mathbb{R}^{2}$.
So you need an extra hypothesis. There are several variants. One that I like
is that $u:U\rightarrow\mathbb{R}$ and $v:U\rightarrow\mathbb{R}$ admit
$\partial_{x}u$ and $\partial_{y}u$ and $\partial_{x}v$ and $\partial_{y}v$ in
$U$, $u$ and $v$ satisfy the Cauchy-Riemann equations, and that the function
$f(z):=u(x,y)+iv(x,y)$, where $z=x+iy$, is continuous in the domain
$A=\{z=x+iy:\,(x,y)\in U\}$. The idea of the proof is to to mollify the
functions $u$ and $v$, taking $u_{\varepsilon}=\varphi_{\varepsilon}\ast u$
and $v_{\varepsilon}=\varphi_{\varepsilon}\ast v$, where $\varphi
_{\varepsilon}$ is a nice kernel, and prove that the $C^{\infty}$ functions
$u_{\varepsilon}$ and $v_{\varepsilon}$ satisfy the Cauchy-Riemann equations.
Using that, you can prove that $\int_{\gamma}f_{\varepsilon}(z)\,dz=0$ for
every rectifiable close curve $\gamma$ contained in $B$, where $f_{\varepsilon
}(z):=u_{\varepsilon}(x,y)+iv_{\varepsilon}(x,y)$, where $z=x+iy$, and then
continue as before to prove that $f_{\varepsilon}$ is analytic. In turn you
get the Cauchy formula (1) for $f_{\varepsilon}$ and then since $f$ is
continuous, you can let $\varepsilon\rightarrow0$ to prove that $f$ satisfies
the Cauchy formula and so is holomorphic.
Hope this answers your first question. Read the paper, it is well-written.
Best Answer
Writing $f(z)=\sum_{n=0}^\infty c_n z^n$, we see that $$u=\frac12(f+\bar f) = \frac12\sum_{n=0}^\infty (c_n z^n+\bar c_n\bar z^n)\tag1$$ The $n$th term in the series in (1) is a homogeneous polynomial in $x,y$ of degree $n$. Therefore, there is no cancellation of monomials $x^ky^m$ between the terms with different $n$.
Conclusion: for $d\ge 1$, $u$ is a harmonic polynomial of degree $d$ if and only if $f$ is a complex polynomial of degree $d$. (The case $d=0$ is a bit special, since zero constant and nonzero constants have different degrees.)
The coefficient of $x^n$ in (1) is $\operatorname{Re}c_n$ while the coefficient of $x^{n-1}y$ is $i\operatorname{Im}c_n$. This gives you $c_n$ for every $n\ge 1$. (Of course, for $n=0$ you don't know $\operatorname{Im}c_0$ for the lack of $x^{n-1}y$.) Having $c_n$, you can write down $$v=\frac12(f-\bar f) = \frac12\sum_{n=0}^\infty (c_n z^n-\bar c_n\bar z^n)\tag2$$ and expand it into $x$ and $y$ monomials if you wish.
(Personally, I never find much need to expand harmonic functions into $x^ky^m$ monomials; the expansion into powers of $z$ and of $\bar z$ is more efficient because it encodes the harmonicity of the function.)