model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    Af <- pmax(0,a*A-k)
    dA <- r*A*(1 - A/K) - a*A*D
    dD <- m*E - d0*D - d1*D/(1+a*A/h)
    dE <- e*Af*D/(H+Af) - m*E
    return(list(c(dA, dD, dE)))  
  }) 
} 

qssa <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    Af <- pmax(0,a*A-k)
    E <- (e/m)*Af*D/(H+Af)
    dA <- r*A*(1 - A/K) - a*A*D
    dD <- m*E - d0*D - d1*D/(1+a*A/h)
    return(list(c(dA, dD)))  
  }) 
}  

p <- c(r=1,K=1,a=1,e=1,m=0.05,k=0.5,h=0.25,H=0.3,d0=0.01,d1=0.2)
s <- c(A=1,D=0.01)
plane(odes=qssa)
run(odes=qssa,traject=TRUE,tstep=0.1)
plane(odes=qssa)
run(state=c(A=1,D=0.01,E=0),traject=TRUE)

# Make Figures
# size<-5#inch
# 
# pdf("daphniaA.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")
# plane(odes=qssa)
# run(odes=qssa,traject=TRUE,tstep=0.1)
# dev.off()
# 
# pdf("daphniaB.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")
# plane(odes=qssa)
# run(state=c(A=1,D=0.01,E=0),traject=TRUE)
# dev.off()