% HW2 code for beetle (lab visit #1)
% barnett 10/4/07

N=2000;
x = zeros(3,N+1);
b=7.48;
cel=0.012; cea=.009;cpa=.004; mul=.267; mup=0;
mua=.96;              % variable parameter

x(:,1) = [100; 100; 100 + 1e-8];
% starting in the quasi-periodic motion for mu_a=0.96: 
%x(:,1) = [13.4314879737091          149.610747643315          13.7574880057382];
for n=1:N
  x(1,n+1) = b*x(3,n) * exp(-cea*x(3,n) - cel*x(1,n));
  x(2,n+1) = x(1,n)*(1-mul);
  x(3,n+1) = x(2,n)*(1-mup)*exp(-cpa*x(3,n)) + x(3,n)*(1-mua);
end

l = 1000:N;
figure; plot3(x(1,l), x(2,l), x(3,l), '+'); axis vis3d;

figure; plot(x(1,:), '+-'); xlabel('n'); ylabel('L_n larvae');

