Another post has an excellent "exact" calculation for the [H+] (or pH) of a generic diprotic acid dissociating in water.
Combining acid dissociation constants to determine pH of diprotic acid
I am interested in what modifications would be needed to do a similar "exact" calculation with the addition of a strong base to effectively create a buffer solution.
When I look at the equations, my only thought is to modify equation (5) by adding in some strong base "counter ion" to the charge balance equation. For example, if NaOH was added to the diprotic acid, some of the acid would be neutralized by OH- to produce some initial conjugate base(s) in the mix. So there would be Na+ to account for in the charge balance. However, I don't think equation (4) would change, since the total amount of H2A, HA-, and A2- would not change -- just the ratios would change, which should fall out of the calculation.
I can't think of any other changes to make to the set of equations. I proceeded to solve the new equations (with major help from the Python Sympy forum). To test the equation, I assumed that the total "I" from equation 4 was 1.0 molar, and the initial concentration of NaOH and therefore the final concentration of Na+ was .5 molar. For Ka1 and Ka2 acid dissociation constants I used the values suggested, Ka1 = 5.9e-2, Ka2 = 6.4e-5, and for water Kw = 1.0e-14.
Using the nroots() function in Python, an equilibrium value for [H+] of 0.04872 was obtained.
If the diprotic acid is approximated as a monoprotic acid, I think the equilibrium should represent a "half-equivalence" point buffer, where ideally pH = pKa, and [H+] = Ka1. Comparing the [H+] of 0.049 with the Ka1 of .059 seems not close enough. I think typically the contribution of the second dissociation to [H+] is small and would only add to the approximation and increase the difference.
I'm not confident that my changes to the original set of equations is valid for the scenario. If anyone else is interested in these "exact solution" cases and this modification I'm interested in discussion. The equations are quite complicated, ending up with a quartic polynomial, but I can try to provide them. My primary interest is in what kind of modifications to the original formulation are valid for the "buffer" case.
For anyone interested, below is the quartic equation (to be set equal to zero) that includes $[Na^{+}{]}$ which essentially represents adding that much strong base, converting that much acid to the same amount of conjugate base, forming a buffer. I'm not clear yet on what ranges of initial conditions it is valid for (such as adding more NaOH beyond the final equivalence point).
$\displaystyle - K_{a1} K_{a2} K_{w} + K_{a1} K_{a2} [H^{+}{]} [Na^{+}{]} + K_{a1} \left([H^{+}{]}\right)^{3} + K_{a1} \left([H^{+}{]}\right)^{2} [Na^{+}{]} + \left([H^{+}{]}\right)^{4} + \left([H^{+}{]}\right)^{3} [Na^{+}{]} + \left([H^{+}{]}\right)^{2} \left(- I K_{a1} + K_{a1} K_{a2} - K_{w}\right) + [H^{+}{]} \left(- 2 I K_{a1} K_{a2} - K_{a1} K_{w}\right)$
(as mentioned in a comment below, I believe I erroneously assumed the Henderson-Hasselbach equation would be accurate for this case, but in fact $K_{a1}$ is sufficiently large so the small-x approximation causes too much error, so $[H^{+}{]} = K_{a1}$ at half-equivalence does not hold accurately.)
