%
% Program 1.7
% Tumor model of eqs. (1.47), (1.48)
% (or eqs. (1.3), (1.4), (1.5))
%
% Global variables
  global lambda alpha ncall;
%
% Model parameters
  lambda=1.0;
  alpha=1.0; 
%
% Select method
  for mf=1:2
%
% Error tolerances
  reltol=1.0e-02;
  abstol=1.0e-02;
  for ncase=1:4
    reltol=1.0e-02*reltol;
    abstol=1.0e-02*abstol;
%
%   Initialize counter for derivative evaluations
    ncall=0;
%
%   Variables for ODE integration
    t0=0.0;
    tf=10.0;
    tout=[t0:1.0:tf]'; 
    nout=11;
%
%   Initial condition
    y0=1.0;
%
%   Call ODE integrator
    options=odeset('RelTol',reltol,'AbsTol',abstol);
    if(mf==1)[t,y]=ode23('ode1p7',tout,y0,options); end    
    if(mf==2)[t,y]=ode45('ode1p7',tout,y0,options); end 
%       
%   Display solution and error
    fprintf('\n\n  mf = %1d\n  case = %1d\n  reltol = %6.2e\n  abstol = %6.2e\n\n',...
            mf,ncase,reltol,abstol); 
    fprintf('    t       ye        y           erry\n');
    for i=1:nout
      ye(i)=y0*exp((lambda/alpha)*(1.0-exp(-alpha*t(i))));
      erry(i)=ye(i)-y(i);
      fprintf('%5.1f%9.4f%9.4f%15.10f\n',t(i),ye(i),y(i),erry(i));
    end
    fprintf('\n  ncall = %5d\n',ncall);
%
% Next case
  end
%
% Next method
  end
%
% Plot last solution
  plot(t,y);
  xlabel('t')
  ylabel('y(t)')
  title(' Program 1.7, dy/dt = \lambda*exp(-\alpha*t)*y)')
  print pro1p7.ps


