romberg.m

matlab程序

function y=romberg(a,b)
for n=1:50
   for i=1:n+1
     h=(b-a)/n;
     x(i)=h*(i-1)+a;
     f(i)=13*(x(i)-x(i)^2)*exp(-2*x(i));
     end;
   t(n,1)=h*(sum(f)-f(1)-f(n+1));
end;
for j=2:50
   for n=j:50
    t(n,j)=(4^(j-1)*t(n,j-1)-t(n-1,j-1))/(4^(j-1)-1);
   end;
end;
y=0;
for k=2:50
    if t(k,k)-t(k-1,k-1)<0.0001;
        y=t(k,k);
    end;
end;