function r=romberg(a,b,n)

r=zeros(n);
h=b-a;

r(1,1) = h*(fun(a)+fun(b))/2;

L=1;

for i=2:n;
    h=h/2;
    L=L+L;
    
    S=0;
    for k=1:2:L-1;
        S = S + fun(a+k*h);
    end
    
    M=1;
    r(i,1) = r(i-1,1)/2 + h*S;
    
    for j=2:i;
        M=4*M;
        r(i,j) = r(i,j-1) + (r(i,j-1)-r(i-1,j-1))/(M-1);
    end
    for j=i+1:n;
        r(i,j)=0;
    end
    
end

end
