clear all

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 
   
cH = 4; 
cT = 6; 
cV = 0.05; 
cE = 3;
cVT = 4; 
aET = .25;

sH = .05; 
sT = .002; 
sV = .12; 
sE = .0012;
sVT = 1; 
sET = .0012;

dH = .01;
dT = .01;  % dying of healthy cell
dV = 0.01; 
dE = 0.03; 
dT = .88; % coexistence.
dT = 1; % clearing of tumor cell



mH = .2;
mT = 5/3; 
mV = .12; 
mE = 12;

D = 1008;
P0 = .00035;

r = .7;

init=[0.0888,0.001,0,.001,.0003];
Tran = 0; %
%Evol = input('Evolution Time:' ); % Length of evolution after transient. 
Evol = 180; %
blah=odeset('RelTol', 1e-16, 'AbsTol', 1e-16, 'BDF', 'on');

N=5; % # of discretization.
std=10; % # of last days steady states.
dT0=0.45;dT1=.95; % initial and end bifurcation values.
z=[];
for k=1:(N+1)
    dT = dT0+(k-1)*(dT1-dT0)/N;
    [t,y]=ode15s('tumormodel',[0 Evol],init,blah);
    for n=1:std
    [aa,bb]=min(abs(t-169-n));
    z=[z; [dT,y(bb,:)]];
    end
end
plot(z(:,1),z(:,2),'r.',z(:,1),z(:,3),'k.')

%plot(t,y(:,1),'r-',t,y(:,2),'k-',t,y(:,3),'m-',t,y(:,4),'b-',t,y(:,5),'g-')


