function s = Romberg(a,b,eps)
%Romberg求積法進行數(shù)值積分，其中a與b為積分區(qū)間
%eps為允許的誤差值
if nargin == 2
    eps = 1.0e-6;
elseif nargin < 2
    error
    return
end
t1 = 10000;
t2 = 10000;
n = 2;
t(1,1) = 0.5*(b-a)*(ff1(a)+ff1(b));
while abs(t2-t1) >= eps
    area = 0.0;
    % n= n+1;
    h = (b-a)/2^(n-1);
    for i = 1:(2^(n-1))
        area = area + 0.5*h*(ff1(a+h*(i-1))+ff1(h*i+a));
    end
    t(n,1)=area;
    for j = 2:n
        for i = 1:(n-j+1)
            t(i,j) = (4^(j-1)*t(i,j-1)-t(i,j-1))/(4^(j-1));
        end
    end
    t1 = t(1,n);
    t2 = t(1,n-1);
    n=n+1;
end
s=t1;
return