Calculating simulatenous equilibrium concentrations given dissociation constants

Reactants A, B, and C form the complexes AB and BC according to:


$$ AB \rightleftharpoons A + B, \quad BC\rightleftharpoons B + C $$


Each reaction has a dissociation constant, $K\_{dA}$ and $K\_{dC}$, respectively, which are defined as:


$$ K\_{dA} = \frac{[A][B]}{[AB]}, \quad K\_{dC} = \frac{[C][B]}{[BC]} $$


Given starting concentrations $[A]\_0, [B]\_0, [C]\_0$, how does one calculate the equilibrium concentrations of the five species in mixture?


I have tried a couple of different approaches to a solution. First off, I identified the three additional equations


$$ [A] = [A]\_0- [AB] $$
$$ [B] = [B]\_0 - [AB] - [BC] $$
$$ [C] = [C]\_0- [BC] $$


Putting these additional expressions for [A], [B], and [C] into the equations for the dissociation constants ultimately yields two second-degree polynomials in [AB] and [BC] with a cross-term [AB][BC] which I have not figured out how to solve analytically.


$$ 0 = [AB]^2 - (K\_{dA} + [A]\_0 +[B]\_0)[AB] + [A]\_0[B]\_0 + [AB][BC] - [A]\_0[BC]$$
$$ 0 = [BC]^2 - (K\_{dB} + [B]\_0 +[C]\_0)[BC] + [B]\_0[C]\_0 + [AB][BC] - [C]\_0[AB]$$


The other attempt I made was to solve the system of five equations for $K\_{dA}$, $K\_{dC}$, [A], [B], and [C] numerically using root-finding algorithms. This works okay, but oftentimes the algorithm finds roots with negative concentrations which are obviously not desirable.


I have seen many solutions online for the classic case with one reaction where a so-called ICE-table is used, but I have found no extension of this for the case I have presented. Is there a way to solve this problem analytically?


