from sympy import *

x, y = var('x, y')

vars = [x, y]
u0 = Matrix([-0.2, 0.1])
f = Lambda(vars, E**(-x**2 - y**2) * (2 * x**2 + 3 * x**2))
# f = Lambda(vars, (x - 2)**4 + (y - 3)**4)
# f = Lambda(vars, 100 * (y - x**2)**2 + (1 - x)**2)

H = Matrix([[f(*vars).diff(x1, x2) for x1 in vars] for x2 in vars])
G = Matrix([f(*vars).diff(x) for x in vars])

def newton(x0, f, H, G, eps):
    print("H = ", H)
    print("G = ", G)
    Hi = H.inv()
    steps = 0
    while True:
        d = dict(zip(vars, x0))
        x1 = x0 - Hi.subs(d) * G.subs(d)
        print(steps, f(*x1), x1)
        if (x0 - x1).norm() < eps:
            break
        x0 = x1
        steps += 1

newton(u0, f, H, G, 0.001)