%
% Program 1.2
% 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
  y2=y0;
  y3=y0;
%
% Print heading
  fprintf('\n\nh = %6.3f\n',h);
  fprintf(...
  '    t       ye       y2          erry2          esty2       y3          erry3\n')
%
% Continue integration
  while t<0.999*tf
% 
% Take nsteps integration steps 
  for i=1:nsteps
%
%   Store solution at base point
    y2b=y2;
    y3b=y3;
    tb=t;
%
%   RK constant k1
    k12=lambda*exp(-alpha*t)*y2*h;
    k13=lambda*exp(-alpha*t)*y3*h;
%
%   RK constant k2
    y2=y2b+(2.0/3.0)*k12;
    y3=y3b+(2.0/3.0)*k13;
     t=tb +(2.0/3.0)*h;
    k22=lambda*exp(-alpha*t)*y2*h;
    k23=lambda*exp(-alpha*t)*y3*h;
%
%   RK integration K3
    y3=y3b+(2.0/3.0)*k23;
     t=tb +(2.0/3.0)*h;
    k33=lambda*exp(-alpha*t)*y3*h;
%
%   RK step
    y2=y2b+(1.0/4.0)*k12+(3.0/4.0)*k22;
    y3=y3b+(1.0/4.0)*k13+(3.0/8.0)*k23+(3.0/8.0)*k33;
    t=tb+h;
  end
%
% Print solutions and errors
  ye=y0*exp((lambda/alpha)*(1.0-exp(-alpha*t)));
  erry2=y2-ye;
  erry3=y3-ye;
  esty2=y3-y2;
  fprintf('%5.1f%9.4f%9.4f%15.10f%15.10f%9.4f%15.10f\n',...
           t,ye,y2,erry2,esty2,y3,erry3);
%
% Continue integration
  end
%
% Next case
  end



