Distribution of $x_4(x_1-x_3)+x_5(x_2-x_1)$ with iid $x_i \sim N(0,1)$ Problem
Given are 5 independent standard normal variables $x_1,x_2,x_3,x_4,x_5$.
What is the pdf of $x_4(x_1-x_3)+x_5(x_2-x_1)$ ?
What I know
$$x_1-x_3\sim \mathcal{N}\left(0,\sqrt{2}\right)\tag{1}$$
$$x_2-x_1\sim \mathcal{N}\left(0,\sqrt{2}\right)\tag{2}$$
$$x_4(x_1-x_3)\sim \frac{1}{\pi \sqrt{2}}K_0\left(\frac{\left|z\right|}{\sqrt{2}}\right)\tag{3}$$
$$x_5(x_2-x_1)\sim \frac{1}{\pi \sqrt{2}}K_0\left(\frac{\left|z\right|}{\sqrt{2}}\right)\tag{4}$$
where $K_0$ is the Bessel function and eq.(3,4) is a normal product distribution.
Remains to add the 2 summands $x_4(x_1-x_3)$ and $x_5(x_2-x_1)$ but they are dependent and cannot be convolved.
Related problem from literature
If two variables are transformed like $x_1(a x_1+b x_2)$ with $a\in \mathbb{R},b\in\mathbb{R}^+$ then the pdf is known $[1]$ but this is not applicable here.
Simulation
Simulation shows that the distribution of $x_4(x_1-x_3)+x_5(x_2-x_1)$ can be approximated with a Laplace distribution with parameter $(0,1.34)$. But the сorrect answer is not a Laplace distribution.

$[1]$ R. Gaunt: A note on the distribution of the product of zero-mean correlated normal random variables, Statistica Neerlandica, 2018

An almost exact expression could be derived using link between $\chi^{2}$, $\Gamma$ and symmetric $VG$ variance gamma distributions. Given very useful results above by @whuber, we can proceed first by the link $\chi^{2}(1)\sim\Gamma(\frac{1}{2},2)$. Therefore, from result that the product in question is distributed as $\frac{1}{2}[X_{1}-X_{2}]+\frac{\sqrt{3}}{2}[Y_{1}-Y_{2}]$ and knowing that $X_{1,2}\sim\chi^{2}(1)$, $Y_{1,2}\sim\chi^{2}(1)$ and $\chi^{2}(1)\sim\Gamma(\frac{1}{2},2)$. We will have,
$$\begin{equation}
\frac{1}{2}[X_{1}-X_{2}]+\frac{\sqrt{3}}{2}[Y_{1}-Y_{2}]\sim\left[\Gamma(\frac{1}{2},1)-\Gamma(\frac{1}{2},1)+\Gamma(\frac{1}{2},\sqrt{3})-\Gamma(\frac{1}{2},\sqrt{3})\right]
\end{equation}$$
which is the sum of difference of $\Gamma(\alpha,\nu)$ rvs. Now, symmmetric $VG$, which was introduced by Madan and Seneta (1990) is a distribution with $\sigma W\sqrt{\Gamma(\alpha_{G},\nu_{G})}$ where $W\sim N(0,1)$ and could be written as a difference of two $\Gamma$ distributed rvs. This has a characteristic function (CF)
$$\begin{equation}
\phi(u,\sigma,\nu)=\left(\frac{1}{1+\frac{\nu u^{2}\sigma^{2}}{2}}\right)^{\frac{1}{\nu}}
\end{equation}$$
Difference of $\Gamma(\alpha,\nu)$ rvs will have CF
$$\begin{equation}
\phi(u,\alpha,\nu_{G})=\left(\frac{1}{1+\nu_{G}^{2}}\right)^{\alpha}
\end{equation}$$
Equating these two will give us
$$\begin{equation}
\sigma=\nu_{G}\sqrt{2\alpha}\\
\nu=\frac{1}{\alpha}
\end{equation}$$
Therefore, $\frac{1}{2}[X_{1}-X_{2}]+\frac{\sqrt{3}}{2}[Y_{1}-Y_{2}] \sim VG(\sigma=1,\nu=2)+VG(\sigma=\sqrt{3},\nu=2)$. Here we have some result, but still we need the convolution of these two $VG$ random variables. There is following convolution rule for $VG$ distribution, $VG(\sigma_{1},\nu_{1})+VG(\sigma_{2},\nu_{2})\sim VG\left(\sqrt{\sigma_{1}^{2}+\sigma_{2}^{2}},\frac{\nu_{1}+\nu_{2}}{\nu_{1}\nu_{2}}\right)$, which has additional constraint $\sigma_{1}^{2}\nu_{1}=\sigma_{2}^{2}\nu_{2}$. Here, we slightly diverge from this additional condition and use the first one, which results an approx. exact density. Therefore, final result is  $\frac{1}{2}[X_{1}-X_{2}]+\frac{\sqrt{3}}{2}[Y_{1}-Y_{2}] \sim VG(\sigma=2, \nu=1)$.
$$\begin{equation}
f_{X}\left(x\right)=\frac{x}{4\sqrt{\pi}}K_{\frac{1}{2}}\left(4\sqrt{2}|x|\right),
\end{equation}$$
After we check this result with MC simulation, we see that the result is fairly successful to have a closed form density and MC simulation is caught almost perfectly by both pdf and cdf.  
Hope this will be useful.
