%
%  Host-Parasitoids Model
%
%  Feb. 22, 2007, Math Biology, Jae-Hun Jung
%
%  H[n+1] = k H[n] exp(-a P[n])
%  P[n+1] = c H[n] [1 - exp(-a P[n])]
%

r = input(' r '); 
a = input(' a ');
c = input(' c ');
K = input(' K ');

nmax=400;
H0=50; P0=10;
subplot(1,2,1)
   plot(1,H0,'bo',1,P0,'rx')
   axis([0 nmax 0 nmax])
   hold on
subplot(1,2,2)
   plot(H0,P0,'go')
   axis([0 nmax 0 nmax])
   hold on

for n=2:nmax
    f = exp(-a*P0);
    kr = exp(r*(1-H0/K));
    H1 = kr*f*H0;
    P1 = c*H0*(1-f);
    subplot(1,2,1)
      plot([n-1;n],[H0;H1],'bo-',[n-1;n],[P0;P1],'rx-')
      axis([0 nmax 0 nmax])
       xlabel('n')
       ylabel('H_n, P_n')
    subplot(1,2,2)
      plot([H0;H1],[P0;P1],'go-')
      axis([0 nmax 0 nmax])
      xlabel('H_n')
      ylabel('P_n')
    H0=H1;
    P0=P1;
    pause
end