midpoint <- function(dydt, t0, y0, h, n, ...) {
  #mispoint's method applied for system of ODEs with grad dydt and
  #initial cond. y at t0
  #n steps of size h
  y <- matrix(0, nrow = n + 1, ncol = length(y0))
  y[1,] <- y0
  t <- t0
  for (i in 1:n) {
    y2 <- y[i,] + h/2*dydt(t, y[i,], ...)
    y[i+1,] <- y[i,] + h*dydt(t + h/2, y2, ...)
    t <- t + h
  }
  return(y)
}

LV <- function(t = NULL, y, a, b, g, e, K = Inf)
  c(a*y[1]*(1 - y[1]/K) - b*y[1]*y[2], g*b*y[1]*y[2] - e*y[2])

par(mfrow = c(2,2))
for (h in c(10, 5, 2)) {
  y <- midpoint(LV, 0, c(100, 50), h, 1000/h,
             a = 0.05, K = Inf, b = 0.0002, g = 0.8, e = 0.03)
  plot(y[,1], y[,2], type = "p",
       xlab = "prey", ylab = "pred", main = paste("midpoint, h=", h))
}
plot(seq(0, 1000, h), y[,1], type = "l", xlab = "time",
     ylab = "prey solid pred dashed", main = paste("midpoint, h=", h))
lines(seq(0, 1000, h), y[,2], lty = 2)