function [myval,Jacob] = evallpa(param)

global mu_l mu_a c_ea c_el c_pa b

% [myval,Jacob] = evallpa(param)
% Beetles, Cannibalism, and Chaos
% This function is used for Problem 2, to determine a fixed 
% point for the values L, P, and A.
% Dianne P. O'Leary 12/2006

% Get current value of L, P, and A from input vector

L = param(1);
P = param(2);
A = param(3);

% Compute values of L, P, and A at the next time step.

Lnew = b * A * exp(-c_ea*A - c_el*L);
Pnew = L * (1 - mu_l);
Anew = P * exp(-c_pa * A) + A * (1 - mu_a);

% Return the change in  L, P, and A.

myval = [Lnew-L; Pnew-P; Anew-A];

% Compute the Jacobian matrix of the function that
% defines new values - old.

ea = exp(-c_ea*A);
ep = exp(-c_pa*A);
el = exp(-c_el*L);
%              L           P             A
Jacob = [-b*A*c_el*ea*el,  0,       b*el*(ea-c_ea*A*ea)         % L
         1-mu_l,           0,       0                           % P
         0,                ep,      -c_pa*P*ep + (1-mu_a)] ...  % A
        - eye(3);