我运行时出错,找不到原因,请各位高手帮助我,非常感激!

function [val,M]=Romberg(f,ab,dalt)
% Romberg算法计 算积分
% [val,M]=Romberg(f,ab,dalt)
% val 返回积分值
% M T数表
% f 被积函数,函数文件名
% ab 积分区间
% dalt 精度
if nargin<3
dalt=1e-7; %设置默认精度1e-7
end
i=1; h=ab(2)-ab(1); t=0; j=0;
T(1,1)=(h/2)*(feval(f,ab)*ones(2,1));
while i<50 %最大跌?0次
i=i+1; h=h/2;
T=[T zeros(i-1,1);zeros(1,i)];
x=ab(1)+h:2*h:ab(2)-h;
y=feval(f,x)*ones(size(x')); % feval(f,x)求f在x点值
T(i,1)=T(i-1,1)/2+h*y; % 细化区间,求得梯形值
for t=2:i
j=t-1;
T(i,t)=(4^j*T(i,j)-T(i-1,j))/(4^j-1); %外推加速
end
if abs(T(i,i)-T(i-1,i-1))<=dalt %控制精度
break
end
end
if nargout==2
M=T; %按需求返回T数表
end
val=T(i,i); %积分值
脚本文件:
y=inline('2*exp(-x).*(x.^(-1/2))');
[I,M]=Romberg(y,[0 1],1.0e-5)
% Romberg算法计 算积分
% [val,M]=Romberg(f,ab,dalt)
% val 返回积分值
% M T数表
% f 被积函数,函数文件名
% ab 积分区间
% dalt 精度
if nargin<3
dalt=1e-7; %设置默认精度1e-7
end
i=1; h=ab(2)-ab(1); t=0; j=0;
T(1,1)=(h/2)*(feval(f,ab)*ones(2,1));
while i<50 %最大跌?0次
i=i+1; h=h/2;
T=[T zeros(i-1,1);zeros(1,i)];
x=ab(1)+h:2*h:ab(2)-h;
y=feval(f,x)*ones(size(x')); % feval(f,x)求f在x点值
T(i,1)=T(i-1,1)/2+h*y; % 细化区间,求得梯形值
for t=2:i
j=t-1;
T(i,t)=(4^j*T(i,j)-T(i-1,j))/(4^j-1); %外推加速
end
if abs(T(i,i)-T(i-1,i-1))<=dalt %控制精度
break
end
end
if nargout==2
M=T; %按需求返回T数表
end
val=T(i,i); %积分值
脚本文件:
y=inline('2*exp(-x).*(x.^(-1/2))');
[I,M]=Romberg(y,[0 1],1.0e-5)