x=c(-5, -4.5, -4,-3.5,-3,-2.5,-2,-1.5,-1,-0.5,0,0.5,1,1.5,2,2.5,3,3.5,4,4.5,5)
y=c(0.16, 0.17, 0.18, 0.2, 0.2, 0.26, 0.3, 0.46, 0.7,  1.4, 2, 1.4, 0.8, 0.5, 0.3, 0.26, 0.2, 0.2, 0.2, 0.17, 0.16)

plot(x,y,col="red")

# Agora cálculos mostram que a função g(u)= |Au-y|^2 / 2
# em que A é a matriz e é um vetor, como montamos no que segue.
n=length(x)
A=matrix(0,n,3); v=0*x;
for ( i in 1:n){A[i,]=c(1,x[i]^2,y[i]*x[i]^2)}

f<-function(u){A%*%u-y} # f(u)=Au-y

g<-function(u){  #  |Au-v|^2 / 2
w=A%*%u-y; t(w)%*%w/2}

# Agora é conhecido que o jacobiano de f(u) é A e o gradiente de g(u) é dado por A^*(Au-z).
#
# 
Jacf<-function(u){ A }   # função jacobiano de f(u).


gradg<-function(u){t(A)%*%f(u)}


#---- Início da rotin para o método de Euler.


ZeroEuler<-function(u0,m){
Y=u0
for ( i in 1:m){
B=Jacf(Y)
w=t(B)%*%f(Y)
bb=abs( sum( (( t(B)%*%B )%*%w )*w ) )
h=sum(w*w)/bb
if (bb!=0){Y=Y-h*w }
}
Y}

# Teste de aproximação de solução do problema de minimização. 
u0=c(0,0,0);print("Condição inicial");u0; 
m=10000
print("Teste de g(u0)");g(u0)                       # Teste de g(u0).

u=ZeroEuler(u0,m) ; print("Aproximação"); u     # Aproximadamente u(t)
print("Teste de minímo aproximado");g(u)                        # Teste de minímo aproximado

# Definição da aproximação de f(u).

fap<-function(s){(u[1]+u[2]*s^2)/(1-u[3]*s^2)} 

curve(fap,-5,5,col="blue",ylim=c(0,2))
points(x,y,col="red")        # Teste de ajuste.

plot(x,y-fap(x),'l')

# Se quiser pode aumentar t e n ou escolher outra condição inicial 
# para se convencer do resultado obtido.