function I=romberg(f,a,b,n);
I=zeros(n);
for i=1:n 
	I(i,1)=trapez(f,a,b,2^(i-1));
end;
for j=2:n 
  for i=j:n
    I(i,j)=(4^(j-1)*I(i,j-1)-I(i-1,j-1))/(4^(j-1)-1);
  end;
end
I=I(end,end);


function T=trapez(f,a,b,n)
h=(b-a)/n;
x=linspace(a,b,n+1);
y=f(x);
w=[1,2*ones(1,n-1),1];
T=h/2*dot(y,w);  
