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 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-}]$$

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^+]$ (called $x$ over there).

I would like to know which is the best way to simulate a titration curve.

And also: are the equations i used correct? do you know what am I missing?

Thanks!