# Truncate the sum at ondex p 
p=2
F<-function(u){ # Vectorial function F(x,y)=(f(x,y),g(x,y)=(0,0).
    x=u[1]
    y=u[2]
    f = x + y^(3^2)/3 + x^(3^5)/(3^2) +y^(3^7)/(3^3) + x^(3^(10))/(3^4) 
    g = y + x^(3^3)/3 + y^(3^5)/(3^2) +x^(3^8)/(3^3) + y^(3^(10))/(3^4) 
    c(f,g)
    }

F(c(-1,-1))
F(c(-1,1))
# Graph location of the roots
x<-seq(-2,2,by=0.01)
y<-seq(-2,2,by=0.01)
zx=matrix(0,length(x),length(y))
zy=zx

for (i in 1:length(x)){ for (j in 1:length(y)){
zx[i,j]=F(c(x[i],y[j]))[1]
zy[i,j]=F(c(x[i],y[j]))[2]    
}}
contour(x,y,zx,level=0,col="blue",xlim=c(-2,2),ylim=c(-2,2))
par(new=TRUE)
contour(x,y,zy,level=0,col="red",xlim=c(-2,2),ylim=c(-2,2))

h=10^{-7} # steo to numerical derivative

JF<-function(u){ # Numerical approximation to Jacobian of F(x,y)
    J= matrix(0,2,2)
    v=u; v[1]=u[1]+h
    w=u; w[1]=u[1]-h
    Gradf=c(F(v)-F(w))/(2*h)
    v=u; v[1]=u[2]+h
    w=u; w[1]=u[2]-h
    Gradg=c(F(v)-F(w))/(2*h)
    J[,1]=Gradf
    J[,2]=Gradg
    J
    }
    
v = c(1,-1) # Initial guess vector (-1,1)
v
alpha=0.0001

for (k in 1:10000){
    v = v +alpha*t(JF(v))%*%F(v)
   }

v
F(v)