function [out1, out2, out3] = PerCapitaRate4(t,X, flag)
%
% x abiotic elements/nutrients
% y plants/producers
% z animals/consumers
%

global el n k tk r w a h b d m c A H B D M C AA HH BB DD MM CC



X=X.*(X>0);
x=(X(1:el).*(X(1:el)>0)); % abiotic elements/nutrients
y=(X((el+1):(el+n)).*(X((el+1):(el+n))>0)); % plants/producers
z=(X((el+n+1):(el+n+k)).*(X((el+n+1):(el+n+k))>0)); % animals/consumers
p=(X((el+n+k+1):(el+n+k+tk)).*(X((el+n+k+1):(el+n+k+tk))>0)); % animals/consumers


v=a./(1+diag(x)*(h.*a));

U=A*diag(1./(1+y'*(H.*A)));

UU=AA*diag(1./(1+z'*(HH.*AA)));

if nargin < 3 | isempty(flag)
   out1 = [r - w.*(x>0).*x - (diag((x>0).*x)*v)*y
           ((min(diag((x>0).*x)*(b.*v),[],1))' - d - m.*y - c*y - U*z)
           ((B.*U)'*y - D - M.*z - C*z - UU*p)
           ((BB.*UU)'*z - DD - MM.*p - CC*p)];
else
   switch(flag)
   case 'init'
     out1 = [0; 50];
     out2 = ones(1,el+n+k+tk);
     out3 = []; %odeset('Vectorized', 1);
     case 'events' 
     sec1 = 1;   
     out1 = [sec1];
     out2 = [0];
     out3 = [1];
   otherwise
      error(['Unknown flag ''' flag ''.']);
   end
end; 

end
