function y=intf(x);
y=1./(1+3*x.^2);
function[s,n]=romberg1(a,b,eps)
if nargin<3
    eps=1e-6;
end
s=10;
s0=0;
k=2;
% intf(b)
t(1,1)=(b-a)*(intf(a)+intf(b))/2;
while(abs(s-s0)>eps)
    h=(b-a)/2^(k-1);
    w=0;
    if (h~=0)
        for i=1:(2^(k-1)-1)
            w=w+intf(a+i*h);
        end
        t(k,1)=h*(intf(a)/2+w+intf(b)/2);
        for p=2:k
            for i=1:(k-p+1)
                t(i,p)=(4^(p-1)*t(i+1,p-1)-t(i,p-1))/(4^(p-1)-1);
            end
        end
            s=t(1,k);
            s0=(t(1,k-1));
            k=k+1;
            n=k;
     else
         s=s0;
         n=-k;
     end    
end