Recall Poisson's integral formula for the half-plane says that if $U : \mathbb{R} \to \mathbb{R}$ is piecewise continuous and bounded, then
$$
f(x,y) = \frac{y}{\pi} \int_{-\infty}^{\infty} \frac{U(t)}{(x-t)^{2} +y^{2}} dt,
$$
is harmonic in the upper half-plane and has boundary values $f(x,0) = U(x)$.
So, given your boundary values, we have
$$
U(x) = \begin{cases} -1 & x < 0 \\ 1 & x > 0. \end{cases}
$$
Note, as we are integrating, we don't care what the value of $U(0)$ is.
Then,
\begin{equation}
f(x,y) = \frac{y}{\pi} \int_{-\infty}^{\infty} \frac{U(t) dt}{(x-t)^{2} +y^{2}} = \frac{y}{\pi} \left[ \int_{0}^{\infty} \frac{dt}{(x-t)^{2} + y^{2}} - \int_{-\infty}^{0} \frac{dt}{(x-t)^{2} + y^{2}}\right].
\end{equation}
I'll compute the first integral, and leave the second as practice. Using the change of variables $u = x-t$, we see that $d u = -dt$. So,
\begin{align*}
\int_{0}^{\infty} \frac{dt}{(x-t)^{2} + y^{2}} &= \int_{0}^{\infty} \frac{-du}{u^{2} + y^{2}}\\
& = \frac{-1}{y} \arctan\left( \frac{x-t}{y} \right) \bigg|_{t=0}^{\infty}\\
& = \frac{1}{y} \arctan\left(\frac{t-x}{y}\right) \bigg|_{t=0}^{\infty} \\
& = \frac{1}{y} \left( \frac{\pi}{2} - \arctan \frac{x}{y} \right).
\end{align*}
Now use similar technique to perform the other integral. Don't forget the $\frac{y}{\pi}$ out in front of Poisson's formula, and be wary of negative signs.
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.
Best Answer
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.