def fullsol(x, DIC, QCharge):

return [x[3] - (DIC/((1 +(K1/x[0]) + (K1*K2/(x[0])**2)))),
        x[1] - (DIC/((1 +(x[0]/K1) + (K2/(x[0]))))),
        x[2] - (DIC/(1 +(x[0]/K2) + ((x[0]**2)/K1*K2))),
        (x[0]*x[4]) - Kw,
        2*QCharge - x[0] + x[1] + 2*x[2] + x[4] - (DIC/0.78)]

x0 = ([0.1, 0.1, 0.1, 0.1, 0.1])

solution = optimize.least_squares(fullsol, x0, args=(0.0035, 0),bounds =(0,np.inf))