model <- function(t, state, parms) {
  with(as.list(c(state,parms)), {
    N <- state
    S <- A %*% N  # R code for matrix x vector multiplication
    dN <- r*N*(1 - S)
    return(list(dN))  
  }) 
}  

makeMatrix <- function(alpha) {
  seqAlpha <- sapply(seq(from=0,n-1),function(i){alpha^(i^2)})
  A <- matrix(0,nrow=n,ncol=n)
  for (i in seq(n)) {
    A[i,i:n] <- seqAlpha[1:(n-i+1)]
    A[i,1:i] <- rev(seqAlpha)[(n-i+1):n]
  }
  return(A)
}

findMaxAlpha <- function(n) {
  n <<- n
  s <- rep(0.1,n)
  names(s) <- paste("N",seq(1,n),sep="")
  for (alpha in seq(0,1,0.01)) {
    A <<- makeMatrix(alpha)
    f <- newton(run(state=s,timeplot=FALSE),atol=1e-20)
    if (min(f) <= 0) return(alpha)
  }
  return(1)
}

p <- c(r=1)
n <- 3
s <- rep(0.1,n)
names(s) <- paste("N",seq(1,n),sep="")
A <- makeMatrix(0.5)
f <- newton(run(),atol=1e-20)

nspecies <- seq(9)+1
maxAlpha <- sapply(nspecies,findMaxAlpha)
cat(maxAlpha)

size <- 12 #inch
pdf("niche.pdf",width=size,height=0.4*size)
par(mar=c(2.6,2.6,1.6,0.2),mgp=c(1.5,0.5,0)) 
par(mfrow=c(1,2)) 
n <- 21
s <- rep(0.1,n)
names(s) <- paste("N",seq(1,n),sep="")
A <- makeMatrix(0.63)
f <- newton(run(timeplot=FALSE),atol=1e-20)
plot(seq(n),f,xlim=c(1,21),ylim=c(0,1),xlab="Species number",ylab="Population size",main="(a)",font.main=1)
lines(seq(n),f)
axis(1,at=seq(21,from=0))
A <- makeMatrix(0.5)
f <- newton(run(timeplot=FALSE),atol=1e-20)
points(seq(n),f,col="red")
lines(seq(n),f,col="red")
legend(14,1,c(expression(alpha==0.50),expression(alpha==0.63)), col=c("red","black"),lty=1,pch=1)


plot(nspecies,maxAlpha,xlim=c(1,10),ylim=c(0.5,1),xlab="Number of species",ylab="Niche overlap",main="(b)",font.main=1)
lines(nspecies,maxAlpha)
axis(1,at=seq(10,from=0))
par(mfrow=c(1,1)) 
dev.off()