F<-function(s,u){
    1+u^2
}

# Condição inicial y(0)=1, y'=1+y^2=F(t,y). Calcular y(-1.5)

beta=1; delta=beta; alpha_2=1/(2*beta); alpha_1=1-alpha_2

n=5
a=0; b=-1.5
h=(b-a)/n 

t=seq(a,b,by=h)
y=0*t 

y[1]=1 # condição incial y(0)=1

for (i in 1:n){ # Runge-Kutta2
k1=F(t[i],y[i])
k2=F(t[i]+beta*h,y[i]+delta*h*k1)
 y[i+1]=y[i]+h*(alpha_1*k1+alpha_2*k2)  
}

v=0*y; v[1]=1 # condição incial y(0)=1
for (i in 1:n){ # Euler
k1=F(t[i],v[i])
v[i+1]=v[i]+h*k1
}

y[n+1]

w=0*t 

w[1]=1 # condição incial y(0)=1

for (i in 1:n){ # Runge-Kutta4
k1=F(t[i],w[i])
k2=F(t[i]+h/3,w[i]+h*k1/3)
k3=F(t[i]+2*h/3,w[i]-h*k1/3+h*k2)
k4=F(t[i]+h,w[i]+h*(k1-k2+k3))
 w[i+1]=w[i]+h*(k1+3*(k2+k3)+k4)/8
}

# Sabemos que y(t)=tan(t+pi/4)
g<-function(s){tan(s+pi/4)}
curve(g,-1.5,0,col="blue")
points(t,y)
points(t,v,col="red")

points(t,w,col="green")