How to calculate pH of a solution when H3PO4, NaH2PO4, Na2HPO4 and Na3PO4 are mixed together in certain amounts to form a solution?

I can provide an example which is useful:

$1,2,3$ and $4$ moles are taken of $ce{H3PO4}$, $ce{NaH2PO4}$, $ce{Na2HPO4}$ and $ce{Na3PO4}$ respectively with $K_{a.1}=10^{-3},K_{a.2}=10^{-7},K_{a.3}=10^{-13}$. Find the resulting $pH$ when these are mixed together to form an aqueous solution.

First write the $3$ ionic equations that are relevant here for the given $K_a$ values and then pair them up in such a way that you get rid of $ce{H^+}$ terms on both sides and the $K$ value is of a very high magnitude such that the reaction shift almost completely in the forward direction.

For eg. in the given question, we can pair up the ionic equation for $K_{a.1}$ and reversed equation for $K_{a.3}$ such that the net $K=10^{10}$ and then we put in the value of the moles of each substance:

$ce{H3PO4 + PO4^{3-}<=>H2PO4^{-} + HPO4^{2-}}$

 1        4            2          3
 _        3            3          4

Similarly we can pair $K_{a.2}$ and reversed $K_{a.3}$ ionic equations, and are finally left with $10$ mole of $ce{HPO4^{2-}}$ for which pH can be evaluated with the given data as $10$.

Usually doing calculations of this kind is not hard. Roughly, you start from some idea of where the pH is going to end up. For example, if you are only adding these phosphates, count the number of phosphate ions and the number of sodium ions. Divide the number of sodium ions by the number of phosphate ions. You (hopefully) get a number $x$ between $0$ and $3$, say $x=2.5$. Almost all phosphates will then be either $ce{NaH2PO4}$ or $ce{Na2HPO4}$. Pretend that these are the only ions and treat the problem like an ordinary two component buffer. Then you are done. This is an approximation but it will work rather well, as long as the $mathrm{p}K_mathrm{a}$ of the various acids are different enough. If this is not good enough, you just treat this as a bunch of acids, with a bunch of $mathrm{p}K_mathrm{a}$ and use the acid-base equilibrium equation over and over, that is $$frac{[ce{H}][ce{B}_i]}{[ce{A}_i]} = mathrm{p}K_mathrm{a}$$ where $[ce{A}_i]$ is the concentration of an acid "$i$" and $[ce{B}_i]$ that of the conjugate base. Given the pH, this gives you the ratio of the concentration of each base to each acid.

But, really, for phosphate buffers, you are very well off if you just use one acid base pair. For acids like terephthalic acid or isophthalic acid, which have two $mathrm{p}K_mathrm{a}$ values that are rather close to each other, that is they differ by less than $1.5$ or so, finding and solving such polynomial equations makes possible sense. Otherwise, other approximations we routinely make are more important.

For an explicit example, suppose that we consider a compound $ce{BH2}$ that can give up two $ce{H+}$ ions, making both $ce{BH-}$ and $ce{B^2-}$. Suppose that we have a total molar concentration of $ce{B}$, $$B_t=[ce{BH2}]+[ce{BH-}]+[ce{B^2-}],$$ and that we have also added strong bases or strong acids to change the $[ce{H+}],$ so that the total possible molar concentration of $[ce{H+}]$ is $$H_t=[ce{H^+}]+[ce{BH^-}]+2[ce{BH_2}].$$

We then also have the relationships begin{align} [ce{H+}][ce{BH-}] &= K_{mathrm{a},ce{BH_2}}[ce{BH_2}]\ [ce{H+}][ce{B^2-}] &= K_{mathrm{a},ce{BH-}}[ce{BH-}]\ end{align}

If we now call $H=[ce{H+}]$ and $B=[ce{B^2-}]$ we find the equations begin{align} H_t &= H + B (K_{mathrm{a},ce{BH-}}(H + 2K_{mathrm{a},ce{BH2}} H))tag{1},\ B_t &= B (1 + K_{mathrm{a},ce{BH^-}} H (1 + K_{mathrm{a},ce{BH2}} H))tag{2}.\ end{align}

We need to solve these equations for $B$ and $H$. Equation $(2)$ is easy to solve for $B$ and we can then substitute into equation $(1)$ to find $$H_t = H + B_t frac{K_{mathrm{a},ce{BH-}}(H + 2K_{mathrm{a},ce{BH_2}}H)} {(1 + K_{mathrm{a},ce{BH-}} H(1 + K_{mathrm{a},ce{BH_2}}H))}tag{3}.$$

We need to solve this formula for $H$, and (of course) $ce{pH}=-log_{10}(H)$. Some algebra from this point gets us to a cubic equation for $H$. Of course this has three possible solutions. But, only one is sensible, provided we insist on "physical sense", e.g. all the equilibrium constants and concentrations are positive real numbers. Somehow, this must be easy to see from the rules about roots of polynomials. I am not currently willing to do that proof.

Physically, however, the total amount of hydrogen bound to the base will increase with increasing (positive) $H$. So the right hand side of $(3)$ increases with increasing $H$ - from $H$ for $H=0$ to $2B_t + H$ as $Htoinfty$. As the right hand side of $(3)$ increases with increasing $H$ there is exactly one positive solution. However, this right hand side is also not a very "nice" function - the rational fraction is approximately constant except when we are near a "buffer" condition, where it changes quickly. So, probably an effective way to solve the equation would be to use a computer and "regula falsi" or something like that. A little thought will allow sensible guesses for any set of $H_t$ and $B_t$. Then regula falsi will converge quickly.
Or, from a good enough initial guess even Newton-Ralphson.

Or, possibly it would be better in many cases to figure out which acid/base pair is closest to being in the "buffer" region. And then, suppose that you have $H$ in that buffer region. Solve for $H$ with only that acid-base pair. Then, figure out what you did wrong - there can be $ce{H+}$'s that are not bound to the base in the buffer pair you have chosen or $ce{H+}$'s that are bound to stronger acids. Assume that $H$ calculated from a single buffer pair is correct. Calculate how many $ce{H+}$'s are unexpectedly bound or unbound from your assumed buffer pairs. Use this to calculate a new $H_t$. If $H_t$ has not changed much, which will be the case unless there is another acid-base pair with very nearly the same $mathrm{p}K_mathrm{a}$ - differing I suppose by less than an integer from that of the pair you are considering, then use this value of $H_t$ to solve again for $H$, and recurse.

For the stated problem, all the $mathrm{p}K_mathrm{a}$ are very different so the $K_mathrm{a}$ values differ by orders of magnitude. And, really, all you need do to is figure out which buffer region you are near. In principle, you can then find tiny errors using the recursion. But, this is silly in most circumstances. This will get a clear, continuous curve which always has the correct slope. But, it will vastly overstate what you know about the problem. The changes will be smaller than (I expect) lots of stuff: deviations from ideal solutions, errors in $K_mathrm{a}$ values, etc.. ).