Transforming a single Gaussian into a mixture of Gaussians Given a standard normal random variable: 
$X\sim \mathcal{N}(0,1)$
and given two real positive numbers $\mu$ and $\sigma$, I am interested in finding a (non-random) function $g:\mathbb{R}\rightarrow\mathbb{R}$ such that the random variable $g(X)$ is a mixture of two Gaussian:
$Y \triangleq g(X)\sim \frac{1}{2}\mathcal{N}(-\mu,\sigma^2)+\frac{1}{2}\mathcal{N}(+\mu,\sigma^2)$
(Where this notation means that the probability density function of $Y$ is half the sum of the probability density functions of $\mathcal{N}(-\mu,\sigma^2)$ and $\mathcal{N}(+\mu,\sigma^2)$ ).
I'm also interested in the more general setting, where we work with Gaussian random vectors, and we want to transform a standard vector into a mixture of $n$ Gaussians with given mean vectors, given covariance matrices and given weights (a weight is the coefficient of a Gaussian in the mixture, $\frac{1}{2}$ in the simpler problem above).
My questions are:


*

*Is this possible?

*If no, why not? If yes, what is the function $g$ that achieves this, or how can it be constructed?



My attempt at solving the simpler version of the problem was to use the change of variables formula for differentiable monotonic functions $g$: 
$f_Y(y) = \left| \frac{d}{dy} \big(g^{-1}(y)\big) \right| \cdot f_X\big(g^{-1}(y)\big)$,
to plug in my required $f_X$ and $f_Y$ and to solve for $g^{-1}$. However, this results in a complicated ordinary differential equation which I don't know how to solve (although I can attempt to solve it numerically).

I can very roughly construct such a function $g$, but not accurately (in the sense that it does not give two Gaussians, only two "blobs"). 
We can take:
$g(x)=\frac{(2\sigma x +\mu)e^{\lambda x}+(2\sigma x - \mu)e^{-\lambda x}}{ e^{\lambda x}+e^{-\lambda x} }$. 
This function look like this (for $\mu=2\,\ \sigma=1,\ \lambda=7$):

The idea is to "stretch" as much as possible the area around $x=0$ (hence the large derivative over there), so that the center of the original Gaussian gets spread along a large area and therefore its mass is reduced, while simply "translating" the mass in the other areas either to the left or to the right, to maintain the Gaussian shape...
This is the resulting distribution:

As you can see, I didn't exactly get Gaussians...
 A: To answer you question (1): Of course it is possible: since your original $X$ has a continuous distribution function $F$, you can use the $Y=G^{-1}(F(X))$ trick to generate a random variable $Y$ with any desired distribution function $G$.  In your case this  $g(t)=G^{-1}(F(t))$ function will be continuous and monotone increasing.
A: The random variable $Y$ is obtained from $X \sim \mathcal{N}\left( \mu,\sigma^{2}\right)$ by random sampling (with probability in your case) from the normal distribution :
\begin{equation*}
Y=\left\{ 
\begin{array}{c}
X \text{ with probability }1/2 \\ 
-X \text{ with probability }1/2
\end{array}
\right.
\end{equation*}
This fully describes the relationship between $X$ and $Y$ and as $g$ is not a one-to-one transformation, I am not convinced that the "change-of-variable technique" is useful here. The "distribution function technique", instead directly gives use the cdf and pdf of $Y$ (which fully describes this random variable):
\begin{eqnarray*}
F\left( y;\mu,\sigma \right) &\equiv & \Pr \left[ Y\leq y\right] \\
&=& 1/2 \Pr \left[ X\leq y\right] + 1/2 \Pr \left[ -X\leq y\right] \\
&=& 1/2 \int_{-\infty }^{y}\phi \left( u;\mu,\sigma^{2}\right) du+
1/2 \int_{-y }^{\infty 
}\phi \left( u;\mu,\sigma^{2}\right) du,
\end{eqnarray*}
\begin{eqnarray*}
f\left( y;\mu,\sigma^{2} \right)
&=&1/2 \phi \left( y;\mu,\sigma^{2}\right) + 1/2 \phi \left( -y;\mu,\sigma^2\right)  \\
&=&\frac{1/2}{\sqrt{2\pi }\sigma}\exp \left[ -\frac{1}{2\sigma^{2}
}\left( y-\mu\right) ^{2}\right] +\frac{1/2 }{\sqrt{
2\pi }\sigma}\exp \left[ -\frac{1}{2\sigma^{2}}\left( -y-\mu
\right) ^{2}\right].
\end{eqnarray*}
