About the Fourier-Legendre series of $f(x)=e^{-x}$ So for the function $f(x) = \exp(-x)$ I have the formula for the coefficients of
$$f(x) = \sum_{n=0}^{\infty}a_n P_n(x)$$
which is(by using Rodrigues formula) 
$$a_n = \frac{2n+1}{2} \int_{-1}^{1}\exp(-x)\frac{1}{2^nn!}\frac{d^n}{dx^n}(x^2-1)^ndx$$
I think my next step is the integrate by parts, 
$$a_n = \frac{2n+1}{2} -[\exp(-x)\frac{1}{2^nn!}\frac{d^{n-1}}{dx^{n-1}}(x^2-1)^n]^{1}_{-1} + \frac{2n+1}{2} \int_{-1}^{1}\exp(-x)\frac{1}{2^nn!}\frac{d^{n-1}}{dx^{n-1}}(x^2-1)^ndx = \frac{2n+1}{2} \int_{-1}^{1}\exp(-x)\frac{1}{2^nn!}\frac{d^{n-1}}{dx^{n-1}}(x^2-1)^ndx  $$
Since the term in the brackets is zero.
It seems to me that if we do this $n$ times (further $n-1$ times). We obtain.
$$a_n = \frac{2n+1}{2}\int_{-1}^{1}\exp(-x)\frac{1}{2^nn!}(x^2-1)^ndx$$
I think I have an idea on how to proceed.
EDIT:
Let
$$\alpha_n = \int_{-1}^{1}\exp(-x)(x^2-1)^ndx$$
Then  $$\alpha_n = 2n\int_{-1}^{1}\exp(-x)x(x^2-1)^{n-1}dx = $$ $$-2n(n-1)\alpha_{n-1} - 2n(n-1)\alpha_{n-2} - 2n(n-1)\int_{-1}^{1}\exp(-x)x(x^2-1)^{n-2}dx$$
I can not seem to rid the $x$ in the last integral.
I have integrated $xe^{-x}$ before integrating the second time then added and subtracted $\alpha_{n-2}$
 A: 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.
A: You have noted that
\begin{align}
   a_n & = \frac{2n+1}{2}\int_{-1}^{1}e^{-x}\frac{1}{2^n n!}\frac{d^{n}}{dx^{n}}(x^2-1)^{n}dx \\
    & = \frac{2n+1}{2^{n+1}n!}\int_{-1}^{1}e^{-x}\frac{d^{n}}{dx^{n}}(x^2-1)^ndx 
\end{align}
And,
\begin{align}
     \int fg^{(n)}dx&=fg^{(n-1)}-f^{(1)}f^{(n-2)}+f^{(2)}{g^{(n-3)}}-\cdots+(-1)^{n}\int f^{(n)}gdx
\end{align}
Because $(x^2-1)^n=(x-1)^n(x+1)^n$, then all derivatives of $(x^2-1)$ of order $k < n$ vanish at the endpoints of $[-1,1]$. So, using a standard binomial expansion of $(x^2-1)^n$,
\begin{align}
      a_n & = \frac{2n+1}{2^{n+1}n!}\int_{-1}^{1}e^{-x}(x^2-1)^ndx \\
  & = \frac{2n+1}{2^{n+1}n!}\int_{-1}^{1}e^{-x}\sum_{k=0}^{n}{{n}\choose{k}}(-1)^{k}x^{2k}dx \\
  & = \frac{2n+1}{2^{n+1}n!}\sum_{k=0}^{n}{{n}\choose{k}}\int_{-1}^{1}e^{-x}x^{2k}dx
\end{align}
Then
\begin{align}
    \int e^{-x}x^{2k}dx & = -e^{-x}x^{2k}+\int e^{-x}2kx^{2k-1}dx \\
      & = -e^{-x}x^{2k}-e^{-x}(2k)x^{2k-1}+\int e^{-x}(2k)(2k-1)x^{2k-2}dx \\
      & = -e^{-x}\{ x^{2k}+(2k)x^{2k-1}+(2k)(2k-1)x^{2k-2}+\cdots+(2k)! \}
\end{align}
It's ugly, but you do get the coefficients once you evaluate the above from $-1$ to $1$.
