2020-11-12

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