I was able to specify the closed-form with Maple.
$$\int_{0}^{a} \frac{1}{\sqrt{(x-a)^2+(x-b)^2}} \ dx$$
equals to
$$\frac {\sqrt 2}2 \left(\ln \left( \left( \sqrt {2}+2\,{\operatorname{csgn}} \left( a-b \right) \right) \left( a-b \right) \right) - \ln \left(2\,\sqrt {{a}^{2}+{b}^{2}}-\sqrt {2}\left(a+b\right) \right) \right),
$$
where $\operatorname{csgn}$ is the complex signum function.
If we assume that $b>a$, then we could simplify it into the form
$$\frac {\sqrt 2}2 \left( \ln \left( b-a \right)+\ln \left( \sqrt {2}-1 \right) -\ln \left( \sqrt {2} \sqrt {{a}^{2}+{b}^{2}}-a-b \right) \right).$$
Thanks in large part to the commenters for getting me back on the right track. Posting the key points of my solution here for completeness.
$$E[XY]=\int_0^\infty\int_0^yxye^{-y}\,\mathrm dx\,\mathrm dy=3$$
$$E[X]=\int_0^\infty\int_0^yxe^{-y}\,\mathrm dx\,\mathrm dy=1$$
$$E[Y]=\int_0^\infty\int_0^yye^{-y}\,\mathrm dx\,\mathrm dy=2$$
$$E[X^2]=\int_0^\infty\int_0^yx^2e^{-y}\,\mathrm dx\,\mathrm dy=2$$
$$E[Y^2]=\int_0^\infty\int_0^yy^2e^{-y}\,\mathrm dx\,\mathrm dy=6$$
$$\boxed{\mathrm{Cov}[X,Y]=E[XY]-E[X]E[Y]=1}$$
$$\boxed{\mathrm{Corr}[X,Y]=\frac{E[XY]-E[X]E[Y]}{\sqrt{(E[X^2]-E[X]^2)(E[Y^2]-E[Y]^2)}}=\frac1{\sqrt2}}$$
$$\boxed{E[X\mid Y=y]=\int_0^y\frac xy\,\mathrm dx=\frac y2}$$
$$\boxed{E[Y\mid X=x]=\int_x^\infty y^2e^{x-y}\,\mathrm dy=x+1}$$
$$\begin{cases}
E[X]=E[E[X\mid Y]]=E\left[\frac Y2\right]=\frac{E[Y]}2\\[1ex]
E[Y]=E[E[Y\mid X]]=E[X+1]=E[X]+1
\end{cases}\implies\boxed{E[X]=1}$$
$$E[X^2\mid Y=y]=\int_0^y\frac{x^2}y\,\mathrm dx=\frac{y^2}3$$
$$\mathrm{Var}[X\mid Y]=E[X^2\mid Y]-E[X\mid Y]^2=\frac{Y^2}3-\left(\frac Y2\right)^2=\frac{Y^2}{12}$$
$$E[\mathrm{Var}[X\mid Y]]=E\left[\frac{Y^2}{12}\right]=\frac{E[Y^2]}{12}=\frac12$$
$$\mathrm{Var}[E[X\mid Y]]=\mathrm{Var}\left[\frac Y2\right]=\frac{\mathrm{Var}[Y]}4=\frac{E[Y^2]-E[Y]}4=\frac12$$
$$\boxed{\mathrm{Var}[X]=E[\mathrm{Var}[X\mid Y]]+\mathrm{Var}[E[X\mid Y]]=1}$$
Best Answer
The problem is that you're not telling WolframAlpha that $\lambda_1$ and $\lambda_2$ are positive. If you don't tell it this, it assumes they could be any real number (or even complex), so it doesn't know if the integral is convergent. To fix this, use the Assumptions option in Integrate like this:
And it will give the correct answer: $(a^2+b^2)/[ab(a+b)] = a^{-1} + b^{-1} - 2(a+b)^{-1}$.