Runge-Kutta in R

runge_kutta.R

# x is xo, y is y0, x1 is x-value for which you want y, n is number of steps,
# f is the function f in the DE dy/dx = f(f,)
rungeKutta <- function(x, y, x1, n, f) {
  h <- (x1 - x) / n
  X <- numeric(n + 1)
  Y <- numeric(n + 1)
  X[1] <- x
  Y[1] <- y
  for ( i in 1:n ) {
    xhalf <- x + 0.5 * h
    xnew <- x + h
    k1 <- f(x, y)
    u1 <- y + 0.5 * h *k1
    k2 <- f(xhalf, u1)
    u2 <- y + 0.5 * h * k2
    k3 <- f(xhalf, u2)
    u3 <- y + h * k3
    k4 <- f(xnew, u3)
    k <- (k1 + 2 * k2 + 2 * k3 + k4) / 6
    x <- xnew
    y <- y + k * h
    X[i + 1] <- x
    Y[i + 1] <- y
  }
  data.frame(x = X, y = Y)
}

# this is the f(x,y) in the DE:  dy/dx = f(x, y)
myDEFunc <- function(x, y) {
  y
}


res <- rungeKutta(x = 0, y = 1, x1 = 1, n = 100, f = myDEFunc)