用matlab 龙贝格求积分

用matlab 龙贝格求积分
令f=1/(x+1),a=0,b=1,eps=10^-4
运行下面算法是总是出错,请问该怎么运行?
function[quad,R]=Romberg(f,a,b,eps)
h=b-a;
R(1,1)=h*(feval(f,a)+feval(f,b))/2;
M=1;J=0;err=1;
while err>eps
J=J+1;
h=h/2;
S=0;
for p=1:M
x=a+h*(2*p-1);
S=S+feval(f,x);
end
R(J+1,1)=R(J,1)/2+h*s;
M=2*M;
for k=1:J
R(J+1,k+1)=R(J+1,k)+(R(J+1,k)-R(J,k))/(4^k-1);
end
err=abs(R(J+1,J)-R(J+1,J+1));
end
quad=R(J+1,J+1);