Calculate pH of a weak acid and weak base neutralization reaction
You have a system of an acid, an ampholyte and a base. As you are interested in the pH value, your main question is what does "produce" or "consume" them. The proton concentration for your specific system is (citric Acid, sodium Bicarbonate, sodium Carbonate):
$$\ce{[H+] = [OH-] + \underbrace{[H2A-] + 2[HA^2-] + 3[A^3-]}_{citric acid's bases} \underbrace{+ [B^2-] - [H2B]}_{bicarbonate's deprotonated\\ and protonated form} \underbrace{- [CH+] - 2 [CH2^2+]}_{carbonate's sinlgy and\\ doubly protonated form}}$$

The multipliers arise from the fact, that a n-times deprotonated/protonated ion yields/consumes n protons.

Now the question is, where you get these concentrations from. It is by solving a linear equation system for each subsystem.



Citric acid, to be deprotonated three times:
\begin{align}
\ce{[H3A] + [H2O] &amp;&lt;=&gt; [H2A-] + [H3O+]\\
[H2A-] + [H2O] &amp;&lt;=&gt; [HA^2-] + [H3O+]\\
[HA^2-] + [H2O] &amp;&lt;=&gt; [A^3-] + [H3O+]}
\end{align}
This gives you the following linear equation system
\begin{align}
k_{a1} &amp;= \ce{\frac{[H2A-][H3O+]}{[H3A]}} \\
k_{a1} \ce{[H3A]} &amp;- \ce{[H2A-][H3O+]} = 0\\
k_{a2} &amp;= \ce{\frac{[HA^2-][H3O+]}{[H2A-]}} \\
k_{a2} \ce{[H2A-]} &amp;- \ce{[HA^2-][H3O+]} = 0\\
k_{a3} &amp;= \ce{\frac{[A^3-][H3O+]}{[HA-]}} \\
k_{a3} \ce{[HA-]} &amp;- \ce{[A^3-][H3O+]} = 0\\
\ce{[H3A] + [H2A-] &amp;+ [HA2-] + [A3-] = [H3A]0} = c_1\\
\end{align}
Which solves to:
\begin{align}
\ce{[H3A]} &amp;= \frac{c_1 x^3}{k_{a1} k_{a2} k_{a3} + k_{a1} k_{a2} x + k_{a1} x^2 + x^3} \\
\ce{[H2A-]} &amp;= \frac{c_1 k_{a1} x^2}{k_{a1} k_{a2} k_{a3} + k_{a1} k_{a2} x + k_{a1} x^2 + x^3} \\
\ce{[HA^2-]} &amp;= \frac{c_1 k_{a1} k_{a2} x}{k_{a1} k_{a2} k_{a3} + k_{a1} k_{a2} x + k_{a1} x^2 + x^3} \\
\ce{[A^3-]} &amp;= \frac{c_1 k_{a1} k_{a2} k_{a3}}{k_{a1} k_{a2} k_{a3} + k_{a1} k_{a2} x + k_{a1} x^2 + x^3}
\end{align}



Sodium bicarbonate, to be protonated or deprotonated:
$$
\ce{[HB^-] + [H2O] &lt;=&gt; [B^2-] + [H3O+]\\
[HB^-] + [H2O] &lt;=&gt; [H2B] + [OH-]}
$$
\begin{align}
k_{a4} &amp;= \ce{\frac{[B^2-][H3O+]}{[HB-]}} \\
k_{a4} \ce{[HB-]} &amp;- \ce{[B^2-][H3O+]} = 0\\
k_{b1} &amp;= \ce{\frac{[H2B][OH-]}{[HB-]}} \\
k_{b1} \ce{[HB-]} &amp;- \ce{[H2B][OH-]} = 0\\
\ce{[H2B] &amp;+ \ce{[HB-] + [B^2-] = [B]0}} = c_2\\
\end{align}
Which solves to:
\begin{align}
\ce{[H2B]} &amp;= \frac{c_2 x k_{b1}}{x k_{b1} + x k_W/x + k_{a4} k_W/x}
       = \frac{c_2 x k_{b1}}{x k_{b1} + k_W + k_{a4} k_W/x} \\
\ce{[HB-]} &amp;= \frac{c_2 x k_W/x}{x k_{b1} + x k_W/x + k_{a4} k_W/x}
       = \frac{c_2 k_W}{x k_{b1} + k_W + k_{a4} k_W/x} \\
\ce{[B^2-]} &amp;= \frac{c_2 k_{a2} k_W/x}{x k_{b1} + x k_W/x + k_{a4} k_W/x}
        = \frac{c_2 k_{a2} k_W/x}{x k_{b1} + k_W + k_{a4} k_W/x}\\
\end{align}



Sodium carbonate, to be doubly protonated:
$$
\ce{[C] + [H2O] &lt;=&gt; [CH+] + [OH-]\\
[CH+] + [H2O] &lt;=&gt; [CH2^2+] + [OH-]}
$$
\begin{align}
k_{b2} &amp;= \ce{\frac{[CH+][OH-]}{[C]}} \\
k_{b2} \ce{[C]} &amp;- \ce{[CH+][OH-]} = 0\\
k_{b3} &amp;= \ce{\frac{[CH2^2+][OH-]}{[CH+]}} \\
k_{b3} \ce{[CH+]} &amp;- \ce{[CH2^2+][OH-]} = 0\\
\ce{[C] &amp;+ [CH+] + [CH2^2+] = [C]0} = c_3\\
\end{align}
Which solves to:
\begin{align}
\ce{[B]} &amp;= \frac{c_3 (k_W/x)^2}{k_{b2} k_{b3} + k_{b2} k_W/x + (k_W/x)^2} \\
\ce{[BH+]} &amp;= \frac{c_3 k_{b2} k_W/x}{k_{b2} k_{b3} + k_{b2} k_W/x + (k_W/x)^2} \\
\ce{[BH^2+]} &amp;= \frac{c_3 k_{b2} k_{b3}}{k_{b2} k_{b3} + k_{b2} k_W/x + (k_W/x)^2}\\
\end{align}



Putting it all together (it looks a bit different, though, as it is copied from a simplified form of Mathematica):
$$x=\frac{k_W}{x} +
c_1 \frac{3 k_{a1} k_{a2} k_{a3}+2 x k_{a1} k_{a2}+x^2 k_{a1}}{k_{a1} k_{a2} k_{a3}+x k_{a1} k_{a2}+x^2 k_{a1}+x^3} +
c_2 \left(\frac{k_{a4} k_W}{k_{a4} k_W+x^2 k_{b1}+x k_W}-\frac{x^2 k_{b1}}{k_{a4} k_W+x^2 k_{b1}+x k_W}\right) -
c_3 \frac{2 x^2 k_{b2} k_{b3}+k_W x k_{b2}}{x^2 k_{b2} k_{b3}+k_W x k_{b2}+k_w^2}$$

Solving this little monster for $x$ and calculating its negative decadic logarithm gives you the wanted pH of your system, which is 5.65.



As you mention, further thinking about the problem makes it possible to simplify it a litle bit.
\begin{align}
\ce{[H+]  &amp;= [OH-] + [H2A-] + 2[HA^2-] + 3[A^3-] + [B^2-] - [H2B] - [CH+] - 2 [CH2^2+]} \\
\ce{[H+] &amp;= [OH-] + [H2A-] + 2[HA^2-] + 3[A^3-] + [CO3^2-]_B - [H2CO3]_B - [HCO3-]_C - 2 [H2CO3]_C} \\
x &amp;= \frac{k_W}{x} + c_1 \frac{x^2 k_{a1} + 2 x k_{a1} k_{a2} + 3 k_{a1} k_{a2} k_{a3}}{x^3 + x^2 k_{a1} + x k_{a1} k_{a2} + k_{a1} k_{a2} k_{a3}} + \frac{-(c_2 + 2 c_3) x^2 - c_3 x k_{a4} + 2 c_2 k_{a4} k_{a5}}{x^2 + x k_{a4} + k_{a4} k_{a5}} \\
\end{align}

This does not change the resulting pH value and thus it is not more accurate in my eyes, it only simplifies (admittedly quite a bit) the equation to be solved.

If we assume your initial pH value to be exactly 6 and that is caused by CO2, then we can add it to the equation.
$$x = \frac{k_W}{x} + c_1 \frac{x^2 k_{a1} + 2 x k_{a1} k_{a2} + 3 k_{a1} k_{a2} k_{a3}}{x^3 + x^2 k_{a1} + x k_{a1} k_{a2} + k_{a1} k_{a2} k_{a3}} + \frac{-(c_2 + 2 c_3) x^2 + (c_4 - c_3) x k_{a4} + 2 (c_2 + c_4) k_{a4} k_{a5}}{x^2 + x k_{a4} + k_{a4} k_{a5}} $$

Where $c_4$ can be calculated by rearranging the proton concentration equation for solely carbonic acid:
\begin{align}
x &amp;= \frac{k_W}{x} + c_4 \frac{x k_{a4} + 2 k_{a4} k_{a5}}{x^2 + x k_{a4} + k_{a4} k_{a5}} \\
c_4 &amp;= \left(x - \frac{k_W}{x}\right) \left(\frac{x^2 + x k_{a4} + k_{a4} k_{a5}}{x k_{a4} + 2 k_{a4} k_{a5}}\right) \\
c_4 &amp;= 3.21119 \cdot 10^{-6}
\end{align}

Without initial pH 6: 5.6467
With initial pH 6: 5.6465