model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    f <- pmax(0,K - e1*N1 - e2*N2)
    dN1 <- N1*(b1*f - d1)
    dN2 <- N2*(b2*f - d2)
    return(list(c(dN1, dN2)))  
  }) 
}  

model2 <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    f <- pmax(0,K - e1*N1 - e2*N2)
    dN1 <- N1*(b1*f/(h1+f+c*N1) - d1)
    dN2 <- N2*(b2*f/(h2+f+c*N2) - d2)
    return(list(c(dN1, dN2)))  
  }) 
}

model3 <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dR <- s - d*R - e1*b1*N1*R/(h1+R) - e2*b2*N2*R/(h2+R) 
    dN1 <- N1*(b1*R/(h1+R) - d1)
    dN2 <- N2*(b2*R/(h2+R) - d2)
    return(list(c(dR, dN1, dN2)))  
  }) 
}

p <- c(K=1,b1=1,b2=0.9,d1=0.1,d2=0.1,e1=0.01,e2=0.009)
s <- c(N1=0.0,N2=0.01)
run()
p["b2"] <- 0.25
plane(xmax=100,ymax=100,portrait=TRUE,grid=3)

p <- c(K=1,b1=1,b2=1,d1=0.1,d2=0.1,e1=0.015,e2=0.01,h1=0.1,h2=0.15,c=0.01)
run(odes=model2,legend=FALSE)
#p["h2"] <- 2
plane(odes=model2,xmax=100,ymax=100,legend=FALSE)

s <- c(R=1, N1=0.01,N2=0.01)
p <- c(s=1,d=1,b1=1,b2=1,d1=0.1,d2=0.1,e1=0.015,e2=0.01,h1=0.1,h2=0.15)
run(odes=model3,legend=FALSE)
p["h2"] <- 2
plane(odes=model3,xmax=100,ymax=100,legend=FALSE)

s <- c(N1=0.0,N2=0.01)

# Make Figures
# size<-5#inch
# pdf("compA.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 <- c(K=1,b1=1,b2=0.9,d1=0.1,d2=0.1,e1=0.01,e2=0.009)
# run(legend=FALSE)
# dev.off()
#
# pdf("compB.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 <- c(K=1,b1=1,b2=0.25,d1=0.1,d2=0.1,e1=0.01,e2=0.009)
# plane(xmax=100,ymax=100,portrait=TRUE,grid=3)
# dev.off()