If we consider $f(z)=f(x,y)$, then this problem is equivalent to the 1st boundary problem of Laplas' equation:
$$\begin{cases}\frac{\partial^2 f}{\partial x^2}+\frac{\partial^2 f}{\partial x^2}=0 \\
|f(x,y)|\le C \text{ for all } x \in \mathbb{R} \text{ and for } y>0\\
f(x,0)=\begin{cases}1,x>0\\
-1,x<0
\end{cases}
\end{cases}$$
It has very simple solution with the help of fundamental operator:
$$f(x,y)=\frac{y}{\pi}\int\limits_{-\infty}^{+\infty}\frac{f_{0}(t)}{\left(x-t\right)^2+y^2}dt= -\frac{y}{\pi}\int\limits_{-\infty}^{0}\frac{1}{\left(x-t\right)^2+y^2}dt+\frac{y}{\pi}\int\limits_{0}^{+\infty}\frac{1}{\left(x-t\right)^2+y^2}dt$$
If you want to represent this function as a complex value function, you can calculate this integrals. The result is
$$f(x,y)=\operatorname{csgn}(\bar{y})\frac{\pi}{y} $$
-- this is that stuff which the program has calculated. In complex it is:
$$ f(z)=\operatorname{csgn}(\overline{\operatorname{Im} z})\frac{\pi}{\operatorname{Im} z}$$
According to the uniqueness of the solution, this solution is unique.
We still need the assumption that $u=0$ on $\partial U^+ \cap \{x_n=0\}$; otherwise $v$ is not continuous.
One option is to use the mean value characterization of harmonic functions: a continuous function $v$ is harmonic iff for every point $x$ there is $r>0$ such that $u(x)$ is equal to the average of $u$ on the sphere $\{x':|x'-x|=\rho\}$, for every $0<\rho<r$.
If you have the above characterization, the problem is easy: for $x$ with $x_n=0$ the spherical average is zero by symmetry.
But you can also use the Poisson integral formula as suggested. Apply it on the ball $B_r\{x: |x|<r\}$ for some $r\in (0,1)$, using $v$ as the boundary values. The Poisson integral defines a harmonic function $w$ in $B_r$ which agrees with $v$ on $\partial B_r$. By symmetry, $w=0$ when $x_n=0$. Thus, on the domain $B_r\cap\{x_n>0\}$ the functions $w$ and $u$ have the same boundary values. Since they are both harmonic, the maximum principle implies $w\equiv u $ in $B_r\cap\{x_n>0\}$. Since both $w$ and $v$ are antisymmetric (odd) in $x_n$, it follows that $w\equiv v$ in $B_r$.
Best Answer
That's not true in $n=2$. For example, $u(x,y)=y$ is harmonic in $\mathbb{R}^2$, and positive for $y > 0$, even though $u(x,0) \le 0$. And it's not true in $\mathbb{R}^n$ for the same reason.
By assuming that $u$ is bounded, then you get the Poisson integral representation of $u$ from its boundary function, and that will give you what you want.