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?