How to increase the precision of scipy.optimize.fsolve function

Question

I want to model PCR reactions, so I have a bunch of algebraic non-linear equations that I and solving using the fsolve function. But after 23 cycles the float overflows and the results come out incorrect. By incorrect i mean that the values of the variables (x[0], x[1], ..) become negative.

from scipy.optimize import fsolve, least_squares
import numpy as np

cycles = 35

KH1 = 5.5531*(10**-13)
KH2 = 8.1493*(10**-11)
KP = 10**-2

P1 = [None] * (cycles+1)
P2 = [None] * (cycles+1)
T1 = [None] * (cycles+1)
T2 = [None] * (cycles+1)
H1 = [None] * (cycles+1)
H2 = [None] * (cycles+1)


P1[0] = 10**-6
P2[0] = 10**-6
T1[0] = 10**-11
T2[0] = 10**-11

p1 = P1[0]
p2 = P2[0]
t1 = T1[0]
t2 = T2[0]


def func(x):
    e1 = (x[4] * KH1) - (x[0] * x[2])
    e2 = (x[5] * KH2) - (x[1] * x[3])
    e3 = p1 - (x[0] + x[4] + x[6])
    e4 = p2 - (x[1] + x[5] + x[6])
    e5 = t1 - (x[2] + x[4])
    e6 = t2 - (x[3] + x[5])
    e8 = (x[6] * KP) - (x[0] * x[1])
    return [e1,e2,e3,e4,e5,e6,e8]

for i in range(cycles):
    sol = fsolve(func, [0.001,0.001,0.001,0.001,0.001,0.001,0.001])

    H1[i] = sol[4]
    H2[i] = sol[5]

    fH1 = H1[i] / T1[i]
    fH2 = H2[i] / T2[i]
    
    efficiency = 0.5 * (fH1 + fH2)
    
    print(i+1)
    print(efficiency)
    
    P1[i+1] = P1[i] - H1[i]
    P2[i+1] = P2[i] - H2[i]
    T1[i+1] = T1[i] + H2[i]
    T2[i+1] = T2[i] + H1[i]

    p1 = P1[i+1]
    p2 = P2[i+1]
    t1 = T1[i+1]
    t2 = T2[i+1]

Is there a way to increase the precision of the results for at least 35 cycles.

The output for 35 cycles:

Cycle: 1; Efficiency: 0.9999589746574444
Pri 1: 9.999900000055537e-07
Pri 2: 9.999900008149532e-07
Temp 1: 1.9999185046783575e-11
Temp 2: 1.999999444636531e-11
Cycle: 2; Efficiency: 0.9999589738335422
Pri 1: 9.999700008316142e-07
Pri 2: 9.999700024504453e-07
Temp 1: 3.999754955458496e-11
Temp 2: 3.999916838595971e-11
Cycle: 3; Efficiency: 0.9999589721929882
Pri 1: 9.999300033042744e-07
Pri 2: 9.99930006542e-07
Temp 1: 7.999345799992356e-11
Temp 2: 7.999669572574256e-11
Cycle: 4; Efficiency: 0.9999589689112148
Pri 1: 9.998500098907068e-07
Pri 2: 9.998500163665432e-07
Temp 1: 1.5998363345666242e-10
Temp 2: 1.599901092934185e-10
Cycle: 5; Efficiency: 0.9999589623463345
Pri 1: 9.99690026346127e-07
Pri 2: 9.996900392996067e-07
Temp 1: 3.1996070039317954e-10
Temp 2: 3.1997365387318973e-10
Cycle: 6; Efficiency: 0.9999589492105335
Pri 1: 9.99370065823541e-07
Pri 2: 9.993700917382614e-07
Temp 1: 6.399082617384076e-10
Temp 2: 6.399341764593294e-10
Cycle: 7; Efficiency: 0.9999589229142298
Pri 1: 9.987301579176373e-07
Pri 2: 9.987302097792151e-07
Temp 1: 1.2797902207846395e-09
Temp 2: 1.2798420823631377e-09
Cycle: 8; Efficiency: 0.9999588702215008
Pri 1: 9.974503684094203e-07
Pri 2: 9.974504722634983e-07
Temp 1: 2.559527736501506e-09
Temp 2: 2.5596315905802006e-09
Cycle: 9; Efficiency: 0.9999587644316008
Pri 1: 9.948908421016906e-07
Pri 2: 9.948910503398154e-07
Temp 1: 5.118949660184387e-09
Temp 2: 5.119157898309939e-09
Cycle: 10; Efficiency: 0.9999585512150795
Pri 1: 9.897718953137677e-07
Pri 2: 9.897723139348765e-07
Temp 1: 1.0237686065123224e-08
Temp 2: 1.023810468623287e-08
Cycle: 11; Efficiency: 0.9999581181067771
Pri 1: 9.7953421505308e-07
Pri 2: 9.795350610263854e-07
Temp 1: 2.0474938973614262e-08
Temp 2: 2.0475784946920725e-08
Cycle: 12; Efficiency: 0.9999572241744039
Pri 1: 9.590592879358986e-07
Pri 2: 9.590610159597535e-07
Temp 1: 4.094898404024612e-08
Temp 2: 4.095071206410206e-08
Cycle: 13; Efficiency: 0.9999553167605227
Pri 1: 9.181103286655014e-07
Pri 2: 9.181139387457051e-07
Temp 1: 8.18960612542945e-08
Temp 2: 8.189967133449928e-08
Cycle: 14; Efficiency: 0.9999509417304484
Pri 1: 8.362143218009405e-07
Pri 2: 8.362222487313798e-07
Temp 1: 1.6378775126861994e-07
Temp 2: 1.6379567819906017e-07
Cycle: 15; Efficiency: 0.9999389959080549
Pri 1: 6.724267058021067e-07
Pri 2: 6.724464196692128e-07
Temp 1: 3.2756358033078687e-07
Temp 2: 3.275832941978939e-07
Cycle: 16; Efficiency: 0.999881092518708
Pri 1: 3.448636529415585e-07
Pri 2: 3.449405021781882e-07
Temp 1: 6.550694978218116e-07
Temp 2: 6.551463470584422e-07
Cycle: 17; Efficiency: 0.5264117174519147
Pri 1: 6.173497903170838e-13
Pri 2: 9.056777506975753e-11
Temp 1: 9.9991943222493e-07
Temp 2: 1.0000093826502103e-06
Cycle: 18; Efficiency: 4.55884719669062e-05
Pri 1: 3.4283980122736808e-19
Pri 2: 7.38063736992815e-15
Temp 1: 1.0000099926193623e-06
Temp 2: 1.0000099999996578e-06
Cycle: 19; Efficiency: 3.6901525033426808e-09
Pri 1: -7.56944722967757e-24
Pri 2: 6.014075684520175e-19
Temp 1: 1.0000099999993982e-06
Temp 2: 1.0000100000000007e-06
Cycle: 20; Efficiency: 3.006739577729103e-13
Pri 1: -4.26443257940018e-24
Pri 2: 5.0334412390713845e-23
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 21; Efficiency: 2.4499055249910828e-17
Pri 1: -4.261015212659181e-24
Pri 2: 1.3323945430461554e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 22; Efficiency: 3.537919550216741e-21
Pri 1: -4.264358355576793e-24
Pri 2: 1.328661776104943e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 23; Efficiency: 6.586578432833721e-24
Pri 1: -4.264651042496745e-24
Pri 2: 1.3289412897362976e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 24; Efficiency: -4.669280931707777e-22
Pri 1: -4.264322421569872e-24
Pri 2: 1.3295465343343283e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 25; Efficiency: -4.033381628412326e-18
Pri 1: -9.643336312842848e-25
Pri 2: 6.096401668505964e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 26; Efficiency: -1.7928462660600777e-22
Pri 1: -9.645834663351964e-25
Pri 2: 6.0970100763957804e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 27; Efficiency: 9.465002751917195e-18
Pri 1: 3.1143527079562413e-24
Pri 2: -1.6912120901785106e-23
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 28; Efficiency: -5.4327442351155336e-18
Pri 1: -4.2573287015899444e-24
Pri 2: 1.325157632876866e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 29; Efficiency: 2.220455013323325e-21
Pri 1: -4.264413609414534e-24
Pri 2: 1.3278015862657087e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 30; Efficiency: -9.181589730610556e-22
Pri 1: -4.264322788059294e-24
Pri 2: 1.3295471012197703e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 31; Efficiency: 5.4326757292416825e-18
Pri 1: 3.1152266062187542e-24
Pri 2: -1.691546240505624e-23
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 32; Efficiency: -9.463688877357565e-18
Pri 1: -9.63949706848734e-25
Pri 2: 6.091280936503942e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 33; Efficiency: 9.461116532403566e-18
Pri 1: 3.116935792777572e-24
Pri 2: -1.6912026850260163e-23
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 34; Efficiency: -5.428891311058792e-18
Pri 1: -4.2643615966129534e-24
Pri 2: 1.3271617390741832e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06
Cycle: 35; Efficiency: 5.429911312197918e-18
Pri 1: -4.2643615966129534e-24
Pri 2: 1.3271617390741832e-24
Temp 1: 1.0000099999999996e-06
Temp 2: 1.0000100000000007e-06

Answer 1

Simply use x**2 instead of x:

...

def opt_func(x):
    """ Function for optimization. """
    # Converts all negative `x` to positive
    # `longdouble` might increase math operations precision
    return func(np.power(x, 2, dtype=np.longdouble))


for i in range(cycles):

    x_init = np.random.random(size = 7) * 0.001
    sol = fsolve(opt_func, x_init)

    sol = np.power(sol, 2, dtype=np.longdouble)
    print(i+1)

    # check that all values are positive
    print(sol)

    # check that solution close to zero
    print(func(sol))
    ...