%四階龍格庫塔法
x(1,1)=620;
x(1,2)=10;
x(1,3)=70;
x1=x(1,1);
x2=x(1,2);
x3=x(1,3);
t=0;
Tf=input('請輸入仿真時間Tf= ');
h=input('請輸入步長h= ');
for i=2:Tf/h+1 
    k(1,1)=-0.001*x(i-1,1)*x(i-1,2);
    k(1,2)=-0.001*[x(i-1,1)+h/2*k(1,1)]*[x(i-1,2)+h/2*k(1,1)];
    k(1,3)=-0.001*[x(i-1,1)+h/2*k(1,2)]*[x(i-1,2)+h/2*k(1,2)];
    k(1,4)=-0.001*[x(i-1,1)+h*k(1,3)]*[x(i-1,2)+h*k(1,3)];
    k(2,1)=0.001*x(i-1,1)*x(i-1,2)-0.072*x(i-1,2);
    k(2,2)=0.001*[x(i-1,1)+h/2*k(2,1)]*[x(i-1,2)+h/2*k(2,1)]-0.072*[x(i-1,2)+h/2*k(2,1)];
    k(2,3)=0.001*[x(i-1,1)+h/2*k(2,2)]*[x(i-1,2)+h/2*k(2,2)]-0.072*[x(i-1,2)+h/2*k(2,2)];
    k(2,4)=0.001*[x(i-1,1)+h*k(2,3)]*[x(i-1,2)+h*k(2,3)]-0.072*[x(i-1,2)+h*k(2,3)];
    k(3,1)=0.072*x(i-1,2);
    k(3,2)=0.072*[x(i-1,2)+h/2*k(3,1)];
    k(3,3)=0.072*[x(i-1,2)+h/2*k(3,2)];
    k(3,4)=0.072*[x(i-1,2)+h*k(3,3)];
    x(i,1)=x(i-1,1)+h*[k(1,1)+2*k(1,2)+2*k(1,3)+k(1,4)]/6;
    x(i,2)=x(i-1,2)+h*[k(2,1)+2*k(2,2)+2*k(2,3)+k(2,4)]/6;
    x(i,3)=x(i-1,3)+h*[k(3,1)+2*k(3,2)+2*k(3,3)+k(3,4)]/6;
    x1=[x1;x(i,1)];
    x2=[x2;x(i,2)];
    x3=[x3;x(i,3)];
    t=[t;t(i-1)+h];
end
plot(t,x1,'g');
hold on
plot(t,x2,'r');
hold on;
plot(t,x3,'b');
grid on;