function R=rk4(f,a,b,f0,h)
%f為函數f(x,y);
%a,b為取值范圍的左右端點；
%f0為函數初值y(0)，h為步長,n為迭代步數；
n=(b-a)/h;
T=zeros(1,n+1);
Y=zeros(1,n+1);
T=a:h:b;
Y(1)=f0;
for j=1:n
    K1=h*feval('f',T(j),Y(j));
    K2=h*feval('f',T(j)+h/2,Y(j)+K1/2);
    K3=h*feval('f',T(j)+h/2,Y(j)+K2/2);
    K4=h*feval('f',T(j)+h,Y(j)+K3);
    Y(j+1)=Y(j)+(K1+2*K2+2*K3+K4)/6;
end
t=[0:.1:5];
y=1./exp(10*t);
figure(1)
plot(T,Y,'*',T,Y,'r-',t,y,'k-');
hold on
R=[T' Y']