Cacher is the code snippet organizer for pro developers

gillespie <- function(N, n, ...) {
  #N is a list of 4 elements: initial markings, pre matrix, post matrix, and hazard function
  #n is the number of events to simulate
  timeInstant <- 0 #A time instant variable, starting as 0
  state <- N$M #system state variable, initialized to initial state
  S <- t(N$Post - N$Pre) #stochiometric matrix
  u <- nrow(S) #number of places
  v <- ncol(S) #number of transitions
  times <- vector("numeric", n)
  stateMat <- matrix(ncol = u, nrow = n + 1)
  stateMat[1, ] <- state
  for (i in 1:n){
    h <- N$h(state, ...)
    timeInstant <- timeInstant + rexp(1, sum(h))
    j <- sample(v, 1, prob=h)
    state <- state + S[ ,j]
    times[i] <- timeInstant
    stateMat[i + 1, ] <- state
  }
  return(list(times=times, stateMat=stateMat))
}
N <- list()
N$M <- c(x1 = 50, x2 = 100)
N$Pre <- matrix(c(1,0,1,1,0,1), ncol=2, byrow=TRUE)
N$Post <- matrix(c(2,0,0,2,0,0), ncol=2, byrow=TRUE)
N$h <- function(state, theta = c(th1=1, th2=0.005, th3=0.6)){
  with(as.list(c(state, theta)), {
    return(c(th1 * x1,
      th2 * x1 * x2,
      th3 * x2))
  })
}

# simulate a realisation of the process and plot it
simOutput <- gillespie(N, 10000)
plot(stepfun(simOutput$times, simOutput$stateMat[ , 2]),
     main = "Instance of Preditor Prey Population Sim",
     xlab = "time", 
     ylab = "population", 
     ylim = c(0, 500), pch = "", col = "darkgreen")
lines(stepfun(simOutput$times, simOutput$stateMat[ , 1]), col = "darkred", pch="")
text(y = 450, x = 25, labels = "prey", col = "darkgreen")
text(y = 425, x = 25, labels = "preditor", col = "darkred")