f <- function(t, y) -2*t*y^2
y.exact <- function(t) 1/(1+t^2)

euler <- function(t, y, h, f) f(t, y)
heun <- function(t, y, h, f) 0.5*(f(t, y) + f(t+h, y + h*f(t, y)))
rk4 <- function(t, y, h, f) {
  K1 <- f(t, y)
  K2 <- f(t + 0.5*h, y + 0.5*h*K1)
  K3 <- f(t + 0.5*h, y + 0.5*h*K2)
  K4 <- f(t + h, y + h*K3)
  1/6*(K1 + 2*(K2+K3) + K4)
}

odesolve <- function(f, t0, T, y0, n, Phi, y.exact) {
  t <- seq(t0, T, length.out=n+1)
  h <- (T-t0)/n
  y <- rep(0, n+1)
  y[1] <- y0
  for (j in 1:n) {
    y[j+1] <- y[j] + h*Phi(t[j], y[j], h, f)
  }
  cbind(t, y, y - y.exact(t))
}

error <- function(n, Phi, y.exact) {
  err <- rep(0, length(n))
  for (k in seq_along(h)) {
    y <- odesolve(f, 0, 1, 1, n[k], Phi, y.exact)
    err[k] <- max(abs(y[,3]))
  }
  err
}

n <- 2^(0:12)
h <- 1/n

svg(point=16)
oldpar <- par(mar=c(5.1,4.1,1.1,1.1))
plot(c(1e-4, 1), c(1e-16, 1), xlab="h", ylab="Fehler", log="xy", type="n")
points(h, error(n, euler, y.exact), col="magenta")
lines(h, error(n, euler, y.exact), lty=2, col="magenta")
points(h, error(n, heun, y.exact), col="darkgreen")
lines(h, error(n, heun, y.exact), lty=2, col="darkgreen")
points(h, error(n, rk4, y.exact), col="gold4")
lines(h, error(n, rk4, y.exact), lty=2, col="gold4")
legend("bottomright", legend= c("expl. Euler", "Heun", "klass. Runge-Kutta"), pch=1, lty=2, col=c("magenta","darkgreen","gold4"))
par(oldpar)
dev.off()