library('Deriv')
e <- exp(1)

### Función ejemplo
f <- function(x, y) 7*x + 4*y + x*y - x^2 - y^2
exp_f = expression(7*x + 4*y + x*y - x^2 - y^2)
d_f <- D(exp_f,"x")
d_f
dx_f <- function(x, y) 7 + y - 2 * x
d_f <- D(exp_f,"y")
d_f
dy_f <- function(x, y) 4 + x - 2 * y

### Iteración 1
x1 <- 8
x2 <- 2
res_f <- f(x1,x2)
res_f

dx <- dx_f(8, 2)
dy <- dy_f(8, 2)

print(cbind(x1, x2, res_f, dx, dy))

Vfx <- function(k) (8 - 7 * k) # x1 + dx * k
Vfy <- function(k) (2 + 8 * k) # x2 + dy * k
# https://www.wolframalpha.com/input/?i=simplify

fxp1 <- function(k) -169 * k^2 + 113*k + 12

exp_fxp1 = expression(-169 * k^2 + 113*k + 12)
d_fxp1 <- D(exp_fxp1,"k")
d_fxp1
d_fxp1 <- function(k) 113 - 338 * k
# igualando a cero, k = 0.3343195
k = 0.3343195

d_fxp1(k)


xa1 <- Vfx(k)
xa2 <- Vfy(k)
fa <- f(xa1, xa2)
print(paste("k:",k,"   X':",xa1,"   Y':",xa2, "   f'", fa))

### Iteración 2
x1 <- xa1
x2 <- xa2
res_f <- f(x1,x2)
res_f

dx <- dx_f(x1,x2)
dy <- dy_f(x1,x2)

print(cbind(x1, x2, res_f, dx, dy))

Vfx <- function(k) (5.659764 + 0.355029 * k) # x1 + dx * k
Vfy <- function(k) (4.674556 + 0.3106515 * k) # x2 + dy * k
# https://www.wolframalpha.com/input/?i=simplify
fxp1 <- function(k) -0.11226 * k^2 + 0.22255 * k + 30.8891

exp_fxp1 = expression(-0.11226 * k^2 + 0.22255 * k + 30.889)
d_fxp1 <- D(exp_fxp1,"k")
d_fxp1
d_fxp1 <- function(k) 0.22255 - 0.11226 * (2 * k)
# igualando a cero, k≈0.991226
k = 0.991226
d_fxp1(k)

xa1 <- Vfx(k)
xa2 <- Vfy(k)
fa <- f(xa1, xa2)
print(paste("k:",k,"   X':",xa1,"   Y':",xa2, "   f'", fa))

### .-.-.-........-.--..-.-.-.-.-.-.-


### Mi función
f <- function(x, y) (0.65-((0.3*e^(-(y-6)/50))/12.5331))*y + (-0.65+((0.3*e^(-(y-6)/50))/12.5331))*x
exp_f = expression((0.65-((0.3*e^(-(y-6)/50))/12.5331))*y + (-0.65+((0.3*e^(-(y-6)/50))/12.5331))*x)
d_f <- D(exp_f, "x")
d_f
dx_f <- function(x, y) (-0.65 + ((0.3 * e^(-(y - 6)/50))/12.5331))
d_f <- D(exp_f, "y")
d_f
dy_f <- function(x, y) 0.3 * (e^(-(y - 6)/50) * (log(e) * (1/50)))/12.5331 * y +
  (0.65 - ((0.3 * e^(-(y - 6)/50))/12.5331)) -
  0.3 * (e^(-(y - 6)/50) * (log(e) * (1/50)))/12.5331 * x

### Iteración 1
x1 <- 5
x2 <- 6
res_f <- f(x1,x2)
res_f

dx <- dx_f(x1, x2)
dy <- dy_f(x1, x2)

print(cbind(x1, x2, res_f, dx, dy))

Vfx <- function(k) (5 - 0.6260634 * k) # x1 + dx * k
Vfy <- function(k) (5 + 0.6265421 * k) # x2 + dy * k
# https://www.wolframalpha.com/input/?i=simplify

fxp1 <- function(k) -1.25261 * (0.0239366 * e^(1/50 - 0.0125308 *k) - 0.65) *k

exp_fxp1 = expression(-1.25261 * (0.0239366 * e^(1/50 - 0.0125308 *k) - 0.65) *k)
d_fxp1 <- D(exp_fxp1,"k")
d_fxp1
d_fxp1 <- function(k) 1.25261 * (0.0239366 * (e^(1/50 - 0.0125308 * k) * (log(e) * 0.0125308))) * k + -1.25261 * (0.0239366 * e^(1/50 - 0.0125308 * k) - 0.65)
# igualando a cero, k≈-170.618
k <- -170.618

d_fxp1(k)

xa1 <- Vfx(k)
xa2 <- Vfy(k)
fa <- f(xa1, xa2)
print(paste("k:",k,"   X':",xa1,"   Y':",xa2, "   f'", fa))

### Iteración 2
x1 <- xa1
x2 <- xa2
res_f <- f(x1,x2)
res_f

dx <- dx_f(x1, x2)
dy <- dy_f(x1, x2)

print(cbind(x1, x2, res_f, dx, dy))

Vfx <- function(k) (111.8177 - 0.4428597 * k) # x1 + dx * k
Vfy <- function(k) (-101.8994 - 0.4425288 * k) # x2 + dy * k
# https://www.wolframalpha.com/input/?i=simplify

fxp1 <- function(k) 0.000215085 * (k - 645866.) - 0.0000685428 * e^(0.00885058 *k) * (k - 645866.)

exp_fxp1 = expression(0.000215085 * (k - 645866.) - 0.0000685428 * e^(0.00885058 *k) * (k - 645866.))
d_fxp1 <- D(exp_fxp1,"k")
d_fxp1
d_fxp1 <- function(k) 0.000215085 - 6.06644* (10^-7)*( e^(0.00885058 *k) )* (k - 645753)
# igualando a cero, k≈645753
k <- 645753

d_fxp1(k)

xa1 <- Vfx(k)
xa2 <- Vfy(k)
fa <- f(xa1, xa2)
print(paste("k:",k,"   X':",xa1,"   Y':",xa2, "   f'", fa))