model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    k2 <- (k1^2*s^3)/(d*(d*k1 + s)*(d*k1 + s - k1*s))
    dN1 <- s*(1-N1/k1) - d*N1
    dN2 <- s - d*(1+N2/k2)*N2
    return(list(c(dN1,dN2)))  
  }) 
}  

s <- c(N1=0,N2=0)
p <- c(s=1,k1=0.25,d=1)
run(2.5,0.01,legend=FALSE)

s <- c(N1=0.2,N2=0.2)
newton(s,jacobian=TRUE)
data <- run(50,0.01,after="if (runif(1)<0.025) state<-abs(state+rnorm(1,0,0.05))",legend=FALSE,lwd=1,table=TRUE)
sd(data$N1)
sd(data$N2)

# Make Figures
size<-5#inch
ymin <- min(c(data$N1,data$N2))
ymax <- max(c(data$N1,data$N2))

pdf("sourceA.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")
s <- c(N1=0,N2=0)
run(2.5,0.01,ymax=ymax,legend=FALSE)
dev.off()

pdf("sourceB.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")
timePlot(data,lwd=1,ymin=ymin,ymax=ymax,legend=FALSE)
lines(c(0,100),c(0.5,0.5))
dev.off()

pdf("sourceC.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")
plot(1,1,type="n",xlim=c(0,0.25),ylim=c(-0.1,1.1),xlab="N",ylab="f(N)")
with(as.list(p),curve(s*(1-N/k1)-d*N,from=0,to=1.1,xname="N",col="red",lwd=2,add=TRUE))
k2 <- (k1^2*s^3)/(d*(d*k1 + s)*(d*k1 + s - k1*s))
with(as.list(p),{
  k2 <- (k1^2*s^3)/(d*(d*k1 + s)*(d*k1 + s - k1*s))
  curve(s-d*(1+N/k2)*N,from=-0.1,to=1.1,xname="N",col="blue",lwd=2,add=TRUE)
  })
lines(c(0,0.25),c(0,0))
dev.off()