h <- 0.1 # timestep
Tmax <- 4

tlist <- seq(0, Tmax, length = Tmax/h + 1)
xlist <- rep(0, length(tlist))

# initialize
x0 <- 1
xlist[1] <- x0

for (n in 1:(length(tlist)-1))
    xlist[n+1] <- xlist[n] + h/2*xlist[n] + h/2*(xlist[n] + h*xlist[n])

plot(tlist, xlist, xlab = "t", ylab = "x(t)")

tlistfine <- seq(0, Tmax, length = 10000)
lines(tlistfine, exp(tlistfine), col = "red",
     xlab = "t", ylab = "y(t)") # exact solution x = exp(t)

legend("topleft", c("Heun Approx", "Exact Solution"),
       col = 1:2, pch = 1:0, lty = 0:1)