Нужно решить систему нелинейных уравнений методом ньютона. Что не так в коде? [python]

Понятно почему выдаёт ошибку, но не понятно почему в экспоненте накапливается огромное число.

Вот такая система

График системы:

Консоль:

import numpy as np

### exp(x*y^2) + (x)^2 - 120 = 0
### (x-y)^2 - x*y - 20 = 0

def jacobian_exercise(x,y):
    return [[np.exp(x*y**2)*(2*x*y+x**2),  2*x],
            [0,  -(x+y)]]

jotinha  = (jacobian_exercise(1,2))

def function_exercise(x,y):
    return [np.exp(x*y**2) + (x)**2 -120,      (x-y)**2 - x*y -20]

bezao = (function_exercise(1,2))

def x_delta_by_gauss(J,b):

    return np.linalg.solve(J,b)

#x_delta_test = x_delta_by_gauss(jotinha,bezao)

def x_plus_1(x_delta,x_previous):

    x_next = x_previous + x_delta

    return x_next

#print (x_plus_1(x_delta_test,[1,2,3]))

def newton_method(x_init):

    first = x_init[0]

    second = x_init[1]

    jacobian = jacobian_exercise(first, second)


    vector_b_f_output = function_exercise(first, second)

    x_delta = x_delta_by_gauss(jacobian, vector_b_f_output)

    x_plus_1 = x_delta + x_init

    return x_plus_1

def iterative_newton(x_init):

    counter = 0

    x_old = x_init

    x_new = newton_method(x_old)


    diff = np.linalg.norm(x_old-x_new)


    while diff>1e-5:

        counter += 1

        x_new = newton_method(x_old)

        diff = np.linalg.norm(x_old-x_new)

        x_old = x_new

    convergent_val = x_new

    return convergent_val

print (iterative_newton([1,2]))