% Script recurrence3.m
% This backward recurrence relation is stable

f = @(x,n) x.^n.*exp(x-1);
E = 0; % E(30)

fprintf(1, '%5s %18s %18s %18s\n', 'n', 'E', 'E_exact', 'E_error');
for n = 30:-1:1
   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