model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dS <- a - d*S - beta*S*I - e*S
    dI <- beta*S*I + e*S - delta*I
    return(list(c(dS,dI)))  
  }) 
}  

p <- c(beta=4e-2,a=0.01,d=0.01,delta=0.02,e=1e-4)
s <- c(S=1,I=1e-3)
with(as.list(p),(beta/delta)*a/d)# R0
with(as.list(p),a/(d+e))
with(as.list(p),delta/beta)
plane()
p["beta"] <- 0.01
plane(xmax=2)

# make figures
size <- 5 #inch
pdf("aidsA.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")
p["beta"] <- 0.04
plane()
dev.off()

size <- 5 #inch
pdf("aidsB.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")
p["beta"] <- 0.01
plane(xmax=2)
dev.off()