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

n <- 5
s <- rep(0.1,n)
r <- rep(1,n)
p <- NULL
names(s) <- paste("N",seq(1,n),sep="")
zmean <- 0.1
z <- rnorm(n*n,zmean,zmean/10)
k <- ifelse(runif(n*n)<0.5,1,-1)
A <- matrix(k*z,nrow=n,ncol=n)
diag(A) <- 1
frun <- run(); cat(frun)
AI <- solve(A)     # Compute the inverse of A
fsol <- AI %*% r   # Use this to solve A N = r
cat(fsol)
print(c(ReturnTime=-1/max(Re(newton(frun,silent=TRUE)$values))))

n <- 144
s <- rep(0.1,n); r <- rep(1,n)
names(s) <- paste("N",seq(1,n),sep="")
nsim <- 20 
print(c(n=n))
for (i in c(0.2,0.25,0.3,0.4,0.5,1)) {
  zmean <- i/sqrt(n)
  div <- c(); fea <- c()
  for (j in seq(nsim)) {
    #r <- abs(rnorm(n,1,0.25))
    z <- rnorm(n*n,zmean,zmean/10)
    k <- ifelse(runif(n*n)<0.5,1,-1)
    A <- matrix(k*z,nrow=n,ncol=n)
    diag(A) <- 1
    #diag(A) <- abs(rnorm(n,1,0.25))
    frun <- run(500,timeplot=TRUE)
    AI <- solve(A) # Compute the inverse of A
    fsol <- AI %*% r 
    feasible <- min(fsol > 0)
    npresent <- sum(frun > 0.01)
    div <- c(div,npresent)
    fea <- c(fea,feasible)
  }
  print(c(zmean=i,meanDiv=mean(div),Pfeasible=100*sum(fea)/nsim))
}