#-----------------------------------------------------------------------------------------------
# Exemplo 4: Modelo Lotka-Volterra D^alpha Y= f(t, Y).

alpha = 0.9

f<-function(t){2*c(t[1]-3*t[1]*t[2]/2,-2*t[2]+t[1]*t[2]/2)}

x0=0
xn=50

y0=c(10,1.50)

n=1000000

h=0.0000000000001 
y=matrix(0,2,n)
y[,1]=y0
y[,2]=alpha*f(y[,1])*h^alpha/gamma(2+alpha)+ f( y[,1]+h^alpha*f(y[,1])/gamma(1+alpha))*h^alpha/gamma(2+alpha)+y[,1]

for (j in 3:n){ 
p=y[,1]+ ( (j-1)^(alpha+1)-(j-alpha-1)*j^(alpha) )*f(y[,1])*h^alpha/gamma(2+alpha)+ f(y[,j-1]+h^alpha*f(y[,j-1])/gamma(1+alpha) )*h^alpha/gamma(2+alpha)

for ( i in 2:(j-1)){
p=p+( (j-i+1)^(alpha+1)-2*(j-i)^(alpha+1)-(j-i-1)^(alpha+1) )*f(y[,i])*h^alpha/gamma(2+alpha)
}
y[,j]=p}
y

y1=y[1,]
y2=y[2,]
# Solução 
plot(y1,y2,col="red",'l')