I would really appreciate if you could direct me to a reference for the following fact.
Given a harmonic function $h$ defined in $R^N\backslash\{0\}$ we can find a holomorphic function $g$ of $N$ complex variables, such that $g$ coincides with $h$ at real points, and $g$ has a holomorphic extension to $C^N\backslash E$ where $E=\{z\in\mathbb{C}^N:z_1^2+z_2^2+\dots+z_N^2=0\}.$
Here holomorphic extention to $C^N\backslash E$ includes the case where there are two functions $g_1$ holomorphic on a domain $D_1$ and $g_2$ holomorphic on a domain $D_2,$ both coinciding with $h$ for real values, and $D_1\cup D_2 =C^N\backslash E.$
Thank you!
Let's try. First $N=2$ and suppose that $h$ is real-valued (it should work for complex valued too with a bit more work). First we complexify $h$ to some neighborhood of $\mathbb{R}^2 \setminus \{ 0 \}$ in $\mathbb{C}^2$. Locally near a point $h(z) = f(z) + \bar{f}(\bar{z})$ ($h$ is a real value of a holomorphic function), where $z = x+iy$. Starting at a point $(x_0,y_0)$ the function $f$ can be path continued to any $(x_1,y_1)$ in $\mathbb{R}^2$ outside the origin. Pick a $(\xi_0,\eta_0) \in \mathbb{C}^2$ outside of $\xi^2+\eta^2=0$. Near $(x_0,y_0)$ $h$ complexifies to $h(\xi + i \eta) = f(\xi+i\eta) + \bar{f}(\xi-i\eta)$, there exists a path $\xi(t),\eta(t)$ space from $(x_0,y_1)$ to $(\xi_0,\eta_0)$ which does not pass through $\xi^2+\eta^2 = 0$ (as the complement of this set is connected). Continue $f$ along the path $\xi(t)+i\eta(t)$ and $\bar{f}$ along the path $\xi(t)-i\eta(t)$ (both paths are in $\mathbb{C} \setminus \{0\} = \mathbb{R}^2 \setminus \{ 0 \}$). So you can continue the complexified $h$ along any path to any point as long as you stay away from $\xi^2+\eta^2=0$. That you cannot do this in a single-valued manner, the counterexample is $\log(x^2+y^2)$ as in the comment above.
For higher $N$, I assume you want to use Laurent series. As long as $N > 2$, then there exist harmonic homogeneous polynomials $p_m$ and $q_m$ such that $$h(x) = \sum_{m=0}^\infty p_m(x) + \sum_{m=0}^\infty \frac{q_m(x)}{|x|^{2m+N-2}}$$ See for example Axler-Bourdon-Ramey (http://www.axler.net/HFT.html), page 209. Convergence is uniform on compact subsets of $\mathbb{R}^N \setminus \{ 0 \}$ as usual. The proof would follow then similarly as above. The $p_m$ and the $q_m$ shouldn't present any problems. I haven't gone through the calculation, so I won't give details, but I'd assume it involves writing the denominator as $(x_1^2 + \cdots + x_N^2)^{(2m+N-2)/2}$, then find a path from a point in the real punctured plane to your desired point (as in the $N=2$ case) where we stay well away from (the complexified) $x_1^2 + \cdots + x_N^2 = 0$ and you prove that along this path the complexified series still converges. Comparison test using the fact that the original series converges uniformly and absolutely on any annulus should do the trick.