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:

[![enter image description here][1]][1]

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][3] 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][4] (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!


  [1]: https://i.sstatic.net/Z7vma.png
  [2]: https://i.sstatic.net/vgaWW.png
  [3]: https://stackoverflow.com/questions/66880202/r-nleqslv-difficulties-solving-for-ph-in-an-acid-base-buffer
  [4]: https://chemistry.stackexchange.com/a/43422