I can show you all examples, not just the simple ones.
Solution
Here they are, schematically:
The bottom panels show how the density function $f$ of a distribution $F$ is split into two parts vertically along a nearly arbitrary curve. The cyan portion of the split is a fraction $\lambda$ of $f;$ the upper left plots its graph. The remaining portion (gray) therefore is a fraction $1-\lambda$ of $f$ whose graph is plotted in the upper right. This is how all mixtures arise.
(Notice that little is assumed about the density $f$ except that it exists.)
Details
The mixture distribution $F$ is Normal, which means there is a mean $\mu$ and variance $\sigma^2$ for which $F$ has a density function $f(z;\mu,\sigma).$ The details of $f$ don't matter!
Let $\lambda:\mathbb{R}\to[0,1]$ be any (measurable) non-negative function. This means the following integrals involving $\lambda$ are defined and non-negative:
$$\pi_\lambda = \int_\mathbb{R} \lambda(z)f(z;\mu,\sigma)\,\mathrm{d}z \le \sup(\lambda)\, \int_\mathbb{R}f(z;\mu,\sigma)\,\mathrm{d}z \le (1)(1)=1;$$
$$1-\pi_\lambda = 1 - \int_\mathbb{R} \lambda(z)f(z;\mu,\sigma)\,\mathrm{d}z = \int_\mathbb{R} (1-\lambda(z))_f(z;\mu,\sigma)\,\mathrm{d}z \le 1.$$
(The first inequality is an easy special case of Holder's Inequality.)
Define two distributions as
$$F_{\lambda}(x) = \frac{1}{\pi_\lambda}\int_{-\infty}^x \lambda(z)f(z;\mu,\sigma)\,\mathrm{d}z;$$
$$F_{1-\lambda}(x) = \frac{1}{1-\pi_\lambda}\int_{-\infty}^x (1-\lambda(z))f(z;\mu,\sigma)\,\mathrm{d}z.$$
It is straightforward to establish that these are distribution functions and, by construction,
$$F = \pi_\lambda F_\lambda + (1-\pi_\lambda) F_{1-\lambda}\tag{*}$$
exhibits the original normal distribution as a mixture of these two.
Conversely, whenever there exist differentiable functions with property $(*),$ then a version of $\lambda$ can be recovered via
$$\lambda(z) = \left\{\begin{aligned}\frac{F^\prime_\lambda(z)}{f(z;\mu,\sigma)} &\quad&f(z;\mu,\sigma)\ne 0\\
0 & &\text{otherwise}\end{aligned}\right.$$
and because $0 \le \pi_\lambda\le 1,$ the range of $\lambda$ is contained in $[0,1],$ QED.
Finally, it is possible for the component distributions to be Normal: for instance, when $\lambda$ is a constant function that will be the case. That is the only possibility, though: see https://stats.stackexchange.com/a/429877/919 for the proof.
Application
As requested in comments, it would be of interest to choose $\lambda$ to meet a set of criteria, such as
Give the components equal weights, which means $$\frac{1}{2}=\pi_\lambda = \int \lambda(z) f(z)\,\mathrm{d}z.$$
Since these are intended to model errors in a regression setting (with $\mu=0,$ we would like each of the components also to have zero mean: $0 = E_{F_\lambda}[X].$ In light of (1), that is equivalent to $$0 = \int z\lambda(z) f(z)\,\mathrm{d}z.$$
Since regression errors are often assumed to be homoscedastic -- of equal variances -- we would like the variances of $F_\lambda$ and $F_{1-\lambda}$ to be equal. Since they have means of zero, when $f$ is a Normal density, this is achieved when $$\sigma^2 = 2\int z^2\lambda(z) f(z)\,\mathrm{d}z.$$
Although there are many solutions to these equations, one simple (striking) solution is obtained by supposing $\lambda$ and $1-\lambda$ are both simple functions: that is, piecewise constant. By making $\lambda$ symmetric around $0$ we can assure that (2) holds. The simplest of such simple functions is zero except on some positive interval $[a,b]$ and its negative $[-b,-a],$ where it equals $1.$
Without any loss of generality take $\sigma^2=1,$ so that $f = \phi$ is the standard Normal density with the property $\phi^(z) = -z\phi(z).$ Using this fact we may compute
$$\int \lambda(z)\phi(z)\,\mathrm{d}z = 2 \int_a^b \phi(z)\,\mathrm{d}z = 2(\Phi(b)-\Phi(a))$$
(where $\Phi$ is the standard Normal distribution function) and
$$\begin{aligned}
\int z^2 \lambda(z)\phi(z)\,\mathrm{d}z &= 2 \int_a^b z^2\phi(z)\,\mathrm{d}z \\
&= 2(\Phi(b) - \Phi(a) + a\phi(b) - b\phi(b)).
\end{aligned}$$
This permits numerical solution of (1) and (3). The work is streamlined by noting from (1) that, given $0 \le a\lt \Phi^{-1}(3/4),$
$$b = b(a) = \Phi^{-1}(\Phi(a) + 1/4).$$
That leaves us to solve (3) for $a \ge 0$. Here is an R
implementation to illustrate:
f <- function(a) {
b <- qnorm(1/4 + q <- pnorm(a))
pnorm(b) - q + a * dnorm(a) - b * dnorm(b) - 1/4
}
uniroot(f, c(0, qnorm(3/4)- 1e-6))$root -> a
qnorm(pnorm(a) + 1/4) -> b
This calculation gives $a \approx 0.508949$ and $b \approx 1.59466.$ Here are plots of the two component densities $f_\lambda$ and $f_{1-\lambda}:$
To illustrate the intended application, here are bivariate data with 150 responses at $X=0$ with errors distributed as $F_\lambda$ and 150 responses at $X=1$ with errors distributed as $F_{1-\lambda}.$ To the right is a quantile plot of the collected residuals.
Although separately neither group of residuals appears Normal, they are both centered at zero, have nearly the same variance, and collectively look perfectly Normal.
Remarks
The basic construction readily generalizes to mixtures with more than two components.
The example in the application can be extended, by using simple (indicator) functions supported on intervals $[a_i,b_i]$ with $0\le a_1 \lt b_1 \le a_2 \lt b_2 \cdots \lt b_k,$ to create component distributions that match the first $2k$ moments of the Normal distribution their mixture creates. With sufficiently large $k,$ the component distributions will be difficult to discriminate even with largish datasets (at which point one might legitimately wonder whether their non-Normality matters at all).