%  Kurt E. Anderson
%
%  December 23, 2014
%
%  DESCRIPTION:	a multipatch model of a two trophic level plant-herbivore community with
%               inducible plant defenses; represents the OBE model in
%               Anderson, Underwood, and Inouye, in review.
%
%  IN:	vector t: tinit:tinterval:tfinal
%		vector y: initial values
%		vector Z: lagged values
%		vector p: parameter values
%
%  OUT:		matrix with t and solutions for plant induction (P) and herbivores (H)

function dy = phfn_obe(~, y, Z, p)	% function input

N = length(y)/2;                    % patch number

%  VARIABLES

	P(1:N)    = y(1:N);	    % initial level of plant inducible defenses for each patch

    H(1:N)    = y(N+1:2*N); % initial density of herbivores
  
    Hlag(1:N) = Z(N+1:2*N); % lagged density of herbivores
    
    
%  PARAMETERS

    b       = p(1);        % induction half-maximum
    theta   = p(2);        % strength of threshold
    D0      = p(3);        % baseline dispersal rate of herbivores
    chi     = p(4);        % sensitivity of herbivore dispersal rate to induction level
    m       = p(5);        % mortality
   
    
% DYNAMIC EQUATIONS

	%plant population dynamics
        
    dP(1:N) = (1+(b^theta)).*(Hlag(1:N).^theta)./((b.^theta) + (Hlag(1:N).^theta))-P(1:N);
      
    % herbivore functions
    
    recr(1:1:N) = m.*(1-H(1:1:N));     % herbivore recruitment rate
    
    diff(1)           = 0.5*(D0+chi*P(2)).*H(2) + 0.5*(D0+chi*P(N-1)).*H(N-1) - (D0+chi*P(1)).*H(1);
    diff(2:1:N-1)     = 0.5*(D0+chi*P(1:1:N-2)).*H(1:1:N-2) + 0.5*(D0+chi*P(3:1:N)).*H(3:1:N) - (D0+chi*P(2:1:N-1)).*H(2:1:N-1);
    diff(N)           = diff(1);
    
    
    % herbivore population dynamics
    
    dH(1:1:N) = recr(1:1:N)+diff(1:1:N);
    
    dy = [dP'; dH'];	% pass to DDE23
    
end