%
% Program 1.3
% 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
  y4=y0;
  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
    y4b=y4;
    y5b=y5;
    tb=t;
%
%   RK constant k1
    k14=lambda*exp(-alpha*t)*y4*h;
    k15=lambda*exp(-alpha*t)*y5*h;
%
%   RK constant k2
    y4=y4b+0.25*k14;
    y5=y5b+0.25*k15;
     t= tb+0.25*h;
    k24=lambda*exp(-alpha*t)*y4*h;
    k25=lambda*exp(-alpha*t)*y5*h;
%
%   RK constant k3
    y4=y4b+(3.0/32.0)*k14...
          +(9.0/32.0)*k24;
    y5=y5b+(3.0/32.0)*k15...
          +(9.0/32.0)*k25;
     t= tb+(3.0/8.0)*h;
    k34=lambda*exp(-alpha*t)*y4*h;
    k35=lambda*exp(-alpha*t)*y5*h;
%
%   RK constant k4
    y4=y4b+(1932.0/2197.0)*k14...
          -(7200.0/2197.0)*k24...
          +(7296.0/2197.0)*k34;
    y5=y5b+(1932.0/2197.0)*k15...
          -(7200.0/2197.0)*k25...
          +(7296.0/2197.0)*k35;
     t= tb+(12.0/13.0)*h;
    k44=lambda*exp(-alpha*t)*y4*h;
    k45=lambda*exp(-alpha*t)*y5*h;
%
%   RK constant k5
    y4=y4b+( 439.0/ 216.0)*k14...
          -(   8.0       )*k24...
          +(3680.0/ 513.0)*k34...
          -( 845.0/4104.0)*k44;
    y5=y5b+( 439.0/ 216.0)*k15...
          -(   8.0       )*k25...
          +(3680.0/ 513.0)*k35...
          -( 845.0/4104.0)*k45;
     t= tb+h;
    k54=lambda*exp(-alpha*t)*y4*h;
    k55=lambda*exp(-alpha*t)*y5*h;
%
%   RK constant k6
    y4=y4b-(   8.0/  27.0)*k14...
          +(   2.0       )*k24...
          -(3544.0/2565.0)*k34...
          +(1859.0/4104.0)*k44...
          -(  11.0/  40.0)*k54;
    y5=y5b-(   8.0/  27.0)*k15...
          +(   2.0       )*k25...
          -(3544.0/2565.0)*k35...
          +(1859.0/4104.0)*k45...
          -(  11.0/  40.0)*k55;
     t =tb+0.5*h;
    k65=lambda*exp(-alpha*t)*y5*h;       
%
%   RK step
    y4=y4b+(    25.0/ 216.0)*k14...
          +(  1408.0/2565.0)*k34...
          +(  2197.0/4104.0)*k44...
          -(     1.0/   5.0)*k54;     
    y5=y5b+(    16.0/ 135.0)*k15...
          +( 6656.0/12825.0)*k35...
          +(28561.0/56430.0)*k45...
          -(    9.0/   50.0)*k55...
          +(    2.0/   55.0)*k65;
     t =tb+h;
  end
%
% Print solutions and errors
  ye=y0*exp((lambda/alpha)*(1.0-exp(-alpha*t)));
  erry4=y4-ye;
  erry5=y5-ye;
  esty4=y5-y4;
  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



