function f=F3m(y,C,F)
f=C*y+F;
%我編寫的Lotka解微分方程組的程序
%M1*y1''+B1*y1'+k1*y1-B1*y2'-y1*y2=F1(t);
%-B1*y1'-k1*y1+M2*y2''+B1*y2'+(k1+k2)*y2-k2*y3=0;
%-k2*y2+M3*y3''+B3*y3'+(k2+k3)*y3=F3(t);
%令y4=y1';y5=y2';y6=y3';寫成矩陣形式y'=f(y,t)
M1=1;M2=1;M3=1;
k1=1;k2=1;k3=1;
F1=0.01;F3=0;
F=[0,0,0,F1/M1,0,F3/M3]';
B1=0.01;B3=0.01;
y(:,1)=[0;0;0;0;0;0];
t(1)=0;n=1;
h=0.01;
C=[0,0,0,1,0,0;...
0,0,0,0,1,0;...
0,0,0,0,0,1;...
-k1/M1,k2/M1,0,-B1/M1,B1/M1,0;...
k1/M2,-(k1+k2)/M2,k2/M2,B1/M2,-B1/M2,0;...
0,k2/M3,-(k2+k3)/M3,0,0,-B3/M3];
while t<=30;
K1=h*F3m(y(:,n),C,F);
K2=h*F3m(y(:,n)+K1/2,C,F);
K3=h*F3m(y(:,n)+K2/2,C,F);
K4=h*F3m(y(:,n)+K3,C,F);
y(:,n+1)=y(:,n)+(1/6)*(K1+2*K2+2*K3+K4);
t(n+1)=n*h;
n=n+1;
end
plot(t,y(1:3,);