Unstable calculation of intergation by parts

Calculate $\lim_{n\to\infty} I_n$ , where

$I_n = \frac{1}{\mathrm{e}}\int_0^1 \mathrm{e}^x x^n \,\mathrm{d}x.$

We can show that

$I_0 = \frac{1}{\mathrm{e}}\int_0^1 \mathrm{e}^x \,\mathrm{d}x = 1-\frac{1}{e},$

and

$I_n = \frac{1}{\mathrm{e}}\int_0^1 \mathrm{e}^x x^n \,\mathrm{d}x = 1-n I_{n-1}.$

Evaluate for $n=1,\dots,25$ :

I = 1 - 1/exp(1);
for n=1:25
    fprintf('Iter = %2d, Approx = %g\n',n,I);
    I = 1 - n * I;
end

Iter =  1, Approx = 0.632121
Iter =  2, Approx = 0.367879
Iter =  3, Approx = 0.264241
Iter =  4, Approx = 0.207277
Iter =  5, Approx = 0.170893
Iter =  6, Approx = 0.145533
Iter =  7, Approx = 0.126802
Iter =  8, Approx = 0.112384
Iter =  9, Approx = 0.100932
Iter = 10, Approx = 0.0916123
Iter = 11, Approx = 0.0838771
Iter = 12, Approx = 0.0773522
Iter = 13, Approx = 0.0717732
Iter = 14, Approx = 0.0669478
Iter = 15, Approx = 0.0627311
Iter = 16, Approx = 0.0590338
Iter = 17, Approx = 0.0554593
Iter = 18, Approx = 0.0571919
Iter = 19, Approx = -0.0294537
Iter = 20, Approx = 1.55962
Iter = 21, Approx = -30.1924
Iter = 22, Approx = 635.04
Iter = 23, Approx = -13969.9
Iter = 24, Approx = 321308
Iter = 25, Approx = -7.7114e+06

Which is clearly numerically unstable.

For $x\in(0,1)$ we can show that $\mathrm{e}^x<(e-1)x+1$ :

x = linspace(0,1,101);
plot(x, exp(x), 'DisplayName', 'e^x');
hold on;
plot(x, (exp(1)-1)*x+1, 'DisplayName', '(e-1)x+1');
legend('Location', 'NorthWest');

Then,

$0 < I_n < \frac{1}{\mathrm{e}}\int_0^1 (\mathrm{e}-1) x^{n+1}+x^n \,\mathrm{d}x = \frac{1}{\mathrm{e}}\left(\frac{\mathrm{e}-1}{n+2}+\frac{1}{n+1}\right) \mathop{\longrightarrow}^{n\to\infty} 0.$