function [i]=romberg(f,a,b,n,e)

 

% function [i]=romberg(f,a,b,m,e)

% integral funkcije f v mejah od a do b z m podintervalov

 

h =b-a;

x=linspace(a,b,2);

y=eval(vectorize(f));

I(1,1)=h*sum(y)/2;

r=3;

for(k=2:n)

      h =h/2;  

      x=linspace(a,b,r);

      y=eval(vectorize(f));

      I(k,1)=h*(sum(y)-y(1)/2-y(r)/2);

      r=2*r-1;

 

    for(j=2:k)

        I(k,j)=(4^(j-1)*I(k,j-1)-I(k-1,j-1))/(4^(j-1)-1);

    end;

    i=I(k,k);

    if(abs((I(k,k)-I(k-1,k-1))/I(k,k))<e)return; end;

end;

return;