a <- 0 #interval domain
b <- 10
h <- 0.01 #step size
t <- seq(a,b,by=h) #time sequence
N <- (b-a)/h #steps needed
t0 <- 0 #initial condition
x0 <- 1

x <- rep(0,N+1)
x[1] <- x0

field <- function(t,x){
  return(x) #exponential
  #return(0.5*(x*x)-2*x+1) #polynomial
  #return(x*(sin(t))^2) #linear trigonometric
  #return(sin(t)*sin(x)) #trigonometric
  #return(1/(1+(x*x))) #arctangent
  #return((cos(6*t))/(1+t+(x*x)))
}

#Euler method
for(n in 1:N){
  x[n+1] <- x[n] + h*field(t0+n*h, x[n])
}

plot(t,x, type="l", col="red")

#Runge-Kutta method
for(n in 1:N){
  k1 <- h*field(t0+n*h, x[n])
  k2 <- h*field(t0+n*h+(h/2), x[n]+k1/2)
  k3 <- h*field(t0+n*h+(h/2),x[n]+k2/2)
  k4 <- h*field(t0+n*h+h, x[n]+k3)
  x[n+1] <- x[n] + (1/6)*(k1+2*k2+2*k3+k4)
}

lines(t,x, type="l", col="blue")