# 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){
    x=u[1]
    y=u[2]
    z=u[3]
    p=matrix(0,3,3)
    df1=c(6*y^4+2+12*x+36*x^2*y^2+24*x^2+30*x^4, -12*y+24*x*y^3+24*x^3*y)
    df2=c(-12*y+24*x*y^3+24*x^3*y, 30*y^4-24*y^2+2-12*x+36*x^2*y^2+6*x^4 )
    df3=c(12*z^2-8,0)
    p[1,]=df1
    p[2,]=df2
    p[3,]=df3
    p
}

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


# Método de Euler-Newton para a Equação Jf(u(t))u'(t)=-alpha f(u(t))
alpha=9.8

ZeroEuler<-function(u0,t,n){
    u=matrix(0,3,n+1)
    print(u)
    u[,1]=u0
    print(u)
    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("Resolução equação ")
u0=c(0,0.5,sqrt(2)); t=1; n = 10
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=10")