我编写的Lotka解微分方程组的程序

				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				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,);