Computation of a Particular Convolution Let $\xi_{1}, \xi_{2}, \xi_{3}$ be i.i.d. $N(0,1)$. I'm attempting to compute the density of $\max \{\xi_{1}, \xi_{2}\} + \xi_{3}$. I know the density of $\max \{\xi_{1}, \xi_{2}\} $ is $2\Phi(y) \phi(y)$ and the density of $\xi_{3}$ is $\phi(s)$ so that the density of interest is the convolution $\int_{\mathbb{R}} 2\Phi(y) \phi(y) \phi(z-y)dy$. Is there any way to get at this expression, say, express it in terms of $\Phi$?
 A: The density of $\max\{\xi_1,\xi_2\}$ is given by 
$$g(t):=\frac 1{2\pi}e^{-t^2/2}\int_{-\infty}^te^{-s^2/2}ds$$
so we want to compute 
$$f(x):=\frac 1{2\pi}\frac 1{\sqrt{2\pi}}\int_{-\infty}^{+\infty}e^{-t^2/2}\int_{-\infty}^te^{-s^2/2}ds\cdot e^{-(x-t)^2/2}dt.$$
We can write 
\begin{align*}
e^{-t^2/2}e^{-(x-t)^2/2}&=\exp\left(-\frac 12(t^2+x^2-2xt+t^2)\right)\\\
&=\exp\left(-(t^2-xt)-\frac{x^2}2\right)\\\
&=\exp\left(-\frac{x^2}2\right)\exp\left(-\left(t-\frac x2\right)^2+\frac{x^2}4\right)\\\
&=\exp\left(-\frac{x^2}4\right)\exp\left(-\left(t-\frac x2\right)^2\right)
\end{align*}
hence 
\begin{align*}f(x)&=(2\pi)^{-3/2}\exp\left(-\frac{x^2}4\right)\int_{-\infty}^{+\infty}\exp\left(-\left(t-\frac x2\right)^2\right)\int_{-\infty}^te^{-s^2/2}dsdt\\\
&=(2\pi)^{-3/2}\exp\left(-\frac{x^2}4\right)\int_{-\infty}^{+\infty}\exp\left(-y^2\right)\int_{-\infty}^{y+x/2}e^{-s^2/2}dsdy.
\end{align*}
Put $h(x):=e^{x^2/4}(2\pi)^{3/2}f(x)$. We have 
\begin{align*}
h'(x)&=\frac 12\int_{-\infty}^{+\infty}\exp\left(-y^2\right)e^{-(y+x/2)^2/2}dy
\end{align*}
and 
\begin{align*}
\exp\left(-y^2\right)e^{-(y+x/2)^2/2}&=\exp\left(-\frac 32 y^2-xy-x^2/4\right)\\\
&=\exp\left(-\frac 32\left(y^2+\frac 23xy+\frac{x^2}6\right)\right)\\\
&=\exp\left(-\frac 32\left(\left(y+\frac 13x\right)^2-\frac{x^2}9+\frac{x^2}6\right)\right)\\\
&=\exp\left(-\frac{x^2}{12}\right)\exp\left(-\frac 32\left(y+\frac 13x\right)^2\right)
\end{align*}
hence 
\begin{align*}
h'(x)&=\frac 12\exp\left(-\frac{x^2}{12}\right)\int_{-\infty}^{+\infty}e^{-3t^2/2}dt\\\
&=\frac 12\frac 1{\sqrt 3}\exp\left(-\frac{x^2}{12}\right)\int_{-\infty}^{+\infty}e^{-s^2/2}ds\\\
&=\frac 12\frac 1{\sqrt 3\sqrt{2\pi}}\exp\left(-\frac{x^2}{12}\right).
\end{align*}
We deduce that 
$$f(x)=\frac 12(2\pi)^{-3/2}\frac 1{\sqrt 3\sqrt{2\pi}}e^{-x^2/4}\int_{-\infty}^x\exp\left(-\frac{t^2}{12}\right)dt.$$
