Solution (source code)

= Solution

Extend the boundary data oddly to the whole real axis. The extended function is precisely the function $f$ from part (a). Taking the <Fourier transform> in $x$, the <Laplace equation> becomes
$$
\partial_y^2\widetilde u-k^2\widetilde u=0.
$$
Decay as $y\to\infty$ selects
$$
\widetilde u(k,y)=\widetilde f(k)e^{-|k|y}.
$$
Part (b), together with the <convolution theorem>, says that the inverse transform is convolution with the <Poisson kernel>
$$
P_y(s)=\frac{y}{\pi(s^2+y^2)}.
$$
For $x,y>0$ this gives
$$
\begin{aligned}
u(x,y)
&=\frac y\pi\int_0^1\left[
\frac1{(x-v)^2+y^2}-\frac1{(x+v)^2+y^2}
\right]dv\\
&=\boxed{\frac{4xy}{\pi}\int_0^1
\frac{v\,dv}{[(x-v)^2+y^2][(x+v)^2+y^2]}}.
\end{aligned}
$$
The odd extension enforces $u(0,y)=0$, while the <Poisson integral> takes the prescribed boundary values at every continuity point and decays at infinity.