model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dX <- X*(1 - X/K) - c*X^n/(1+X^n)
    return(list(c(dX)))  
  }) 
}  

model2 <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dX <- X*(1 - X/K) - C*X^n/(1+X^n)
    dC <- eps
    return(list(c(dX,dC)))  
  }) 
}  

p <- c(K=10,c=1,n=2,eps=0)
s <- c(X=10)
s <- newton(run())

continue(state=s,x="c",y="X",xmin=0.1,xmax=3,ymax=10)

after <- "state[1]<-abs(state[1]*rnorm(1,1,0.1))"

# start at steady state, try this for various values of c
p["c"] <- 2; s <- newton(run())
data <- run(750,after=after,table=TRUE)
plot(data$X[1:nrow(data)-1],data$X[2:nrow(data)],type="p")
cor(data$X[1:nrow(data)-1],data$X[2:nrow(data)])

# Now change c slowly by using model2()
s <- c(X=10,C=1)
p["eps"] <- 0.001
s <- run(odes=model2)

p["eps"] <- 0.001
run(2000,odes=model2,after=after)
s <- c(X=0.4,C=3)
p["eps"] <- -0.001
run(2000,odes=model2,after=after)

# Noise on K:
s <- c(X=10,C=1)
s <- run(odes=model2)
after <- "parms[1]<-abs(10*rnorm(1,1,0.1))" 
p["eps"] <- 0.001
run(2000,odes=model2,after=after)
s <- c(X=0.4,C=3)
p["eps"] <- -0.001
run(2000,odes=model2,after=after)

# Make one plot with noise on K
s <- c(X=10,C=1)
s <- run(odes=model2)
after <- "parms[1]<-abs(10*rnorm(1,1,0.1))" 
p["eps"] <- 0.001
data1 <- run(2000,odes=model2,after=after,table=TRUE)
s <- c(X=0.4,C=3)
p["eps"] <- -0.001
data2 <- run(2000,odes=model2,after=after,table=TRUE)
plot(data1$C,data1$X,type="l",col="red",xlab="Harvest C",ylab="Density X")
lines(data2$C,data2$X,type="l",col="blue")