You could look into the skew-normal distribution (see wikipedia, estimation for skew normal) and you could use it in the same way you used the normal distribution.
But, lacking any knowledge of how the $(x_i, y_i)$ pairs were obtained, there is no principled statistical way of estimating parameters. It doesn't look like you have IID data! So this is probably more a problem of function approximation, more numerical analysis than statistics (unless you tell us some context).
First, the two-piece normal$^*$ in Mudholkar and Huston (2000) is not the same as the skew-normal in Azzalini (1985). The density function of the former is given by:
$$g(x;\mu,\sigma,\epsilon) = \phi\left(\frac{x-\mu}{\sigma(1+\epsilon)}\right)I(x<\mu) + \phi\left(\frac{x-\mu}{\sigma(1-\epsilon)}\right)I(x\geq \mu),$$
$\epsilon \in (-1,1)$; while the density function of the latter is given by
$$g(x;\mu,\sigma,\lambda) = \frac{2}{\sigma}\phi\left(\dfrac{x-\mu}{\sigma}\right)\Phi\left(\lambda\dfrac{x-\mu}{\sigma}\right),$$
$\lambda\in(-\infty,\infty)$.
The estimation of the two-piece normal based on three quantiles is very popular in finance, where they construct "Fan plots" using three quantiles reported by the Bank of England. Basically, you need to solve the system of equations $q_1 = Q(0.1;\mu,\sigma,\epsilon)$; $q_2 = Q(0.5;\mu,\sigma,\epsilon)$, and $q_3 = Q(0.9;\mu,\sigma,\epsilon)$, where $Q$ is the quantile function given in Mudholkar and Hutson, equation (2.3).
First, you need to look at $(q_2-q_1,q_3-q_2)$. If $q_2-q_1>q_3-q_2$, then the density is left-skewed and $\epsilon>0$. Analogously, if $q_3-q_2 > q_2-q_1$, then the density is right-skewed and $\epsilon<0$. Once you identify this relationship, you can solve the system of equations using equation (2.3) from Mudholkar and Huston (2000). If you send me a cheque, I can do these calculations for you :).
$^*$The distribution in Mudholkar and Hutson belongs to a family of distributions called "two-piece", whose invention goes back to Fechner (1897). See:
The Two-Piece Normal, Binormal, or Double Gaussian Distribution: Its Origin and Rediscoveries
Best Answer
You can use MLE (maximum likelihood estimation), see also the links in comments, where there are given some formulas for method of moments. The easy way out is to use the
snpackage in R (on CRAN).snuse MLE (maximum likelihood estimation).Since you didn't give any data or context, I will just simulate some data, and show the use of the package: