Solution (source code)

= Solution

Fourier transformation gives
$$
(k^2+m^2)\widehat G(k)=1,
\qquad
\widehat G(k)=\frac1{k^2+m^2}.
$$
For $x>0$, close the inverse-transform contour in the upper half-plane and take the residue at $k=im$; for $x<0$, close it in the lower half-plane. This gives the <one-dimensional modified Helmholtz Green function>
$$
\boxed{G(x)=\frac{e^{-m|x|}}{2m}}.
$$
The same formula works for complex $m$ with $\operatorname{Re}m>0$: it decays at both ends, is continuous at zero, and its <derivative> has the jump $G'(0+)-G'(0-)=-1$ required by the delta source.

Convolution therefore yields
$$
\boxed{
u(x)=\frac1{2m}\int_{-\infty}^{\infty}
 e^{-m|x-y|}f(y)\,dy}.
$$