pred2=function(r1, r2, b1, b2, x0, y0, n) { tdeath = rexp(n) x = y = rep(0, n + 1) time = rep(0, n + 1) x[1] = x0 y[1] = y0 time[1] = 0 for(i in (2:(n + 1))) { B1 = r1 * x[i - 1] B2 = b2 * x[i - 1] * y[i - 1] D1 = b1 * x[i - 1] * y[i - 1] D2 = r2 * y[i - 1] r = B1 + D1 + B2 + D2 if(r < min(c(r1, b2, b1, r2))) { x = x[1:i - 1] y = y[1:i - 1] time = time[1:i - 1] return(list(x = x, y = y, time = time)) } else { u = runif(1) x[i] = x[i - 1] y[i] = y[i - 1] if(u < B1/r) { x[i] = x[i - 1] + 1 } else if(u < (B1 + D1)/r) { x[i] = x[i - 1] - 1 } else if(u < (B1 + D1 + B2)/r) { y[i] = y[i - 1] + 1 } else { y[i] = y[i - 1] - 1 } time[i] = time[i - 1] + tdeath[i - 1]/r } } return(list(x = x, y = y, time = time)) }