Romberg积分法(MATLAB)

% fprintf('请输入区间最小值:\n')
% a=input('');
% fprintf('请输入区间最大值:\n')
% b=input('');
% fprintf('请输入N:\n')
% N = input('');
% fprintf('请输入精度:\n')
% ac = input('');
a=0;
b=1;
N=25;
ac=0.000001;
T = [];
S = [];
C = [];
R = [];
T1=0;
T2=0;
S1=0;
S2=0;
C1=0;
C2=0;
R1=0;
R2=0;
tol = 0;
t = 0;
m=1;
n=1;
h =(b-a)/n;
T1 = h*(fun_rom(b)-fun_rom(a))/2;
T(1)=T1;
n=2;
h =(b-a)/n;
% i = 2;
for i=2:N
% while tol<ac
    ii=2^(i-1);
    sum=0;
    for k=1:ii
       sum = sum + fun_rom(a+(k-1/2)*h);
    end
    T2 = T1/2+h*sum/2;
    T(i) = T2;
    S2 = (4*T2-T1)/3;
    S(i)=S2;
    if m~=1
        C2 = (16*S2-S1)/15;
        C(i)=C2;
    end
    if m~=2
        R2 = (64*C2-C1)/63;
        R(i) = R2;
    end
    if m~=3
       tol = abs(R2-R1);
       if tol~=0
       if tol<ac
           t = t+1;
           break;
       end
       end
    end
        R1=R2;
        C1=C2;
        S1=S2;
        T1=T2;
        h=h/2;
        m=m+1;
%         if i>=N
%            break;
%         end
        i=i+1;
end

函数

function y=fun_rom(x)
% y = exp(x)*(x^2);
% y= exp(x)*sin(x);
% y = 4/(1+x^2);
% y = 1/(x+1);
% y=sin(x)/x;
y=cos(x)/((1-x)^2);
% y=cos(x)/(1-x)^2;
% y=(cos(x)/(1-x))^2;
end