How to (numerically) model a phosphoric acid titration curve
Much easier is to calculate the inverse function $[\ce{Na+}]=f([\ce{H+}],K_\mathrm{a1},K_\mathrm{a2},K_\mathrm{a3})$.
Calculate fractions of respective phosphate forms as the function of $\ce{[H+]}$ and acidity constants.
First, calculate the common denominator $CD$:
$$a_1 = [\ce{H+}]^3$$
$$a_2 = [\ce{H+}]^2 \cdot K_\mathrm{a1} = a_1 \cdot \frac {K_\mathrm{a1}}{[\ce{H+}]}$$
$$a_3 = [\ce{H+}] \cdot K_\mathrm{a1} \cdot K_\mathrm{a2} = a_2 \cdot \frac {K_\mathrm{a2}}{[\ce{H+}]}$$
$$a_4 =  K_\mathrm{a1} \cdot K_\mathrm{a2} \cdot K_\mathrm{a3} = a_3 \cdot \frac {K_\mathrm{a3}}{[\ce{H+}]}$$
$$CD = a_1 + a_2 + a_3 + a_4$$
.. and then reuse it and its additive terms(*):
$$[\ce{H3A}]=c_0 \cdot \frac {a_1}{CD}$$
$$[\ce{H2A-}]=c_0 \cdot \frac {a_2}{CD}$$
$$[\ce{HA^2-}]=c_0 \cdot \frac {a_3}{CD}$$
$$[\ce{A^3-}]=c_0 \cdot \frac {a_4}{CD}$$
$c_0$ is the total phosphate concentration. From the fractions, you can get $\ce{[Na+]}$:
$$\ce{[Na+]} = \ce{[H2A-]} + 2 \cdot \ce{[HA^2-]} + 3 \cdot \ce{[A^3-]} + \frac{K_\mathrm{w}}{[\ce{H+}]} - \ce{[H+]}$$

There is one important thing to consider: Titration dilutes the acid.
If  there is the volume $V_0$ of $\ce{H3PO4}$ with molar concentration $c_\mathrm{0,init}$, and if we spend the volume $V_1$ of $\ce{NaOH}$ solution of concentration $c_1$, then:
$$[\ce{Na+}]=c_1 \cdot \frac {V_1}{V_0 + V_1 }$$
$$c_0=c_\mathrm{0,init} \cdot \frac {V_0}{V_0 + V_1}$$
So perhaps the best approach can be:

There are given concentrations of respective solutions
For given volumes of respective solutions calculate $c_0$ and $\ce{[Na+]}$
For given $c_0$ and $\ce{[Na+]}$, calculate $\ce{[H+]}$ by finding the solution of the equation below ( $f([\ce{H+}]) =  0$ ) by numerical methods, like Regular falsi, or Newton one.

$$ \ce{[H2A-]} + 2 \cdot \ce{[HA^2-]} + 3 \cdot \ce{[A^3-]} + \frac{K_\mathrm{w}}{[\ce{H+}]} - \ce{[H+]} - \ce{[Na+]} = 0$$
For $y_1=f(x_1)$, $y_2=f(x_2)$, $y_1 \cdot y_2 \lt 0$, there is the Regula falsi iteration:
$$x_3 =  \frac{x_2 \cdot y_1 - x_1 \cdot y_2  }{  y_1 - y_2}$$

(*) Note that the formulas are quite common, I have just rearrange them. They should be able to be found on internet. I have once derived them myself, but not this time, because I remember them.