There is an error of one order for concentration in Borek calculation (should be 0.05 and 0.1 approximately for concentration of H2PO4(-) and PO4(3-), respectively).
[PO4(3-)] + [HPO4(2-)] + [H2PO4(-)] + [H3PO4] = 0.15
I think it is consistent with the question posted. Correct me if I am wrong.
Moreover, from the information on program he used, the activity, not concentrations, were used in his calculations, and this may give substantial differences of calculations with using multicharged ions.
I haven't posted any information about the program I am using so you are just guessing. And you are wrong.
Results shown were for concentration calculations. Here are results for activities (substance formula, substance concentrations, activity coefficient, substance activity):
H3PO4 2.692e-013 * 1.000 = 2.692e-013 (p = 12.57)
H2PO4^- 3.614e-005 * 0.559 = 2.019e-005 (p = 4.69)
HPO4^2^- 1.386e-001 * 0.097 = 1.349e-002 (p = 1.87)
PO4^3^- 1.135e-002 * 0.005 = 6.005e-005 (p = 4.22)
H^+ 1.690e-010 * 0.559 = 9.442e-011 (p = 10.02)
OH^- 1.896e-004 * 0.559 = 1.059e-004 (p = 3.98)
NH4+ 1.156e-002 * 0.559 = 6.458e-003 (p = 2.19)
NH3 3.844e-002 * 1.000 = 3.844e-002 (p = 1.42)
K+ 2.999e-001 * 0.559 = 1.675e-001 (p = 0.78)
KOH 5.615e-005 * 1.000 = 5.615e-005 (p = 4.25)
Ionic strength is 0.968 so it is too high for a Debye-Huckel theory anyway.