% Program 1.1
% Tumor model of eqs. (1.47), (1.48)
%
% Model parameters
  V0=1.0;
  lambda=1.0;
  alpha=1.0; 
%
% Step through cases
  for ncase=1:4
%
% Integration step
  if(ncase==1)h=1.0  ;nsteps=1   ;end
  if(ncase==2)h=0.1  ;nsteps=10  ;end
  if(ncase==3)h=0.01 ;nsteps=100 ;end
  if(ncase==4)h=0.001;nsteps=1000;end
%
% Variables for ODE integration
  tf=10.0;  
  t=0.0;
%
% Initial condition
  V1=V0;
  V2=V0;
%
% Print heading
  fprintf('\n\nh = %6.3f\n',h);
  fprintf(...
  '    t       Ve       V1          errV1          estV1       V2          errV2\n')
%
% Continue integration
  while t<0.999*tf
% 
% Take nsteps integration steps
  for i=1:nsteps
%
%   Store solution at base point
    V1b=V1;
    V2b=V2;
    tb=t;
%
%   RK constant k1
    k11=lambda*exp(-alpha*t)*V1*h;
    k12=lambda*exp(-alpha*t)*V2*h;
%
%   RK constant k2
    V1=V1b+k11;
    V2=V2b+k12;
     t=tb+h;
    k22=lambda*exp(-alpha*t)*V2*h;
%
%   RK step
    V2=V2b+(k12+k22)/2.0;
     t=tb+h;
  end
%
% Print solutions and errors
  Ve=V0*exp((lambda/alpha)*(1.0-exp(-alpha*t)));
  errV1=V1-Ve;
  errV2=V2-Ve;
  estV1=V2-V1;
  fprintf('%5.1f%9.4f%9.4f%15.10f%15.10f%9.4f%15.10f\n',...
           t,Ve,V1,errV1,estV1,V2,errV2);
%
% Continue integration
  end
%
% Next case
  end
