I'm trying to find the quantile function of the two-term gaussian.
From https://statproofbook.github.io/P/norm-qf.html, I've got that I can take the inverse of the CDF of the two-term gaussian. I've got the two-term CDF as $${1\over2}\left[2+erf\left({(x-\mu_1)\over\sqrt2\sigma_1}\right)+erf\left({(x-\mu_2)\over\sqrt2\sigma_2}\right)\right]$$
Sadly, $erf\left({x-\mu_1\over\sqrt2\sigma_1}\right)+erf\left({x-\mu_2\over\sqrt2\sigma_2}\right)$ is not equal to $erf\left({x-\mu_1\over\sqrt2\sigma_1}+{x-\mu_2\over\sqrt2\sigma_2}\right)$. (I examined it numerically.)
I tried to make a sum function $$erf\_sum(x,\mu_1,\sigma_1,\mu_2,\sigma_2) = {2\over\sqrt\pi}\left[\int_{0}^{x-\mu_1\over\sqrt{2}\sigma_1}e^{-t_1^2}dt_1+\int_{0}^{x-\mu_2\over\sqrt{2}\sigma_2}e^{-t_2^2}dt_2\right]$$ with the hope I could do some change-of-variables trickery to bring the two integrals together into an $erf$-like form so I could still express the result using $erfinv$, but I wasn't able to get very far. The domains of the integrals look like they'll forever be different. So, I'm stuck here. Any help is appreciated.