function [tout,yout] = rk_4(odefile,t0,h,th,y0)
%RK_4 實現了 4 階 Runge-Kutta 定步長微分方程求解。
%
%  [tout,yout] = rk_4(odefile,t0,h,th,y0)
%
%  odefile 為描述微分方程的 M 函數， t0, h 和 th 分別為
%  初始時間、步長和終止時間， y0 為狀態變量初值向量，返回
%  的 tout 和 yout 分別為時間向量和狀態變量解矩陣。

%Designed by Prof D Xue (c) 2000
tout=[t0:h:th]'; yout=[];
for t=tout'
   k1=h*eval([odefile '(t,y0)']);
   k2=h*eval([odefile '(t+h/2,y0+0.5*k1)']);
   k3=h*eval([odefile '(t+h/2,y0+0.5*k2)']);
   k4=h*eval([odefile '(t+h,y0+k3)']);
   y0=y0+(k1+2*k2+2*k3+k4)/6;
   yout=[yout; y0'];
end