fun.m新增微分方程，通過RK遞推下一時刻的函式值

主函式如下：

n=90;

x=zeros(1,n+1);
y=zeros(1,n+1);
x(1)=0;
y(1) =1;   %初值
h=0.1;

for i =1:n
    
    x(i+1) = x(i) + h;
    
    y(i+1) = RK(@fun,x(i),y(i),h);
    
    
end

plot(x,y,'-o')

fun.m如下：y' = 2*(1-y/20)*y -x;  

function dy= fun(x,y)

dy  = 2*(1-y/20)*y -x;   

RK.m如下：

function y = RK(F_xy,x,y,h)


k_1 = F_xy(x,y);

k_2 = F_xy(x+0.5*h,y+0.5*h*k_1);

k_3 = F_xy((x+0.5*h),(y+0.5*h*k_2));

k_4 = F_xy((x+h),(y+k_3*h));


y = y + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;  

 結果如下：