def newton_method(x,f_x,gradient_f_x,hessian_f_x):
try:
    inverse = np.linalg.inv(hessian_f_x(x))
except np.linalg.LinAlgError:
    print('Hessian is non invertible')
    pass
else:
    iters = 0
    f_values = []
    gradient_f = gradient_f_x(x)
    f = f_x(x)
    f_values.append(f)
    while (np.sqrt(gradient_f.dot(gradient_f))/(1+abs(f))) >= epsilon and iters<=itmax: # condition to achieve optimality
        iters+=1
        p = -np.dot(inverse,gradient_f_x(x)) # search direction
        x = x+p # update of solution point
        gradient_f = gradient_f_x(x)
        f = f_x(x)
        f_values.append(f)
    return (x,iters,f_values)