#
# 2-ODE patient model
#
# Delete previous workspaces
  rm(list=ls(all=TRUE))
#
# Access ODE integrator
  library("deSolve");
#
# Access functions for numerical solution
  setwd("f:/comp3/wiley_download/chap6");
  source("ode_1.R");
#
# Level of output
#
#   ip = 1 - graphical (plotted) solutions 
#            (y1(x,t), y2(x,t) only)
#
#   ip = 2 - numerical and graphical solutions
#
  ip=2;
#
# Parameters
    VP=5000; VE=(5/3)*VP; lam12=0.0047; lam21=0.0017;  QB=150;
    QE=0.01;     alpha=1;    y2in=0.01;        y10=1;   y20=1;
  lam21=0.017;
  cat(sprintf("\n  VP = %4.1f  VE = %4.1f\n",VP,VE));  
  cat(sprintf("\n  lam12 = %8.5f  lam21 = %8.5f\n",lam12,lam21));
  cat(sprintf("\n  QB = %4.1f  QE = %8.5f  alpha = %5.2f\n",
              QB,QE,alpha));
  cat(sprintf("\n  y2in = %5.3f  y10 = %5.3f  y20 = %5.3f\n",
              y2in,y10,y20));
#
# Independent variable for ODE integration
  nout=37;t0=0;tf=360;
  tout=seq(from=t0,to=tf,by=(tf-t0)/(nout-1));
# 
# Initial conditions
  y0=rep(0,2);
  y0[1]=y10;
  y0[2]=y20;
  ncall=0;
#
# ODE integration
  out=lsodes(func=ode_1,y=y0,times=tout,parms=NULL)
#
# Save and display numerical solution
  y1=rep(0,nout);y2=rep(0,nout);
  for(it in 1:nout){
    y1[it]=out[it,2];y2[it]=out[it,3];
  }
  if(ip==2){
    for(it in 1:nout){
      if(it==1){
        cat(sprintf("\n      t     y1(t)     y2(t)\n"));
      }
    cat(sprintf("%7.1f%10.4f%10.4f\n",tout[it],y1[it],y2[it]));
    }
  }
#
# Calls to ODE routine
  cat(sprintf("\n\n  ncall = %3d\n\n",ncall));
#
# Plot y1, y2
  par(mfrow=c(1,1))
  plot(tout,y1,
     xlab="t (min)",ylab="y1,y2 (micromol/ml)",
     main="y1(t),y2(t)",type="l",lwd=2,
     xlim=c(0,400),ylim=c(0.5,1));
   points(tout,y1, pch="1",lwd=2);
    lines(tout,y2,type="l",lwd=2);
   points(tout,y2, pch="2",lwd=2);