%
% Program 1.6
% Tumor model of eqs. (1.47), (1.48)
% (or eqs. (1.3), (1.4), (1.5))
%
% Model parameters
  y0=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
  y5=y0;
%
% Print heading
  fprintf('\n\nh = %6.3f\n',h);
  fprintf(...
  '    t       ye       y4          erry4          esty4       y5          erry5\n')
%
% Continue integration
  while t<0.999*tf
% 
% Take nsteps integration steps
  for i=1:nsteps
%
%   Store solution at base point
    yb=y5;
    tb=t;
%
%   RK constant k1
    k1=lambda*exp(-alpha*t)*y5*h;
%
%   RK constant k2
    y5=yb+0.25*k1;
     t=tb+0.25*h;
    k2=lambda*exp(-alpha*t)*y5*h;
%
%   RK constant k3
    y5=yb+(3.0/32.0)*k1...
         +(9.0/32.0)*k2;
     t=tb+(3.0/8.0)*h;
    k3=lambda*exp(-alpha*t)*y5*h;
%
%   RK constant k4
    y5=yb+(1932.0/2197.0)*k1...
         -(7200.0/2197.0)*k2...
         +(7296.0/2197.0)*k3;
     t=tb+(12.0/13.0)*h;
    k4=lambda*exp(-alpha*t)*y5*h;
%
%   RK constant k5
    y5=yb+( 439.0/ 216.0)*k1...
         -(   8.0       )*k2...
         +(3680.0/ 513.0)*k3...
         -( 845.0/4104.0)*k4;
     t=tb+h;
    k5=lambda*exp(-alpha*t)*y5*h;
%
%   RK constant k6
    y5=yb-(   8.0/  27.0)*k1...
         +(   2.0       )*k2...
         -(3544.0/2565.0)*k3...
         +(1859.0/4104.0)*k4...
         -(  11.0/  40.0)*k5;
     t=tb+0.5*h;
    k6=lambda*exp(-alpha*t)*y5*h;       
%
%   RK step
    y4=yb+(    25.0/ 216.0)*k1...
         +(  1408.0/2565.0)*k3...
         +(  2197.0/4104.0)*k4...
         -(     1.0/   5.0)*k5;     
    y5=yb+(    16.0/ 135.0)*k1...
         +( 6656.0/12825.0)*k3...
         +(28561.0/56430.0)*k4...
         -(    9.0/   50.0)*k5...
         +(    2.0/   55.0)*k6;
    esty4=y5-y4;
    t=tb+h;
  end
%
% Print solutions and errors
  ye=y0*exp((lambda/alpha)*(1.0-exp(-alpha*t)));
  erry4=y4-ye;
  erry5=y5-ye;
  fprintf('%5.1f%9.4f%9.4f%15.10f%15.10f%9.4f%15.10f\n',...
           t,ye,y4,erry4,esty4,y5,erry5);
%
% Continue integration
  end
%
% Next case
  end



