%% Unstable calculation of integration 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

%%
% 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.
% $$
