function R = romberg(f,a,b,n)
%
% Performs Romberg extrapolation to approximate the 
% definite integral of a function f over an interval 
% [a,b]. 
% 
% Inputs: 
%   f     = function handle for the integrand
%   a, b  = endpoints of the interval of integration (a < b) 
%   n     = number of extrapolation levels 
%
% NOTE: f must accept a vector argument x and produce a vector 
%       of function values. 
% 
% Output: 
%   R     = vector of Romberg extrapolates: R(1,1) is the simple 
%           trapezoid rule approximation and, for 1 <= n <= N, 
%	    R(n+1,1) is the nth approximate integral determined by 
%	    Romberg extrapolation.
% 

if b<=a, error('\n Must have a < b.\n'); end
if n<=0, error('\n N must be positive.\n'); end 

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

h = b-a;
R(1,1) = 2\h*(feval(f,a) + feval(f,b)); 

for j = 2:np1 
  h_old = h;
  h = 2\h; 
  xvals = a+h:h_old:b-h;
  fvals = feval(f,xvals); 
  R(j,1) = 2\R(j-1,1) + h*sum(fvals);
end 


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

