function r=Romberg(a,b,n)
s=diag(zeros(n));      % n維列向量，初始化為0
h=b-a;
m=0;
s(1,1)=(feval('f',a)+feval('f',b))*h/2;
for k=2:n
    z=0;
    for i=1:2^(k-2)
        z=z+feval('f',a+(2*i-1)*h/2^(k-1));
    end
    s(k,1)=(s(k-1,1)+z*h/2^(k-2))/2;
end
r=zeros(n);
r(:,1)=s;    % 先求出r的第一列
for k=2:n
    for j=2:k
        r(k,j)=r(k,j-1)+(r(k,j-1)-r(k-1,j-1))/(4^(j-1)-1);
    end
end