Htykuut

How to calculate the concentration of all species present in a solution with $mathrm{0.3~M}$ $mathrm{NaH_2PO_4}$?

The whole truth is that any time that you add any phosphate ion into an aqueous solution, then you will have all four phosphate species. The $mathrm{NaH_2PO_4}$ solution contains the ions $mathrm{Na^+}$, $mathrm{H^+}$, $mathrm{OH^−}$, $mathrm{H_2PO_4^{-}}$, $mathrm{HPO_4^{2-}}$, $mathrm{PO_4^{3-}}$ and undissociated acid $mathrm{H_3PO_4}$.

I form seven independent equations in order to specify the unknown concentrations of the seven species present in an aqueous solution of $mathrm{NaH_2PO_4}$. these equations are:

$$
K_1 = frac{[mathrm{H^+}][mathrm{H_2PO_4^-}]}{[mathrm{H_3PO_4}]} = 7.6 times 10^{−3} tag{1}$$
$$K_2 = frac{[mathrm{H^+}][mathrm{HPO_4^{2-}}]}{[mathrm{H_2PO_4^-}]} = 6.2 times 10^{−8} tag{2}$$
$$K_{3} = frac{[mathrm{H^+}][mathrm{PO_4^{3-}}]}{[mathrm{HPO_4^{2-}}]} = 2.1 times 10^{−13} tag{3}$$
$$K_mathrm{w} = {[mathrm{H^+}][mathrm{OH^-}]} = 1 times 10^{−14} label{eq:4}tag{4}
$$
$$C_mathrm{P} =[mathrm{Na^+}]= [mathrm{NaH_2PO_4}]_0 ={0.3}tag{5}$$
$$C_mathrm{P} = [mathrm{PO_4^{3-}}] + [mathrm{HPO_4^{2-}}] + [mathrm{H_2PO_4^-}] + [mathrm{H_3PO_4}] tag{6}$$
$${[mathrm{H_3PO_4}] + [mathrm{H^+}] = [mathrm{HPO_4^{2-}}] + 2[mathrm{PO_4^{3-}}] +frac{K_mathrm{w}}{[mathrm{H^+}]}}tag{7}$$

Neglecting $[mathrm{Na^+}]$ and equation $eqref{eq:4}$, then given $C_mathrm{P} = {0.3}$ and $K_mathrm{w}=1times10^{−14}$, how to solve five equations to find five variables using wolframalpha for solving quintic equation?

Too long for a comment. This just deals with the numerical aspects of the problem.

Starting from @Batominovski answer, let us define both sides of the charge balance.
$$text{lhs}=y^5+(C_P+K_1),y^4+K_1K_2,y^3+K_1K_2K_3,y^2$$
$$text{rhs}=K_w,y^3+(K_1K_w+C_PK_1K_2),y^2+(K_1K_2K_w+2C_PK_1K_2K_3),y+K_1K_2K_3K_w$$
Instead of solving $text{lhs}=text{rhs}$ for $y$, let $y=10^{-text{pH}}$ and consider that we look for the zero of function
$$f(text{pH})=log(text{lhs})-log(text{rhs})$$ which, if plooted, looks quite linear.

As a proof, using a Taylor expansion built at $text{pH}=7$ (corresponding to $ y=sqrt{k_w}$), the first estimate is $text{pH}=4.65444$ while the exact solution is $text{pH}=4.66886$. One more iteration of Newton method and we are done (as shown below for a ridiculous number of significant figures)
$$left(
begin{array}{cc}
n & text{pH}_n \
0 & 7.00000000000 \
1 & 4.65444138638 \
2 & 4.66886040691 \
3 & 4.66886042380
end{array}
right)$$

For simplicity, write $x_0:=[texttt{H}_3texttt{PO}_4]$, $x_1:=[texttt{H}_2texttt{PO}_4^-]$, $x_2:=[texttt{H}texttt{PO}_4^{2-}]$, and $x_3:=[texttt{PO}_4^{3-}]$ so that $x_i$ is the concentration of the phosphate species with negative-$i$ charge, for $i=0,1,2,3$. Let $y$ denote $[texttt{H}^+]$. Therefore,
$$x_0+x_1+x_2+x_3=C_P,,$$
$$x_0+y=x_2+2x_3+frac{K_w}{y},,$$
$$K_1=frac{yx_1}{x_0},,$$
$$K_2=frac{yx_2}{x_1},,$$
and
$$K_3=frac{yx_3}{x_2},.$$
This gives $x_1=dfrac{K_1x_0}{y}$, $x_2=dfrac{K_2x_1}{y}=dfrac{K_1K_2x_0}{y^2}$, and $x_3=dfrac{K_3x_2}{y}=dfrac{K_1K_2K_3x_0}{y^3}$. Hence, we end up with two equations in $x_0$ and $y$:
$$x_0+frac{K_1x_0}{y}+frac{K_1K_2x_0}{y^2}+frac{K_1K_2K_3x_0}{y^3}=C_Ptag{*}$$
and
$$x_0+y=frac{K_1K_2x_0}{y^2}+frac{2K_1K_2K_3x_0}{y^3}+frac{K_w}{y},.tag{#}$$
From (*), we easily get
$$x_0=frac{C_P}{1+frac{K_1}{y}+frac{K_1K_2}{y^2}+frac{K_1K_2K_3}{y^3}},.$$
Plugging this into (#) yields
$$frac{C_P}{1+frac{K_1}{y}+frac{K_1K_2}{y^2}+frac{K_1K_2K_3}{y^3}},left(1-frac{K_1K_2}{y^2}-frac{2K_1K_2K_3}{y^3}right)+y-frac{K_w}{y}=0,.$$
That is,
$$y^5+(C_P+K_1)y^4+(K_1K_2-K_w)y^3-(K_1K_w+C_PK_1K_2-K_1K_2K_3)y^2-(K_1K_2K_w+2C_PK_1K_2K_3)y-K_1K_2K_3K_w=0,.$$
I think you need a numerical solver to solve for $y$ from the last equation. Once you know $y$, you can find the values of $x_0$, $x_1$, $x_2$, and $x_3$.

Mathematica tells me that there is only one positive real solution $y$ (or using Descartes's Rule of Signs, you know that there exists only one positive real solution) to the previous equation, which is $yapprox 2.1times 10^{-5}$. This means $x_0approx 8.4times 10^{-4}$, $x_1approx 0.30$, $x_2approx 8.6times 10^{-4}$, and $x_3approx 8.5times 10^{-12}$. In a sense, you should expect $x_1approx C_P=0.30$. After all, the dissociative reactions into other phosphate species are quite weak.