Double integral $\int_0^\infty \int_0^\infty \frac{e^{-(x+y)/2}x^j y^k}{ax+by+c} dx dy$ In an ideal scenario I would like to have an analytic expression for the integral
$$
\operatorname{I}\left(\,{a,b,c,j,k}\,\right) =
\int_{0}^{\infty} \int_{0}^{\infty}
\dfrac{\operatorname{e}^{-\left(\,{x + y}\,\right)/2}\,\,x^{j}\, y^{k}}{ax + by + c}\,{\rm d}x\,{\rm d}y
$$
for positive real $a,b,c$ and non-negative integers $j,k$, which can then be evaluated for very large values of $c$. The reason I mention the latter circumstance is because while Mathematica ($13.1$ for sure) is able to take this integral for small values of $j$ and $k$, the resulting expression is numerically unstable for large values of $c$ (meaning that it's a product of terms exponentially large and small in $c$). For example, in the case $j = k = 0$ one gets:
$$
\frac{2 e^{\frac{c}{2 b}} \operatorname{Ei}\left(-\frac{c}{2 b}\right)-2 e^{\frac{c}{2 a}} \operatorname{Ei}\left(-\frac{c}{2 a}\right)}{a-b}
$$
Would appreciate any advice on how to:

*

*Take the original integral analytically,

*Efficiently evaluate $\operatorname{e}^\alpha \operatorname{Ei}(-\alpha)$ for large values of $\alpha$.

Given that Mathematica can take this integral for small $j$ and $k$, I would guess that there exists a recursive relation.
 A: $\newcommand{\bbx}[1]{\,\bbox[15px,border:1px groove navy]{\displaystyle{#1}}\,}
 \newcommand{\braces}[1]{\left\lbrace\,{#1}\,\right\rbrace}
 \newcommand{\bracks}[1]{\left\lbrack\,{#1}\,\right\rbrack}
 \newcommand{\dd}{\mathrm{d}}
 \newcommand{\ds}[1]{{\displaystyle #1}}
 \newcommand{\expo}[1]{\,\mathrm{e}^{#1}\,}
 \newcommand{\ic}{\mathrm{i}}
 \newcommand{\on}[1]{\operatorname{#1}}
 \newcommand{\pars}[1]{\left(\,{#1}\,\right)}
 \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}}
 \newcommand{\root}[2][]{\,\sqrt[#1]{\,{#2}\,}\,}
 \newcommand{\sr}[2]{\,\,\,\stackrel{{#1}}{{#2}}\,\,\,}
 \newcommand{\totald}[3][]{\frac{\mathrm{d}^{#1} #2}{\mathrm{d} #3^{#1}}}
 \newcommand{\verts}[1]{\left\vert\,{#1}\,\right\vert}$
\begin{align}
\on{I}\pars{a,b,c,j,k} & \equiv\color{#44f}
{\int_{0}^{\infty}\int_{0}^{\infty}
{\expo{-\pars{x + y}/2}\,\,\,x^{j}\,y^{k} \over ax + by + c}
\,\dd x\,\dd y}
\\[5mm] & \sr{\substack{\,\,\,\,x/c\ \mapsto\ x\\[0.5mm]
\atop y/c\ \mapsto\ y\\[0.5mm] \mbox{}}}{=}
c^{j + k + 1}\int_{0}^{\infty}\int_{0}^{\infty}
{\expo{-c\pars{x + y}/2}\,\,\,x^{j}\,y^{k} \over ax + by + 1}
\,\dd x\,\dd y
\\[5mm] & \sr{{\rm as}\ c\ \to\ \infty}{\sim}
c^{j + k + 1}\pars{\int_{0}^{\infty}\expo{-cx/2}\,x^{j}\,\dd x}
\pars{\int_{0}^{\infty}\expo{-cy/2}\,y^{k}\,\dd y}\tag{1}\label{1}
\\[5mm] & = c^{j + k + 1}\
\bracks{\pars{2 \over c}^{j + 1}\Gamma\pars{j + 1}}
\bracks{\pars{2 \over c}^{k + 1}\Gamma\pars{k + 1}}
\end{align}
(\ref{1}) involves the $\underline{Laplace's\ Method}$ because the main contribution to the integral, as $\ds{c \to \infty}$, occurs in a neighborhood of "$\ds{x = 0\ \mbox{and}\ y = 0}$".  Therefore,
\begin{align}
\on{I}\pars{a,b,c,j,k} & \equiv\color{#44f}
{\int_{0}^{\infty}\int_{0}^{\infty}
{\expo{-\pars{x + y}/2}\,\,\,x^{j}\,y^{k} \over ax + by + c}
\,\dd x\,\dd y}
\\[5mm] & \sr{{\rm as}\ c\ \to\ \infty}{\sim} \bbx{\color{#44f}{2^{j + k + 2}\,\,\,\Gamma\pars{j + 1}\Gamma\pars{k + 1}} \over c}
\end{align}
A: For the second question
$$e^{\alpha } \,\text{Ei}(-\alpha )=\sum _{n=1}^{\infty } (-1)^n \,\frac{(n-1)!}{\alpha^n}$$ You will not need many terms for a very good approximation.
You could also use the $[k,k+1]$ Padé approximant $P_k$. For example
$$P_2=-\frac{\alpha  (3 \alpha +13)}{3 \alpha ^3+16 \alpha ^2+10 \alpha   -4}$$ whose error is $\frac {44}{3 \alpha^6 }$.
For $\alpha=100$, this would give
$$e^{100 } \,\text{Ei}(-100)=-\frac{7825}{790249}=-0.009901942299$$ while the exact value is $-0.009901942287$.
Using $P_2$ the relative error is smaller than $0.01$% as soon as $\alpha \geq 8$ and smaller than $0.001$% as soon as $\alpha \geq 14$.
