$X_{1},X_{2},X_{3}overset{i.i.d.}{sim}N(0,1)$, find m.g.f. of $Y=X_{1}X_{2}+X_{1}X_{3}+X_{2}X_{3}$

In this situation it helps to keep the computations a little abstract, because that exposes their essential simplicity. Let, then, $Q$ be the homogeneous quadratic form $$Q(mathbf{x}) = x_1 x_2 + x_2 x_3 + x_3 x_1 = frac{1}{2}mathbf{x}^primepmatrix{0&1&1\1&0&1\1&1&0}mathbf{x} = frac{1}{2}mathbf{x}^primemathbb{Q}mathbf{x}.$$

By definition the mgf of any function $f$ is $E[exp(tf(X))],$ which for this form is

$$Eleft[expleft(t Q(X)right)right] = C intcdotsint expleft(t Q(mathbf{x})right)expleft(-frac{1}{2} mathbf{x}^primemathbf{x}right),mathrm{d}mathbf{x}.$$

The product of the exponentials is, of course, the exponential of the sum of the arguments, so let's focus on that sum, still keeping the details of $Q$ unexposed:

$$t Q(mathbf{x})-frac{1}{2} mathbf{x}^primemathbf{x} = -frac{1}{2}mathbf{x}^primeleft(mathbb{I} - tmathbb{Q}right)mathbf{x}.$$

It is immediate, just by inspection, that this gives an integrand proportional to the multivariate density of a Normal distribution of zero mean and covariance $left(mathbb{I} - tmathbb{Q}right)^{-1}.$ The constant of proportionality is the square root of the determinant of the covariance. Thus (at least for $|t|$ sufficiently small),

$$Eleft[expleft(t Q(X)right)right] = left|mathbb{I} - tmathbb{Q}right|^{-1/2}.$$

As a quick check, when $mathbb{Q}=2mathbb{I}$ represents the sum of squares $X_1^2+X_2^2+X_3^2,$ we should obtain the mgf for that sum of squares, which defines the $chi^2(3)$ distribution. Indeed, $$left|mathbb{I} - t(2mathbb{I})right|^{-1/2} = left|(1-2t)mathbb{I}right|^{-1/2} = ((1-2t)^3|mathbb{I}|)^{-1/2} = (1-2t)^{-3/2}$$ is the correct result according to https://en.wikipedia.org/wiki/Chi-square_distribution.

The answer for the particular form $Q$ in the question is

$$begin{aligned} left|mathbb{I} - tmathbb{Q}right|^{-1/2} &= left|pmatrix{1&-t&-t\-t&1&-t\-t&-t&1}right|^{-1/2} \ &= left(1 - 3t^2 - 2t^3right)^{-1/2} \ &= (1-2t)^{-1/2}(1+t)^{-2/2}. end{aligned}$$

At https://stats.stackexchange.com/a/318908/919, Kjetil Halvorsen diagonalizes $Q$ to generalize this calculation to inhomogeneous quadratic forms (giving non-central chi-squared distributions). I have offered the foregoing analysis because it is complete and more elementary (not requiring any matrix decomposition), thereby perhaps helping you to find where your original work went awry.

The post by StubbornAtom at https://stats.stackexchange.com/a/478682/919 presents the same generalization in a different way.

In your first line you write the expression as the sum of three quadratic forms. You then multiply the three mgfs. But, this is only possible if the variables (quadratic forms) are independent. Quadratic forms are independent if and only if the product of their matrices is the zero matrix. This is not true of the first and second quadratic forms. However, it is true of the first and third and it is true of the second and third. That is why you could multiply two of the mgfs and get a correct result for those.