% foxes and rabbits
deltat=1;          	% time in days
t=0;			% starting time
r=zeros(20000,1); 	% array to hold number of rabbits as function of time
f=zeros(20000,1);	% array to hold number of foxes as function of time

betar=.01;		% see class notes for definitions of the
deltar=.005;		%     constants

betaf=.005;
deltaf=.01;

epsf=.0001;
epsr=.0001;

r(1)=100;		% initial rabbit population
f(1)=100;		% initial fox population

for i=2:20000 		% loop through the coupled equations 20,0000 times
  r(i)=r(i-1)+betar*r(i-1)*deltat-deltar*r(i-1)*deltat-epsf*f(i-1)*r(i-1)*deltat;
  f(i)=f(i-1)+betaf*f(i-1)*deltat-deltaf*f(i-1)*deltat+epsr*r(i-1)*f(i-1)*deltat;
end

plot(r,'b')		% plot the rabbit population in blue
xlabel('time (days)')
ylabel('population')
hold on
plot(f,'r')		% plot the fox population in red