10.9 Using LSODA algorithm

require(wnl) ; require(deSolve)
dPK12 = read.csv("C:/G/Rt/PK/PK12.csv", skip=1) ; colnames(dPK12) = c("TIME", "DV")
dPK12
Dpo = 2500

PKde = function(t, y, p)
{
  if (t >= 60 & t < 75 ) RateIn = 33.33333 # 500/15
  else RateIn = 0

  if (t < p["Tlag"]) Ka = 0
  else Ka = p["Ka"]

  dy1dt = -Ka*y[1]                                # gut
  dy2dt = (RateIn + Ka*y[1] - p["Cl"]*y[2] - p["Cld"]*y[2] + p["Cld"]*y[3])/p["Vc"]
  dy3dt = (p["Cld"]*y[2] - p["Cld"]*y[3])/p["Vt"] # peripheral
  return(list(c(dy1dt, dy2dt, dy3dt)))
}
# Test PKde
Times = c(0, dPK12[,"TIME"])
iTime = 2:length(Times)
y = lsoda(y=c(0.0464*Dpo, 0, 0), times=Times, func=PKde, 
          parms=c(Ka=0.104, Vc=0.12, Vt=0.2758, Cl=0.014566, Cld=0.02, Tlag=4.8694))
plot(dPK12[,"TIME"], dPK12[,"DV"]) ; lines(y[,1], y[,3])