par(mar=c(2.6,2.6,1.6,0.2),mgp=c(1.5,0.5,0))

model <- function(t, state, parms){
  state <- ifelse(state < 0, 0, state)
  with(as.list(c(state,parms)),{
    R <- state[1:nr]
    LV <- A %*% R
    dR <- 1e-5 + r*R*(1-LV/K)
    N <- state[(nr+1):(nr+nn)]
    SR <- S %*% R
    SN <- sapply(seq(nn),function(i){S[i,]*N[i]/(h+SR[i])})
    dR <- dR - g*R*rowSums(SN)
    dN <- (e*g*SR/(h+SR) - d)*N
    return(list(c(dR,dN)))  
  }) 
} 

nr <- 12; nn <- 8; eps <- 0; w <- 0.1
p <- c(r=0.5,K=10,d=0.15,g=0.4,h=2,e=0.6)
A <- matrix(runif(nr*nr,1+eps-w,1+eps+w),nrow=nr,ncol=nr)
diag(A) <- 1
R <- runif(nr); names(R) <- paste("R",seq(1,nr),sep="")
S <- matrix(runif(nr*nn,0,1),nrow=nn,ncol=nr)
SR <- S %*% R
N <- runif(nn); names(N) <- paste("N",seq(1,nn),sep="")
s <- c(R,N)
f <- run(1000)

Rpresent <- which(f[1:nr] > 0.01)
Npresent <- which(f[(nr+1):(nr+nn)] > 0.01)
print(c(Rpresent=length(Rpresent),Npresent=length(Npresent),nr=nr,nn=nn))
print(f[Rpresent])
print(f[nr+Npresent])

Rabsent <- which(f[1:nr] < 0.01)
Nabsent <- which(f[(nr+1):(nr+nn)] < 0.01)
print(c(Apresent=mean(A[Rpresent,]),Aabsent=mean(A[Rabsent,])))
print(c(Spresent=mean(S[Npresent,]),Sabsent=mean(S[Nabsent,])))