function [T] = romberg(@f,a,b,n)

% This routine performs a Romberg Integration


T = zeros(n+1,n+1);

h = b-a;

T(1,1) = 0.5*h*(f(a)+f(b));

for k=2:(n+1)
    h=h/2; % successive halving of step size
    x_odd = (a+h):(2*h):(b-h);
    f_odd = feval(f,x_odd);
    T(k,1) = T(k-1,1)/2 + h*sum(f_odd);

    for l=2:k % construct current row from previous
        T(k,l) = ( 4^(l-1)*T(k,l-1)-T(k-1,l-1) )/(4^(l-1)-1);
    end
end



