Find independent $X$ and $Y$ that are not normally distributed s.t. $X+Y$ is normally distributed. It's a rather elementary fact that a sum of two independent normally distributed random variables $X$ and $Y$ is normally distributed (or, if you will, the convolution of two normal densities is a normal density). To what extent does this go the other way around? It seems that is $X$ and $Y$ are independent and $X$ is normally distributed, then $Y$ is normally distributed or constant.

If you drop both independence of $X$ and $Y$ and $X$ being normally distributed, it's pretty easy to e.g. take $X \sim \mathcal{N}(0,1)$, $Y = 1(X \geq 0)$ and consider $U=XY$, $V=X \,1(Y = 0)$. Now neither $U$ or $V$ are normally distributed but $U+V = X$.

However, if you require $X$ and $Y$ to be independent, but drop the requirement of $X$ being normally distributed, it's seems more difficult. Is there a counterexample in that case, where $X$ and $Y$ are not normally distributed but $X+Y$ is?
According to a theorem conjectured by P. Lévy and 
proved by H. Cramér (see Feller, Chapter XV.8, Theorem 1), 


  If $X$ and $Y$ are independent random variables and
  $X+Y$ is normally distributed, then both $X$ and $Y$ 
  are normally distributed.


I assume that $Y$ being a constant (and hence independent of $X$)
is considered as being accounted for in this theorem by
considering $Y$ to be a degenerate normal random variable with
variance $0$.
Consider the stable distribution, which I will denote as $\text{Stable}(\alpha, \beta, c, \mu)$. Let $X_1 \sim \text{Stable}(2, -1, \frac{1}{2}, 0)$ and $X_2 \sim \text{Stable}(2, 1, \frac{1}{2}, 0)$. Then the characteristic function of $X_1 + X_2$ is
\begin{align*}
\varphi_{X_1 + X_2}(t) &amp;= \varphi_{X_1}(t)\varphi_{X_2}(t) \\ &amp;= \exp\left(-\frac{1}{4}t^2(1+i\operatorname{sgn}(t)\Phi_2)\right)\exp\left(-\frac{1}{4}t^2(1-i\operatorname{sgn}(t)\Phi_2)\right) \\&amp;= \exp\left(-\frac{1}{2}t^2\right)
\end{align*}
which is the characteristic function of $N(0, 1)$. Note that both $X_1$ and $X_2$ are with $\alpha = 2$, which corresponds to the normal distribution. In general, such example doesn't exist for the statement to hold (see Feller's $\textit{An Introduction Probability Theory and Its Applications}$ Volume II Chapter XV.8 Theorem 1 for a proof) .

Edits and Remarks: Another technique I considered was the cumulant generating function of a normal, which takes the form $K(t) = at + bt^2$ (a quadratic form without a constant term). Since the cumulant of a sum of RV's = the sum of the cumulants, it suffices to consider a RV's such that $K_{X_1}(t) = a_1t + b_1t^2 + h(t)$ and $K_{X_2}(t) = a_1t + b_1t^2 - h(t)$, where $h(t)$ must satisfy $h(0) = 0$. If we restrict $h(t)$ to be the set of polynomials this suffices to just letting $h(t)$ be a degree 3 or higher polynomial. Unfortunately, such a distribution is then impossible. We could then let $h(t)$ include fractional powers &lt; 2, but the resulting distribution (inverse Laplace transform) would be super weird and I'm not sure if we can even define them in terms of elementary functions.