model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    # For Z=1 the model uses the Hill function
    # For Z=0 the exponential function
    dR <- r*R*(1 - R/K) - Z*a*R*N/(h+R) - (1-Z)*a*N*(1-exp(-log(2)*R/h)) 
    dN <- Z*c*a*R*N/(h+R) + (1-Z)*c*a*N*(1-exp(-log(2)*R/h)) - delta*N
    return(list(c(dR, dN)))  
  }) 
}  

p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.35)
s <- c(R=1,N=0.01)
Z <- 1;plane()
Z <- 0;plane(add=TRUE)