function R=Rungkuta4(f,a,b,N,ya)
%  f是微分方程右端函數句柄
%  a,b自變量的取值范圍[a,b]的左右端點。
%  N為區間等分個數或迭代步數。
%  h為步長。
%  ya為初值。
%  R=[x' y']是自變量x和y所組成的矩陣。
h=(b-a)/N;
x=zeros(1,N+1);
y=zeros(1,N+1);
x=a:h:b;
y(1)=ya;
for j=1:N
k1=feval(f,x(j),y(j));
k2=feval(f,x(j)+h/2,y(j)+h/2*k1);
k3=feval(f,x(j)+h/2,y(j)+h/2*k2);
k4=feval(f,x(j)+h,y(j)+h*k3);
y(j+1)=y(j)+(h/6)*(k1+2*k2+2*k3+k4);
end
R=[x' y']