Solving system of nonlinear equations numerically in python

I want to solve the following problem numerically and first defined it as a function:

def F(x, na=1, nb=2, L=6, A=0.09, V=0.54, dh=0.03, k=10**-2, dB=0.005, fB=38, ur=0.39, rhoW=998, rhoG=1.2, g=9.81, pu=10**5, myW = 10**-3):
    
    #variables
    uWa = x[0]
    uWb = x[1]
    ta = x[2]
    tb = x[3]
    ea = x[4]
    eb = x[5]
    rhoa = x[6]
    rhob = x[7]
    ua = x[8]
    ub = x[9]
    pa0 = x[10]
    pb0 = x[11]
    Rea = x[12]
    Reb = x[13]
    lambdaa = x[14]
    lambdab = x[15]

    #nonlinear equations in form f=0:
    f = np.zeros(16)
    f[0] = ta - L/(ur+uWa)
    f[1] = tb - L/(ur-uWb)
    f[2] = ea - na*ta*fB*np.pi*dB**3/(6*V)
    f[3] = eb - nb*tb*fB*np.pi*dB**3/(6*V)
    f[4] = rhoa - rhoW + ea*(rhoW-rhoG)
    f[5] = rhob - rhoW + eb*(rhoW-rhoG)
    f[6] = rhoa*ua - ea*rhoG*(ur+uWa) - (1-ea)*rhoW*uWa
    f[7] = rhob*ub + eb*rhoG*(ur-uWb) - (1-eb)*rhoW*uWb
    f[8] = (-rhoa*ua + rhob*ub + rhoG*(na+nb)*fB*np.pi*dB**3/(6*A**2))*A
    f[9] = pb0 - pa0 - rhoa*rhob*(ua**2-ub**2)/(rhoa+rhob)
    f[10] = pa0 - pu - rhoa*((ua**2*lambdaa*dh*np.pi)/(8*A)+g)*L    
    f[11] = pb0 - pu + rhob*((ub**2*lambdab*dh*np.pi)/(8*A)+g)*L
    f[12] = Rea - rhoa*ua*dh/myW
    f[13] = Reb - rhob*ub*dh/myW  
    f[14] = lambdaa - fixpunktiteration(colebrook, 0.05, args=(Rea,k,dh)) 
    f[15] = lambdab - fixpunktiteration(colebrook, 0.05, args=(Reb,k,dh))

    return f

Following functions I defined for the equations f[14] and f[15]:

def colebrook(lam,args):
    Re, k, d = args
    return ( 1/(4*np.log10( 2.51/(Re*lam**0.5) + (k/d)/3.71 )**2) )

def fixpunktiteration(f,xalt,args,kmax=500,tol=1.e-7):
    for k in range(1,kmax):
        xneu = f(xalt,args)
        xdiff = xneu - xalt
        if abs(xdiff/xneu) < tol:
            break
        xalt = xneu
    else:
        xneu = None
    return xneu

Then I defined all the parameters and initial values of the function F and tried to solve with fsolve (after importing it from scipy):

xinit = np.array([uWa_0, uWb_0, ta_0, tb_0, ea_0, eb_0, rhoa_0, rhob_0, ua_0, ub_0, pa0_0, pb0_0, Rea_0, Reb_0, lambdaa_0, lambdab_0]) 
loesung = fsolve(func=F,x0=xinit)

It works without error but it returns all initial values as solution which I don't understand. If I substitute f[14]=64/Rea the program gives a solution which is different from the initial values but physical nonsense. So I think the problem must have something to do with equations f[14] and f[15] or colebrook but I don't find it.