% Script recurrence2.m
% An unstable recurrence relation

n_max = input('Enter number of iterations ');
f = @(x,n) x.^n.*exp(x-1);
E = 1 - exp(-1);

fprintf(1, '%5s %18s %18s %18s\n', 'n', 'E', 'E_exact', 'E_error');
for n = 0:n_max
   E_exact = quadl(@(x) f(x,n), 0, 1, 1E-15);
   E_error = E_exact - E;
   fprintf(1, '%5d %18.10e %18.10e %18.10e\n', n, E, E_exact, E_error);
   E = 1 - (n+1)*E;
end