model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    N <- state[1:ncomp]
    Nleft <- rep(0,ncomp); Nright <- rep(0,ncomp);
    Nleft[2:ncomp] <- N[1:(ncomp-1)]; Nleft[1] <- N[ncomp]
    Nright[1:(ncomp-1)] <- N[2:ncomp]; Nright[ncomp] <- N[1]
    dtN <- r*N*(1-N/K) + D*(Nleft+Nright-2*N)
    return(list(dtN))  
  }) 
}  

ncomp <- 100
N <- rep(0,ncomp)
names(N) <- paste("N",seq(1,ncomp),sep="")
N[ncomp/3] <- 1
s <- N
p <- c(r=1,K=1,D=0.1)
f <- run(10)
plot(seq(ncomp),f)

data <- run(100,table=TRUE)
plot(1,1,type="n",xlim=c(1,ncomp),ylim=c(0,1),xlab="Position",ylab="Density")
for (i in seq(100)) {
  lines(seq(ncomp),data[i,2:(ncomp+1)])
  Sys.sleep(0.1)
}

# Make figure
size <- 5 #inch
pdf("fisher.pdf",width=size,height=0.5*size)
par(mar=c(2.6,2.6,0.2,0.2),mgp=c(1.5,0.5,0)) 
plot(1,1,type="n",xlim=c(1,ncomp),ylim=c(0,1),xlab="Position",ylab="Density")
for (i in seq(60))
  lines(seq(ncomp),data[i,2:(ncomp+1)])
dev.off()