I came across a problem on this site (which I now cannot find) asking about the distribution of sums of i.i.d. inverse gamma random variables for paramters $\alpha=\beta=1$. Defining $Z=X+Y$ with $X,Y\sim\operatorname{InvGamma}(1,1)$ we can apply the usual formula to obtain the density $$ f_Z(z)=\int_0^z(z-x)^{-2}\exp\left(-\tfrac{1}{z-x}\right)x^{-2}\exp\left(-\tfrac{1}{x}\right)\,\mathrm dx,\qquad z>0. $$ This integral seems deceptively hard which is mainly due to the exponential terms. Now we do have the useful fact that $$ \int(z-x)^{-2}\exp\left(-\tfrac{1}{z-x}\right)\,\mathrm dx=\exp\left(-\tfrac{1}{z-x}\right), $$ however, I don't think integrating by parts gets us very far. Another observation is that the integrand is symmetric about $x=z/2$ so that $$ f_Z(z)=2\int_0^{z/2}(z-x)^{-2}\exp\left(-\tfrac{1}{z-x}\right)x^{-2}\exp\left(-\tfrac{1}{x}\right)\,\mathrm dx. $$ But again I'm not sure this accomplishes much. How can we deal with the exponential terms to put this integral in a more useful form for evaluating? Taylor series expansions are certainly out. Maybe a clever substitution? The integrand is quite interesting and reminds me of a bump function. Here it is for $z=1$:
Substitute $y=z^2/[4x(z-x)]$. Then $0<x<z/2$ maps to $1<y<\infty$, and $$f_Z(z)=\frac{8}{z^3}\int_1^\infty\left(1-\frac1y\right)^{-1/2}\exp\left(-\frac{4y}{z}\right)\,{\rm d}y=-\frac{4}{z^3}g'\left(\frac2z\right),$$ where $$g(a)=\int_1^\infty\frac{e^{-2ay}\,{\rm d}y}{\sqrt{y(y-1)}}\ \underset{y=(1+t)/2}{\phantom{\Big[}=\phantom{\Big]}}\ \int_1^\infty\frac{e^{-a(1+t)}}{\sqrt{t^2-1}}\,{\rm d}t=e^{-a}K_0(a)$$ uses the modified Bessel function $K_0$. Hence $f_Z(z)=(4/z^3)e^{-2/z}\big(K_0(2/z)+K_1(2/z)\big)$.