model <- function(t, state, parms) {
  state <- ifelse(state < 0, 0, state)
  with(as.list(c(state,parms)), {
    dH <- r*H*(1 - H/K) - Q  
    return(list(dH))  
  }) 
}  

p <- c(r=0.1,K=1,Q=0)
s <- c(H=0.1)
run(200,tstep=1/52,after="parms[\"K\"]<-abs(rnorm(1,mean=1,sd=0.2))",xlab="Time in years")

p <- c(r=0.1,K=1,Q=0.025)
s <- c(H=1)
run(200,xlab="Time in years")
lines(c(0,200),c(0.5,0.5),lty=2)
run(200,tstep=1/52,after="parms[\"K\"]<-abs(rnorm(1,mean=1,sd=0.2))",add=TRUE)

model <- function(t, state, parms) { 
  state <- ifelse(state < 0, 0, state) 
  with(as.list(c(state,parms)), {
    dH <- r*H*(1 - H/K) - f*H
    return(list(dH)) })
}
p <- c(r=0.1,K=1,f=0.05) 
run(400,tstep=1/52,after="parms[\"K\"]<-abs(rnorm(1,mean=1,sd=0.5))",xlab="Time in years")

# size <- 5
# pdf("parabola1.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")
# curve(N*(1-N),xname="N",ylim=c(0,0.3),lwd=2,col="red")
# dev.off()
# pdf("parabola2.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")
# curve(N*(1-N)-N/2,xname="N",from=0,to=1,ylim=c(0,0.1),lwd=2,col="red")
# dev.off()