Solving non-linear equations using Numerical methods

Question

I have been trying to solve this for some time but I cannot come up with the answer using Python: This is the task

And this is the code I have come up with but I think it's wrong because I am getting the same values for the vector x_k

import numpy as np
from scipy.optimize import fsolve


N = 10
K = 5
a_n = np.random.uniform(0, 1, (N, 1))
b_nk = np.random.uniform(0, 1, (N, K))
beta = -1
# make sure that a_n add up to 1
a_n = a_n / np.sum(a_n)


part_1_lis = []
part_2_lis = []
def solve_equation(x_k):
    for k in range(K):
        for n in range(N):
            part_1 = (b_nk[n, k] * x_k[k] **beta)
            part_2 = np.array([np.sum(np.array([b_nk[n, j] * x_k[j] 
                                       ** beta for j in range(K)])) + 1])

            part_3 = ((a_n[n] * (part_1/part_2)))
            part_1_lis.append(part_3)
            print(part_3)
            part_4 = np.sum(np.array([b_nk[n, j] * x_k[j] ** beta for j in range(K)]
                                    ))
            part_5 = (a_n[n] * (part_4 / part_2))
    #         print(part_5)
            part_2_lis.append(part_5)
    part_2_array = np.array(part_2_lis)
    part_1_array = np.array(part_1_lis)
    part_1_array = part_1_array.reshape(10,5)
    part_2_array = part_2_array.reshape(10,5)
    numerator = np.concatenate([part_1_array[i:i+2, :].reshape(-1, 1) for i in range(0, part_1_array.shape[0], 2)], axis=1)
    denominator = np.concatenate([part_2_array[i:i+2, :].reshape(-1, 1) for i in range(0, part_2_array.shape[0], 2)], axis=1)
    return x_k - (numerator / (1 - denominator))

        

#             return g

# Initial guess for x_k

x0 = np.array(np.arange(-6, -1 ), dtype = 'float32')


x_k = fsolve(func, x0)

print(x_k)