
%
%

clear all

parameters;

cH=8;
sH=.05;
dH=.01;
mH=.2;

cH=4;


cT = 8^(2/3); 
cV = 0.05; 
cVT = 8; 

cT = 6; 
cV = 0.05; 
cVT = 4; 

sT = .002; % capacityT = .006 kg guessing % capacityVT = 3 kg % Kuang04
sV = .12; % capacityV = .3 kg guess
sVT = 1; % capacityT = 3 kg % Kuang04

dT = .01;  
dV = 0.01; 

dT = .88; 


c = [cH, cT, cV,cE,cVT,aET];
s = [sH, sT, sV,sE,sVT,sET];
d = [dH, dT, dV,dE];
m = [mH, mT, mV,mE];

P0 = .000035; % 1\% phosphorus Kuang04
P0 = .00035;

tauC =5;

dosageC=0;
frequencyC=[1;0;50];
dosageI=0;
frequencyI=[0;0;0];

dosage=[dosageC,dosageI];
frequency=[frequencyC frequencyI];



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

% 
% 
% Generate the return points 'z' in matrix 'b' for the Poincare map.
%
[t,HTVEP]=ode15s('tumormodel',[0 Evol],init,blah,gamma,q,Q,c,s,d,m,D,P0,tauC,dosage,frequency,r);

%init=HTVEP(end,:);
%figure(1)
%plot(t,HTVEP(:,1),'r-',t,HTVEP(:,2),'k-',t,HTVEP(:,3),'m-',t,HTVEP(:,4),'b-',t,HTVEP(:,5),'g-')

%figure(2)
%plot3(HTVEP(:,1),HTVEP(:,2),HTVEP(:,4))

%[t,HTVEP]=ode15s('test',[0 Evol],init,blah,gamma,q,Q,c,s,d,m,D,P0,tauC,dosage,frequency,r);

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

%figure(1)
%plot(t,HTVEP(:,2),'k-',t,HTVEP(:,4),'b-',t,HTVEP(:,5),'g-')
