To derive a recurrence relation for such coefficients is the same as finding the structure of the continued fraction of $e^2$ (explanation at the bottom of this answer): that can be done by exploiting the Gauss continued fraction for the hyperbolic cotangent, but it is easier to compute such coefficients directly.
If we start from
$$ g(t)=\int_{-1}^{1}\frac{e^{-x}}{\sqrt{1+t^2-2t x}}\,dx = \sum_{n\geq 0}\left(\int_{-1}^{1}P_n(x)e^{-x}\,dx \right)t^n\tag{1} $$
we may recover our coefficients through the Taylor series of $g(t)$ at $t=0^+$.
If we define $F(x)$ as the Dawson integral
$$ F(x)=e^{-x^2}\int_{0}^{x}e^{y^2}\,dy \tag{2}$$
we have:
$$ g(t) = \sqrt{\frac{2}{t}}\left[e\cdot F\left(\frac{1+t}{\sqrt{2t}}\right)-e^{-1}\cdot F\left(\frac{1-t}{\sqrt{2t}}\right)\right] \tag{3} $$
but the asymptotic expansion of $F(x)$ for large $x$ is given by:
$$ F(x) = \sum_{n\geq 0}\frac{(2n-1)!!}{2^{n+1} x^{2n+1}}\tag{4} $$
hence:
$$\begin{eqnarray*}[t^n]\,g(t) &=& [t^n]\left[e\cdot\sum_{m\geq 0}\frac{(2n-1)!! t^m}{(1+t)^{2m+1}}-e^{-1}\cdot\sum_{m\geq 0}\frac{(2n-1)!! t^m}{(1-t)^{2m+1}}\right]\end{eqnarray*}\tag{5}$$
but due to stars and bars we have:
$$ [t^n]\frac{t^m}{(1-t)^{2m+1}} = \binom{m+n}{2m},\qquad [t^n]\frac{t^m}{(1+t)^{2m+1}}=\binom{m+n}{2m}(-1)^{n+m} \tag{6}$$
as soon as $m\leq n$, hence by $(5)$ and $(6)$ we have:
$$ \int_{-1}^{1}P_n(x)\,e^{-x}\,dx =\color{red}{\sum_{m=0}^{n}(2n-1)!!\binom{n+m}{n-m}\left(e-\frac{(-1)^{m+n}}{e}\right)}\tag{7}$$
where the RHS is a linear combination of $e$ and $e^{-1}$ with integer coefficients. This has a very interesting consequence: since $e^{-x}$ is a $C^{\infty}$ function over $(-1,1)$, the coefficients
$$ c_n = \int_{-1}^{1}P_n(x)\,e^{-x}\,dx $$
decay very fast to zero as $n\to +\infty$ by the Riemann-Lebesgue/Paley-Wiener theorem, so $(7)$ provides extremely good rational approximation of $e^2$. For instance:
$$ \int_{-1}^{1}P_4(x)e^{-x}\,dx = 36e-\frac{266}{e}\quad\Longrightarrow\quad e^2\approx\color{red}{\frac{133}{18}}$$
where the approximation error is $\leq 1.673\cdot 10^{-4}$. With $n=5$ we get $e^2\approx\color{red}{\frac{2431}{329}}$ whose approximation error is $\leq 1.652\cdot 10^{-6}$. And these approximation are so good that they are convergents of the continued fraction of $e^2$ by Lagrange's theorem.
Now back to the original question about the recursion. Since
$$(x^2-1) P_n'(x) = nx P_n(x)-nP_{n-1}(x),\qquad P_{n+1}'(x)-P_{n-1}'(x)=(2n+1)P_n(x)\tag{8}$$
and
$$ c_n = \int_{-1}^{1}e^{-x}P_n(x) = \left. -e^{-x}P_n(x)\right|_{-1}^{1}+\int_{-1}^{1}e^{-x} P_n'(x)\,dx \tag{9} $$
we have:
$$ \color{red}{c_{n+1}-c_{n-1} = (2n+1)\,c_n} \tag{10} $$
since $P_n(1)=1$ and $P_n(-1)=(-1)^n$.
Now you may prove $(7)$ through $(10)+\text{induction}$, if you like it.
In the meanwhile, I hope you discovered many interesting facts.
