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.