Look at the second-order ODE on the real line
$$ -u'' + u = \delta_0 $$
in the distributional sense, where $\delta_0$ denotes the Dirac delta at $0.$ Taking the Fourier transform of both sides, we can easily find out that $$ u(x) = {1 \over 2} e^{-|x|}$$ is a solution.
I am just reading Lieb-Loss Analysis, where it is written (on page 166) that the Yukawa potential
$$ u(x) = \int_0^\infty {1 \over \sqrt{4\pi t}} e^{-{x^2 \over 4t} - t} \: dt $$
solves the above ODE. Indeed, taking Fourier transforms of both sides, we have
$$ \int_0^\infty {1 \over \sqrt{4\pi t}} e^{-{x^2 \over 4t} - t} \: dt = {1 \over 2} e^{-|x|},$$
or one might use some uniqueness argument to get the above equality. But how on earth did someone come up with such an integral formula for the solution?
One possible way to get there: consider the diffusive PDE $$v_t -v_{xx}+v = \delta_0$$ with initial condition $v(x,0)=0$. Define the steady state solution as the $t\to\infty$ limit: $$u(x) = \lim_{t \to \infty} v(x,t).$$ Formally, the steady state solution should solve the ODE $$-u_{xx} + u = \delta_0.$$ On the other hand if you apply Duhamel’s formula to the above PDE using the fundamental solution to the heat equation, you get $$v(x,t) = \int^t_0 \frac 1{\sqrt{4 \pi s}} e^{-x^2/4s-s} ds$$ whereupon taking the limit as $t\to\infty$ yields the formula for $u(x)$ that you have.