I've been learning how to simulate a 1M phosphoric acid titration curve using numerical methods in R.
So far has this has been the best curve:
Notice how it flattens out (artificially) towards the left.
I am not certain if this is a programming and /or or a "chemistry" issue.
The system is defined by charge balance, mass balance, and equilibrium equations:
$$[H_3A] + [H_2A^-] + [HA^{2-}] + [A^{3-}] = P_{CA}$$ $$[H^+] + [Na^+] = [H_2A^-] +2[HA^{2-}] +3[A^{3-}] + Kw/[H^+]$$ $$Ka_1 = [H_2A^-][H^+]/[H_3A]$$ $$Ka_2 = [HA^{2-}][H^+]/[H_2A^-]$$ $$Ka_3 = [A^{3-}][H^+]/[HA^{2-}]$$$$[\ce{H3A}] + [\ce{H2A-}] + [\ce{HA^{2-}}] + [\ce{A^{3-}}] = P_{CA}$$
$$[\ce{H+}] + [\ce{Na+}] = [\ce{H2A-}] + 2[\ce{HA^{2-}}] + 3[\ce{A^{3-}}] + Kw/[\ce{H+}]$$
$$K_{a, 1} = [\ce{H2A-}] [\ce{H+}] / [\ce{H3A}]$$
$$K_{a, 2} = [\ce{HA^{2-}}] [\ce{H+}] / [\ce{H2A^-}]$$
$$K_{a, 3} = [\ce{A^{3-}}] [\ce{H+}] / [\ce{HA^{2-}}]$$
I've posted a very similar question at SO and shared code for the simulation over there, using R's nleqslv::nleqslv (go have a look if you'd like to).
I thought this is also a good place for this issue. For instance, this answer (and the ones it links to) do seem to be useful, but I'm not sure how to solve that equation for $[H_3O^+]$$[\ce{H3O+}]$ (called $x$ over there).
I would like to know which is the best way to simulate a titration curve.
And also: are the equations iI used correct? do you know what am I missing?
Thanks!
