A variant of the exponential integral Consider the following integral (for $x,y\in \mathbb{R}_{>0}$)
$$E(x,y) =  \int_0^1 \frac{\mathrm{e}^{-x/s-ys}}{s}\,\mathrm{d}s,$$
which is a variant of the usual exponential integral $E_1(x)$ to which it reduces when $y=0$.
I am interested in the efficient evaluation of $E(x,y)$, either by numerical or analytical means. One example would be e.g. a series representation; to give an example, consider expansion around $|ys|\ll 1$, which yields a series representation in generalized exponential integrals $E_n(x)$
\begin{align}
E(x,y) &=  \int_0^1 \frac{\mathrm{e}^{-x/s}}{s}\sum_{n=0}^\infty \frac{(-ys)^n}{n!} \,\mathrm{d}s
= \sum_{n=0}^\infty \frac{(-1)^n y^n}{n!}\int_0^1 \mathrm{e}^{-x/s} s^{n-1}\,\mathrm{d}s\\
&= \sum_{n=0}^\infty \frac{(-1)^n y^n}{n!}E_{n+1}(x).
\end{align}
This series, however, is rather poorly convergent when $y>1$. Can we do better?  
Context: This integral arises in the Ewald summation technique for 1D-periodic systems embedded in a 2D coordinate system; in my case, $x$ and $y$ take values from $0$'ish to $10$'ish. Presently, I'm forced to evaluate the integral by numerical quadrature, which is painstakingly slow.

$\int_0^1\dfrac{e^{-\frac{x}{s}-ys}}{s}~ds$
$=\int_\infty^1se^{-xs-\frac{y}{s}}~d\left(\dfrac{1}{s}\right)$
$=\int_1^\infty\dfrac{e^{-xs-\frac{y}{s}}}{s}~ds$
$=K_0(x,y)$ (according to https://core.ac.uk/download/pdf/81935301.pdf)
