%4級runge-kutta法編寫成子函數
%函數的輸入變量為：T為x的取值上限，n為把[0，T]區(qū)間等份的數目，u0為方程的初始值
function u=rk(T,n,u0)
h=T/n;
u=[];
u(1)=u0;
for i=1:n-1
    w1=f(i*h,u(i));
    w2=f(i*h+0.5*h,u(i)+0.5*w1*h);
    w3=f(i*h+0.5*h,u(i)+0.5*w2*h);
    w4=f(i*h+0.5*h,u(i)+0.5*w3*h);
    u(i+1)=u(i)+(1/6)*h*(w1+2*w2+2*w3+w4);
end