function R = romberg(func,a,b,n)
%  R = romberg(func,a,b,n) tries to evaluate the integral of the
%  function func from a to be using n level of Romberg integration, 
%  see the pseudocode at page 224.
% 
%  Example: R = romberg(@sin,0,pi,5)
h = (b-a);
R(1,1) = 0.5*h*(feval(func,a)+feval(func,b));
for i=1:n-1
    h = 0.5*h;
    sum = 0;
    for k=1:2:2^i-1
	sum = sum+feval(func,a+k*h);
    end
    R(i+1,1) = 0.5*R(i,1) + sum*h;
    for j=1:i
	R(i+1,j+1) = R(i+1,j)+(R(i+1,j)-R(i,j))/(4^j-1);
    end
end
