I have a problem when solving a system of equations. To understand the code,

the term: (1. – Ni_oct + delta_Ni_mod) used to be a variable. But I have to implement it in order to remove one equation that is not seen here. This equation is what is seen in the variable Q. Q should equal to iterable j from the loop. The goal is to find for which delta_Ni_mod the diff satisfies the tolerance. But for some reason the code goes not find any solutions and is infinitely running. Increasing the tolerance is not an option, as then I get solutions that logically not satisfy a suitable difference. Btw I should have mentionen that values from the iterable j are very small numbers like between 10^-170 and 10^-21

Would be grateful if someone helps out. What I tried so far is below

import numpy as np
from scipy.optimize import root

# Initialize variables
delta_Ni_mod_max = np.float64(-0.002)
delta_Ni_mod_min = np.float64(0.0)
tolerance = 1e-21
diff_list = []
delta_Ni_mod_values = []
resultlist_def1_3_4_5_6 = []

for j, l, Kp_def6_subarray in zip(Kp_def3, Kp_def5_neutr, Kp_def6_neutr):
    for m in Kp_def6_subarray:
        def equations(vars):
            Ni_oct, V_oct, O_o, V_o, Al_tet, Al_oct, V_tet = vars

            eq1 = Ni_oct + V_oct + (1. - Ni_oct + delta_Ni_mod) + Al_oct - 2.0  # Oktaederplätze
            eq2 = O_o + V_o - 4.0  # Sauerstoffplätze
            eq3 = -1 * Ni_oct - 3 * V_oct + 2 * V_o + Al_tet - 2 * V_tet  # Ladungsbilanz
            eq7 = Al_tet + V_tet - 1  # Tetraederplätze
            eq8 = Al_tet + Al_oct - 2.0  # Massebilanz Al
            eq10 = (Al_oct * V_tet) / (V_oct * Al_tet) - l  # Kp_def5_neutr Site exchange
            eq11 = ((V_o * (Ni_oct ** 2)) / (O_o * ((1. - Ni_oct + delta_Ni_mod) ** 2))) - m  # K Bildung Sauerstoffvakanzen

            return [eq1, eq2, eq3, eq7, eq8, eq10, eq11]

        initial_guess = [0.999,  # Ni_oct # 0.999
                         0.0001,  # V_oct
                         3.99,  # O_o
                         0.0001,  # V_o
                         1.,  # Al_tet
                         1.,  # Al_oct
                         0.001]  # V_tet

        while True:
            Q = (V_oct * (1 - Ni_oct + delta_Ni_mod) ** 2) / ((Ni_oct) ** 3)
            diff = j - Q

            if abs(diff) <= tolerance:
                sol = root(equations, initial_guess, method='lm', tol=1e-40)
                diff_list.append(diff)
                delta_Ni_mod_values.append(delta_Ni_mod)
                resultlist_def1_3_4_5_6.append(sol.x)
                break

            if diff > 0:
                # Q muss positiver gemacht werden, also delta_Ni muss richtung negative Zahl gebtacht werden

                delta_Ni_mod_min = delta_Ni_mod
            else:
                # Hier gilt K < Q, also muss Q negativer gemacht werden
                # Q muss negativer gemacht werden, also delta_Ni Ri 0 gebracht werden
                delta_Ni_mod_max = delta_Ni_mod

            delta_Ni_mod = (delta_Ni_mod_min + delta_Ni_mod_max) / 2```