% Script recurrence1.m

f = @(x,n) x.^n./(4*x+1); % integrand
y = log(5)/4; % first integral y0

fprintf('%5s %23s %23s %23s\n', 'n', 'y', 'yexact', 'yerror');
for n = 0:20
   yexact = quadl(@(x) f(x,n), 0, 1, 1E-15);
   yerror = yexact - y;
   fprintf('%5d %23.15e %23.15e %23.15e\n', n, y, yexact, yerror);
   y = 0.25*(1/(n+1) - y);
end

% numeric integration where integrand has a parameter,
% n in this case, can be done with anonymous functions

% f = @(x,n) x.^n./(4*x+1);
% quadl(@(x) f(x,20), 0, 1, 1E-15)