I am working on an exercise, the teacher has given the answer to the exercise, but I can’t understand it.
Exercise Let $U$ be an open set of $\mathbb{C}$, $u: U\to \mathbb{R}$ be a harmonic function, $z_0\in U$. On note $x_0 = \operatorname{Re} z_0, y_0 = \operatorname{Im} z_0$ and $r_0 = d(z_0, \mathbb{C}\setminus U)$. Prove that for all $z \in D(z_0, r_0)$ (open disk of centered at $z_0$ with radius $r_0$) , $u(z) = \sum_{p, q\in\mathbb{N}} a_{p, q} (x-x_0)^{p}(y-y_0)^{q}$ where $a_{p, q} = \frac{\partial^{p+q} ~~u}{p!q!\partial x^{p}~\partial y^{q}}(z_0)$.
What I understand : $D(z_0, r_0)$ is a connected and simply connected open set of $\mathbb{C}$, so there exist a holomorphic function $f : D(z_0, r_0) \to \mathbb{C}$ such that $u_{|D(z_0, r_0)} = \operatorname{Re} f$. By Cauchy's integral formula, we have for all $z\in D(z_0, r_0)$
$$f(z) = \sum_{n=0}^{+\infty} \frac{f^{(n)}(z_0)}{n!}(z-z_0)^{n}. \quad (1)$$
For all $(x, y)\in\mathbb{R}^{2}$ such that $|x-x_0|+|y-y_0| < r_0$, by Fubini's theorem
$$(1) = \sum_{k = 0}^{+\infty}\sum_{n = k}^{+\infty} \mathrm{i}^{k} \binom{n}{k} \frac{f^{(n)}(z_0)}{n!} (x-x_0)^{n-k} (y-y_0)^{k} = \sum_{q = 0}^{+\infty}\sum_{p = 0}^{+\infty} \mathrm{i}^{q} \frac{f^{(n)}(z_0)}{p!q!} (x-x_0)^{p} (y-y_0)^{q}$$
Denote $c_{p, q} = \mathrm{i}^{q} \frac{f^{(n)}(z_0)}{p!q!}$ for every $p, q \in\mathbb{N}$. Then $(c_{p, q} + \overline{c_{p, q}}) /2 = a_{p,q}$ Until now, we prove the exercise for $z = x+\mathrm{i}y \in D(z_0, r_0)$ such that $|x-x_0|+|y-y_0| < r_0$.
What I don't understand : Then the anwser use Abel's theorem for two variables to conclude.
So, can someone help me solve this exercise (any alternative solution is welcome) and tell me the statement of Abel's theorem for two variables and how to conclude by Abel's theorem？
W.L.O.G, let $x_0 = 0$ and $y_0 = 0$. Denote $T := \{x+iy \mid |x| + |y| < r\}$ For $z = x + iy \in D(0, r_0)$ not in $T$, I choose $\theta \in \mathbb{R}$ such that $\mathrm{e}^{\mathrm{i}\theta} z \in T$. Then I do the same calcul above for $u_{\theta} : z\in \mathrm{e}^{\mathrm{i}\theta} U \mapsto u(\mathrm{e}^{\mathrm{-i}\theta}z)$, I find
$$u(z) = \sum_{p, q\in\mathbb{N}} a_{p, q} (x\cos\theta - y\sin\theta)^{p}(y\cos\theta + x\sin\theta)^{q-1}y \quad (2)$$