T = 300;                                     %仿真時間  
t=1:T+1;
f=zeros(1,T);                                %幼年狐貍數
F=zeros(1,T);                                %成年狐貍數
r=zeros(1,T);                                %幼年嚙齒動物數
R=zeros(1,T);                                %成年嚙齒動物數
fox=zeros(1,T);                              %狐貍總數                  
rody=zeros(1,T);                             %嚙齒動物總數
f(1)=6;
F(1)=10;
r(1)=5800;
R(1)=10000;                                  %設定初值
rody(1)=r(1)+R(1);
r1=4;                                        %狐貍出生率
r2=1.2;                                      %嚙齒動物出生率
g=0.2;
for i=1:T 
  temp = 10*f(i)+60*F(i)-g*R(i);             %假定每月被吃的嚙齒動物為其g倍
  if(temp>0)                                 %有餓死狐貍的情況
     f(i+1)=f(i)-temp/(10*f(i)+60*F(i))*f(i);%餓死的小狐貍數   
    if(i<=120)                                         %獵人出現之前
     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    %獵人出現后
       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);                            %被吃的成年嚙齒動物
   else                                       %沒有出現餓死的情況 
       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;                           %成年嚙齒動物
   r(i+1)=r(i)-eat/R(i)*r(i);                 %幼年嚙齒動物
   if(mod(i,3)==0)                            %幼年長成成年
       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)                         %數量達到十萬
       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]);