model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {  
    vR <- ifelse(R-f < 0, 0, R-f)
    dR <- r*R*(1 - R/K) - a*vR*N/(h+vR)
    dN <- c*a*vR*N/(h+vR) - delta*N    
    return(list(c(dR, dN)))  
  }) 
}  

p <- c(r=1,K=1,h=0.1,a=0.5,c=1,delta=0.4,f=0.05)
s <- c(R=1,N=0.01)
plane()