model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    a2 <- a1 + (1+d)/h - d/(a1*h)
    dR1 <- R1*(1 - R1) - a1*R1*N1
    dN1 <- a1*R1*N1 - d*N1
    dR2 <- R2*(1 - R2) - a2*R2*N2/(1+R2/h+N2/h)
    dN2 <- a2*R2*N2/(1+R2/h+N2/h) - d*N2
    return(list(c(dR1, dN1, dR2, dN2)))  
  }) 
}  

s <- c(R1=1,N1=0.01,R2=1,N2=0.01)
p <- c(a1=0.5,d=0.1,h=0.1)
f <- run(100)
newton(f)

plane(ymax=2)
plane(x=3,y=4,ymax=2,add=TRUE)
# The value of a2:
with(as.list(p),a1 + (1+d)/h - d/(a1*h))

# Make Figures
size<-5#inch
pdf("dampenA.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")
run(100)
dev.off()

pdf("dampenB.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(ymax=2,legend=FALSE)
plane(x=3,y=4,ymax=2,add=TRUE)
legend("topright",legend=names(s),col=colors,lty=1,lwd=2,cex=sizeLegend)
dev.off()