# differential equation y'=sin(t)^2 * y
dy <- function(t, y) sin(t)^2 * y

# exact solution 
exact <- function(t) 2 * exp(0.5*(t - sin(t)*cos(t)))

# euler's method
euler <- function(t, y, h, fun) {
  y1 <- y + h*fun(t, y)
  return(c(t + h, y1))
}

# heun's method
heun <- function(t, y, h, fun) {
  yp <- y + h*fun(t, y)
  y1 <- y + 0.5*h * (fun(t, y) + fun(t+h, yp))
  return(c(t + h, y1))
}

# classical Runge–Kutta method
runge <- function(t, y, h, fun) {
  y0 <- fun(t, y)
  ya <- fun(t+h/2, y + h/2*y0)
  yb <- fun(t+h/2, y + h/2*ya)
  yc <- fun(t+h, y + h*yb)
  
  y1 <- y + h/6*(y0 + 2*(ya+yb) + yc)
  return(c(t + h, y1))
}

# step size = 0.5, last value = 5
h <- 0.5
niter <- 5/h
run <- eul2 <- eul <- heu <- data.frame(t=0, y=exact(0))

for(i in seq_len(niter)+1) {
  eul[i, ] <- euler(t=eul$t[i-1], y=eul$y[i-1], h=h, fun=dy)
  heu[i, ] <- heun (t=heu$t[i-1], y=heu$y[i-1], h=h, fun=dy)
  run[i, ] <- runge(t=run$t[i-1], y=run$y[i-1], h=h, fun=dy)
}

# euler's method with reduced step size
h <- 0.25
niter <- 5/h
for(i in seq_len(niter)+1) {
  eul2[i, ] <- euler(t=eul2$t[i-1], y=eul2$y[i-1], h=h, fun=dy)
}

# evaluating exact solution at 
t <- seq(0, 5, 0.1)

# concatenating the methods into a data.frame
odesolve <- rbind(data.frame(t=t, y=exact(t), method="Exact Solution"),
                  data.frame(run,  method="Runge-Kutta method"),
                  data.frame(heu,  method="Heun's method"),
                  data.frame(eul2, method="Euler's method (reduced step size)"),
                  data.frame(eul,  method="Euler's method"))

# translating into german
odesolve$method <- factor(odesolve$method, 
                          levels=c("Exact Solution", "Runge-Kutta method", 
                                   "Heun's method", 
                                   "Euler's method (reduced step size)", 
                                   "Euler's method"),
                          labels=c("Exakte Lösung", "Klassisches Runge-Kutta", 
                                   "Heun", "Euler (halbe Schrittweite)", 
                                   "Euler"))

library(ggplot2)
p <- ggplot(odesolve, aes(x=t, y=y, col=method)) +   geom_line() + 
  geom_point(data=subset(odesolve, as.numeric(method)!=1)) +
  scale_color_discrete("") + 
  theme_bw() + theme(legend.position=c(0.02, 1), legend.justification=c(0, 1))

ggsave("runge-kutta.svg", width=8, height=6, plot=p)