# Resolver a equação f(x,y)=(f1(x,y),f2(x,y),f3(x,y))=(0,0.5,sqrt(2)).

f1<-function(u){   # Definição de f1
    x=u[1]
    y=u[2]
    z=u[3]
    return(2-6*y^2+6*x*y^4+2*x+6*x^2+12*x^3*y^2+8*x^3+6*x^5)
}

f2<-function(u){   # Definição de f2
    x=u[1]
    y=u[2]
    z=u[3]
    p=6*y^5-8*y^3+2*y-12*x*y+12*x^2*y^3+6*x^4*y
    return(p)
}

f3<-function(u){    # Definição de f3
    x=u[1]
    y=u[2]
    z=u[3]
    return(4*z^3 - 8*z)
}

#--- Minimo de |f(u)|^2 para resolver equação f(x)=0.

f<-function(u){      # Definição de função f(u), para resolver a equação f(u)=0
    c( f1(u), f2(u),f3(u))
}

g<-function(u){
    t(f(u)) %*% f(u)/2
} 

Jf<-function(u){ #Calcular o Jacobiano
    h = 10^(-5)
    v=u; v[1]=u[1]+h; v1=u; v1[1]=u[1]-h; 
    dfx= f(v)-f(v1)
    v=u; v[2]=u[2]+h; v1=u; v1[2]=u[2]-h; 
    dfy= f(v)-f(v1)
    v=u; v[3]=u[3]+h; v1=u; v1[3]=u[3]-h; 
    dfz= f(v)-f(v1)
    p=matrix(0,3,3)
    p[1,]=dfx
    p[2,]=dfy
    p[3,]=dfz
    p = p/(h*3)
    return(p)
}

gradg<-function(u){
    p=t(Jf(u))%*%f(u)
    t(p)
}

# Metodo de Euler-Newton para a Equacao Jf(u(t))u'(t)=-alpha f(u(t))
alpha=9.8

ZeroEuler<-function(u0,t,n){
    u=matrix(0,3,n+1)
    u[,1]=u0
    tj=t/n
    for (i in 1:n){
        h=solve( Jf(u[,i]),-tj*alpha*f(u[,i]))    # Resolve o sistema linear Jf(u)h=-tj alpha f(u)
        u[,i+1]=u[,i]+h
    }
    u
}

# Teste 

print("Resolucao equacao ")
u0=c(0,0.5,sqrt(2)); t=1; n = 9
u=ZeroEuler(u0,t,n) ; u [,n+1]     # Aproximadamente u(1)
#f(u[,n+1])
#g(u[,n+1])

#x=u[1,]; y=u[2,]
#plot(x,y,'l',col="blue", main="Euler-Newton para Jf(u(t))u'(t)=-alpha f(u(t)), com n=9")