%
% Program 1.4
% 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;
%
% Print heading
  fprintf('\n\nh = %6.3f\n',h);
  fprintf(...
  '    t       ye       y1          erry1          esty1       y2          erry2\n')
%
% Continue integration
  while t<0.999*tf
% 
% Take nsteps integration steps
  for i=1:nsteps
%
%   Store solution at base point
    yb=y2;
    tb=t;
%
%   RK constant k1
    k1=lambda*exp(-alpha*tb)*y2*h;
%
%   RK constant k2
    y2=yb+k1;
    t=tb+h;
    k2=lambda*exp(-alpha*t)*y2*h;
%
%   RK step
    y1=yb+k1;
    y2=yb+(k1+k2)/2.0;
    esty1=y2-y1;
  end
%
% Print solutions and errors
  ye=y0*exp((lambda/alpha)*(1.0-exp(-alpha*t)));
  erry1=y1-ye;
  erry2=y2-ye;
  fprintf('%5.1f%9.4f%9.4f%15.10f%15.10f%9.4f%15.10f\n',...
           t,ye,y1,erry1,esty1,y2,erry2);
%
% Continue integration
  end
%
% Next case
  end



