function [out1, out2, out3] = tumormodel(t,y,flag)

global cH cT cV cE cVT aET sH sT sV sE sVT sET dH dT dV dE ... 
       mH mT mV mE D P0 r dsgC 

H=y(1).*(y(1)>0);
T=y(2).*(y(2)>0);
V=y(3).*(y(3)>0);
E=y(4).*(y(4)>0);
P=y(5).*(y(5)>0);

BH=H2(P,cH,sH)*H;
DH=dH*H+mH*(H+T+V+E)*H;
%chemoH=(cH-tauC>0)*chemo(t,dosageC,frequencyC);

BT=H2(P,cT,sT)*T+H2(V,cVT,sVT)*T;
DT=dT*T+mT*(T+V)*T+H2(T,aET,sET)*E; % last term for antigenicity.
%chemoT=(cT/gammaT-tauC>0)*chemo(t,dosageC,frequencyC);

BV=H2(P,cV,sV)*T;
DV=dV*V;
%chemoV=(cV/gammaV-tauC>0)*chemo(t,dosageC,frequencyC);

BE=H2(P,cE,sE)*E*T;
DE=dE*E+mE*(H+V+E)*E;
%chemoE=(cE/gammaV-tauC>0)*chemo(t,dosageC,frequencyC);
%immunoE=immune(t,dosageI,frequencyI);

BP=r*(DH+DT+DV+DE); %+chemoE+chemoT+chemoV+chemoE;
DP=BH+BT+BV+BE;

if nargin < 3 | isempty(flag)
   out1 = [
      BH-DH; 
      BT-DT; 
      BV-DV;
      BE-DE;
      D*(P0-P)+BP-DP 
          ];
else
   switch(flag)
   case 'init'
     out1 = [0; 50];
     out2 = [.5,.5,.5,.5,.5];
     out3 = []; %odeset('Vectorized', 1);
     case 'events' 
     sec1 = P-P0/2;
     %sec2 = y(2)-.5;   
     out1 = [sec1];
     out2 = [0];
     out3 = [1];
   otherwise
      error(['Unknown flag ''' flag ''.']);
   end
end; 
