model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    dN1 <- r1*N1*(1-(N1/K)^m1)
    dN2 <- r2*N2*(1-(N2/K)^m2)
    dN3 <- r3*N3*(1-(N3/K)^m3)
    return(list(c(dN1,dN2,dN3))) 
  }) 
}  

p <- c(r1=1,r2=1,r3=1,m1=1,m2=1,m3=1,K=1)
s <- c(N1=1,N2=1,N3=0)

run(tmax=10,tstep=0.1)

# Make Figures
size<-5#inch

pdf("logist3A.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["r2"] <- 0.1
run(50,0.01,after="if (runif(1)<0.025) state[1:2]<-abs(state[1:2]+rnorm(1,0,0.05))",ymin=0.25,legend=FALSE,lwd=1)
dev.off()

pdf("logist3B.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(50,0.01,after="if (runif(1)<0.025)parms[\"K\"]<-abs(rnorm(1,1,0.1))",ymin=0.5,legend=FALSE,lwd=1)
dev.off()

pdf("logist3C.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["r2"] <- 1
s["N3"] <- 1
p["m1"] <- 2; p["m3"] <- 0.5
data <- run(50,0.01,after="if (runif(1)<0.025) state<-abs(state+rnorm(1,0,0.05))",ymin=0.6,legend=FALSE,lwd=1,table=TRUE)
dev.off()