model <- function(t, n, parms) {
  with(parms, {
    dn <- r * n  + n * (A %*% n)
    return(list(c(dn)))
  })
}

parms <- list(
  r = c(r1 = 0.1, r2 = 0.1, r3 = -0.1, r4 = -0.1),
  ## only pairwise interactions:
  A = matrix(c(0.0, 0.0, -0.2, 0.0,      # prey 1
               0.0, 0.0, 0.0, -0.1,      # prey 2
               0.2, 0.0, 0.0, 0.0,       # predator 1; eats prey 1
               0.0, 0.1, 0.0, 0.0),      # predator 2; eats prey 2
               nrow = 4, ncol = 4, byrow=TRUE)
)

times = seq(from=0, to=500, by = 0.1)

n0  = c(n1=1, n2=1, n3=2, n4=2)

system.time(
out <- ode(n0, times, model, parms)
)

plot(out)

################################################################################
windows()
model <- function(t, n, parms) {
  with(as.list(c(n, parms)), {
    dn1 <- r1 * n1 - a13 * n1 * n3
    dn2 <- r2 * n2 - a24 * n2 * n4
    dn3 <- a13 * n1 * n3 - r3 * n3
    dn4 <- a24 * n2 * n4 - r4 * n4
    return(list(c(dn1, dn2, dn3, dn4)))
  })
}

parms <- c(r1 = 0.1, r2=0.1, r3=0.1, r4=0.1,
          a13=0.2, a24 = 0.1)

system.time(
out <- ode(n0, times, model, parms)
)

plot(out)