model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dN1 <- N1*(b1+beta1*N2/(h1+N2) - d1*(1+e1*N1))
    dN2 <- N2*(b2+beta2*N1/(h2+N1) - d2*(1+e2*N2))
    return(list(c(dN1, dN2)))  
  }) 
}  

p <- c(b1=0.1,b2=0.1,beta1=1,beta2=1,d1=0.5,d2=0.5,e1=1,e2=1,h1=0.1,h2=0.1)
s <- c(N1=1,N2=1)
plane(xmax=2,ymax=2)