T = 300;                                     %仿真時(shí)間  
t=1:T+1;
f=zeros(1,T);                                %幼年狐貍數(shù)
F=zeros(1,T);                                %成年狐貍數(shù)
r=zeros(1,T);                                %幼年嚙齒動(dòng)物數(shù)
R=zeros(1,T);                                %成年嚙齒動(dòng)物數(shù)
fox=zeros(1,T);                              %狐貍總數(shù)                  
rody=zeros(1,T);                             %嚙齒動(dòng)物總數(shù)
f(1)=6;
F(1)=10;
r(1)=5800;
R(1)=10000;                                  %設(shè)定初值
rody(1)=r(1)+R(1);
r1=4;                                        %狐貍出生率
r2=1.2;                                      %嚙齒動(dòng)物出生率
g=0.2;

for i=1:T
   
  temp = 10*f(i)+60*F(i)-g*R(i);             %假定每月被吃的嚙齒動(dòng)物為其g倍
  if(temp>0)                                 %有餓死狐貍的情況
     f(i+1)=f(i)-temp/(10*f(i)+60*F(i))*f(i);%餓死的小狐貍數(shù)
       
    if(i<=120)                                         %獵人出現(xiàn)之前
     F(i+1)=F(i)-temp/(10*f(i)+60*F(i))*F(i);
     elseif (F(i)-temp/(10*f(i)+60*F(i))*F(i)-4)>=0    %獵人出現(xiàn)后
       F(i+1)=F(i)-temp/(10*f(i)+60*F(i))*F(i)-4;
       else F(i+1)=0;
     end
       eat=g*R(i);                            %被吃的成年嚙齒動(dòng)物
   else                                       %沒有出現(xiàn)餓死的情況 
       f(i+1)=f(i);
       F(i+1) = F(i);
       eat=10*f(i)+60*F(i);
   end
   if(mod(i,8)==0)                            %幼年狐貍變成成年狐貍
       F(i+1)=F(i)+f(i);
       f(i+1)=r1*F(i);
      
   end
   if f(i+1)<=0
           f(i+1)=0;
   end
   
   R(i+1)=R(i)-eat;                           %成年嚙齒動(dòng)物
   r(i+1)=r(i)-eat/R(i)*r(i);                 %幼年嚙齒動(dòng)物
  
   if(mod(i,3)==0)                            %幼年長(zhǎng)成成年
       r(i+1)=(R(i)-eat)*r2;
       R(i+1)=(R(i)-eat)+r(i)-eat/R(i)*r(i);
   end
   if (r(i)+R(i)>10^5)                         %數(shù)量達(dá)到十萬(wàn)
       R(i)=10^5-r(i);
       R(i+1)=10^5-r(i);
   end
    
end


fox=f+F;
figure(1)
plot(t,fox,'r');
title('fox');xlabel('Time/Month'),ylabel('Numbers')

hold on;
figure(2)
rody=R+r;
plot(t,rody,'b:');
title('rody');xlabel('Time/Month'),ylabel('Numbers')

ylim([0,120000]);
fox,rody