model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dS <- p*S*(1-S/K) - d*S
    dD <- (p/K)*S^2 - delta*D
    return(list(c(dS,dD)))  
  }) 
}  

s <- c(S=1,D=0)
p <- c(p=1,d=0.01,delta=0.1,K=1)

run()
plane(xmax=2,ymax=20)

# Make Figures
# size<-5#inch
# 
# pdf("stem.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(xmax=2,ymax=20)
# newton(c(S=1,D=10),plot=TRUE)
# dev.off()