Solution (source code)

= Solution

For integrable functions $u,v$, the <convolution theorem> in the chosen convention is
$$
\widehat{u*v}(s)=\widehat u(s)\widehat v(s),\qquad (u*v)(x)=\int_{\mathbb R}u(x-t)v(t)\,dt.
$$
Let $F(s,y)$ be the <Fourier transform> of $f(\cdot,y)$ in the horizontal coordinate. Transforming <Laplace's equation> gives
$$
F_{yy}-s^2F=0,\qquad F(s,0)=\frac2{1+s^2}.
$$
For $s\ne0$, exclude the growing exponential and obtain $F(s,y)=2e^{-|s|y}/(1+s^2)$. At $s=0$, continuity in frequency gives $F(0,y)=2$; a constant zero-frequency transform is compatible with spatial decay of $f$.

By the transform pair in part (a), the inverse transform of $e^{-|s|y}$ is the <Poisson kernel for the upper half-plane>, $P_y(x)=y/[\pi(x^2+y^2)]$. Thus
$$
\boxed{f(x,y)=(P_y*e^{-|\cdot|})(x)=\frac y\pi\int_{\mathbb R}\frac{e^{-|t|}}{(x-t)^2+y^2}\,dt\quad(y>0).}
$$
The <Poisson extension of exponential boundary data> is harmonic, approaches the continuous boundary data and decays at spatial infinity. For uniqueness, apply the <maximum principle for harmonic functions> to the difference of two solutions on increasingly large half-discs: it is zero on the real boundary and uniformly small on the distant semicircle.