function I=romberg(a,b)
h=b-a;
T(1)=h/2*(fun(a)+fun(b));
m=1;
while 1
h=h/2;
S(1)=1/2*T(1)+h*sumf(2^(m-1),a,h);
for j=1:m
S(j+1)=S(j)+(S(j)-T(j))/(4^j-1);
end
if abs(S(m+1)-T(m))<1e-6
break;
end
T=S;
m=m+1;
end
I=S(m+1);
end
function f=sumf(m,a,h)
for j=1:m
y(j)=fun(a+(2*j-1)*h);
end
f=sum(y);
end
function f=fun(x)
f=1/(1+x^2);
end