function D = richardson(phi,h,n)
%
% Performs Richardson extrapolation to approximate a 
% limit L using values of an approximating function phi. 
% It is assumed that 
%
%   L = phi(h) + a_2*h^2 + a_4*h^4 + ... 
%
% Inputs: 
%   phi  = function handle for the approximating function routine 
%   h    = initial approximation step 
%   n    = number of extrapolation levels 
%
% Output:
%   D    = vector of extrapolation values. The jth component of D
%          is the (j-1)st extrapolation value, i.e., D(1) = phi(h). 
%

np1 = n+1; 
D = zeros(np1,1); 

for j = 1:np1
  D(j,1) = feval(phi,(2^(j-1))\h); 
end

for j = 2:np1 
  D(j:np1,1) = D(j:np1,1) + (4^(j-1) - 1)\(D(j:np1,1) - D((j-1):(np1-1),1));
end


