model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    V <- I
    dtC <- s - dC*C - beta*C*V
    dtI <- beta*C*V - (dI + k*E)*I
    dtE <- a*E*V - dE*E
    return(list(c(dtC, dtI, dtE)))
  }) 
}  

model2 <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    V <- I
    dtC <- s - dC*C - beta*C*V/(1+C/hb+V/hb)
    dtI <- beta*C*V/(1+C/hb+V/hb) - dI*I - k*E*I/(1+E/hk+I/hk)
    dtE <- a*E*V/(1+E/ha+V/ha) - dE*E
    return(list(c(dtC, dtI, dtE))) 
  }) 
}  

# To prevent problems with T meaning TRUE we call target cells C

p <- c(s=0.1,dC=0.1,dI=0.1,k=1,beta=1.6,a=1,dE=0.01,hb=1,hk=1,ha=1)
s <- c(C=1,I=0,E=1e-3)
newton(s)
s["I"] <- 0.01
run(log="y")
run(odes=model2,log="y")

# Make HIV figures
size<-5#inch

pdf("hivA.pdf",width=size,height=size)
par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
run(log="y")
dev.off()

pdf("hivB.pdf",width=size,height=size)
par(mar=c(0.1,0.1,0.1,0.1),xaxt="n",yaxt="n",ann=F,xaxs="i",yaxs="i")
run(odes=model2,log="y")
dev.off()