# Gráfico em 3 dimensões

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

f1<-function(u){   # Definição de f1
x=u[1]
y=u[2]
f1=exp(x+y)-2
f1 
}

f2<-function(u){   # Definição de f1
x=u[1]
y=u[2]
f2=sin(x)+cos(y)-1
f2
}

f<-function(u){   # Definição de f=(f1,f2)
c(f1(u),f2(u))
}

Jacf<-function(u){ # Jacobiano de f(x,y)
x=u[1]
y=u[2]
dfx= exp(x+y)*c(1,1)
dfy= c(cos(x),-sin(y))

A=matrix(0,2,2)
A[1,]=dfx; A[2,]=dfy
p=A
return(p)
} 
 

#---------------------------Método de Euler para resolver a equação u'(t)=-[Jacf(u(t))]^(-1) x f(u(0)), e(0)=(2,-2)

#---------------------------Continuacao-homotopia



t0=0  # tempo inicial
tf=1  # t final 
e0=c(1,1) # condição inicial 
n=100
h=1/n  # Tamanho do passo


tt=seq(t0,tf,by=h)

Y=matrix(0,2,length(tt))
Y[,1]=e0

for ( i in 1:(length(tt)-1)){
k1=solve(Jacf(Y[,i]),-f(Y[,1]))
k2=solve(Jacf(Y[,i]+h*k1),-f(Y[,1]))
Y[,i+1]=Y[,i]+h*(k1+k2)/2
}



print("Valores a´rpximados por homotopia, com u(0)=(1,1)")
print("u(1) aproximado"); Y[,length(tt)]
print("f(u(1)) aproximado"); f(Y[,length(tt)])

plot(Y[1,],Y[2,],col="blue",main = "u'(t)=-[Jacf(u(t))]^(-1) x f(u(0)), u(0)=(1,1)")

G<-function(s,u){f(u)-f(e0)+s*f(e0)}


Gg=0*tt
for (i in 1:length(tt)){Gg[i]=sum(abs(G(tt[i],Y[,i])))}
plot(tt,Gg,xlab="t",ylab="G(t,u(t))")