model <- function(t, state, parms) {
  state <- ifelse(state < 0, 0, state)
  N <- state
  S <- A %*% N
  dN <- b*N*(1 - S) - d*N
  return(list(dN))  
}  

n <- 10
b <- rep(1,n)
d <- rep(0,n)
p <- NULL
s <- rep(0.1,n)
names(s) <- paste("N",seq(1,n),sep="")
A <- matrix(0,nrow=n,ncol=n)
diag(A) <- 1
A[lower.tri(A)] <- runif((n*n-n)/2)*0.5
tA <- t(A)
A[upper.tri(A)] <- tA[upper.tri(tA)]
run()

# Here generating the data for the table starts

n <- 32
s <- rep(0.1,n); b <- rep(1,n); d <- rep(0,n)
names(s) <- paste("N",seq(1,n),sep="")
nsim <- 10
print(c(n=n))
for (i in c(0.3,0.4,0.5,0.75,1)) {
  barZ <- i/sqrt(n)
  div <- c()
  for (j in seq(nsim)) {
    A <- matrix(0,nrow=n,ncol=n)
    diag(A) <- 1
    A[lower.tri(A)] <- runif((n*n-n)/2)*barZ*2
    tA <- t(A)
    A[upper.tri(A)] <- tA[upper.tri(tA)]
    b <- abs(rnorm(n,1,0.1))
    d <- abs(rnorm(n,0.5,0.05))
    while (min(b-d) < 0) d <- abs(rnorm(n,0.1,0.01))
    frun <- run(1000,timeplot=TRUE)
    present <- (frun > 0.01)
    div <- c(div,sum(present))
  }
  print(c(barZ=i,meanDiv=mean(div)))
}

# Results from last simulation:

present <- (frun > 0.01)
meanN <- mean(frun[present])
meanR0present <- mean(b[present]/d[present])
meanR0absent <- mean(b[!present]/d[!present])
Apresent <- A[present,present]
Aabsent <- A[!present,!present]
meanApresent <- mean(Apresent[lower.tri(Apresent)])
meanAabsent <- mean(Aabsent[lower.tri(Aabsent)])
meanF <- (sum(present)-1)*meanApresent*meanN
print(c(maxDiversity=n, Diversity=sum(present),meanN=meanN,meanF=meanF))
print(c(meanR0present=meanR0present,meanR0absent=meanR0absent,meanApresent=meanApresent,meanAabsent=meanAabsent))
barN <- (1-mean(d)/mean(b))/(1+(n-1)*barZ)
barF <- (n-1)*barZ*barN
print(c(barZ=barZ,barN=barN,barF=barF))