model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    f = R*W/(h*c+h*W+W*R)
    dR <- R*(1-R/K) - b*W*f
    dW <- b*W*f - d*W
    return(list(c(dR, dW)))  
  }) 
} 

p <- c(b=1,c=1,d=0.5,h=0.15,K=1)
s <- c(R=1,W=1)

#size <- 5
#pdf("wolvesA.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")
plane(xmax=1.1,main="(a)")
newton(state=c(R=0.5,W=0.5),plot=TRUE,jacobian=TRUE)
newton(state=c(R=0.9,W=0.2),plot=TRUE,jacobian=TRUE)
newton(state=c(R=1,W=0),plot=TRUE,jacobian=TRUE)
newton(state=c(R=0,W=0),plot=TRUE,jacobian=TRUE)
#dev.off()

#pdf("wolvesB.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["K"] <- 2
plane(xmax=2.1,main="(b)")
newton(state=c(R=0.5,W=0.5),plot=TRUE,jacobian=TRUE)
newton(state=c(R=2,W=0.2),plot=TRUE,jacobian=TRUE)
newton(state=c(R=1,W=0),plot=TRUE,jacobian=TRUE)
newton(state=c(R=0,W=0),plot=TRUE,jacobian=TRUE)
#dev.off()