model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dA <- r*A*(1-A/K) - p*Z*A/(1+A)
    dZ <- p*Z*A/(1+A) - m*Z - F*Z^2/(h^2+Z^2);
    return(list(c(dA, dZ)))  
  }) 
}  

p <- c(r=1,K=10,h=1,p=0.5,m=0.4,F=0.15)
s <- c(A=10,Z=0.01)
plane(xmax=20,ymax=10)