I want to find the density function of the random variable $Y=X_1+X_2+X_3$, where the random variables $X_1, X_2$ and $X_3$ have a joint density function $$f_{X_1, X_2, X_3}(x_1, x_2, x_3) = (2\pi)^{-\frac{3}{2}}e^{-\frac{1}{2}(x_1^2+x_2^2+x_3^2)},\qquad -\infty<x_1, x_2,x_3<\infty.$$ It is a known result that the sum of three independent normally distributed random variables with mean 0 and variance 1 is normal, with its mean being 0 and variance being 3, so that we know $f_Y(y)=\frac{1}{\sqrt{6\pi}}e^{-\frac{y^2}{6}},\quad -\infty<y<\infty.$
However, I want to verify this by finding the derivative of $F_Y(y):=P(X_1+X_2+X_3\leq y)$ with respect to $y$. In other words, I need to find $$\frac{d}{dy}\Bigg[\int_{-\infty}^{\infty}\int_{-\infty}^{y-x_3}\int_{-\infty}^{y-x_2-x_3}(2\pi)^{-\frac{3}{2}}e^{-\frac{1}{2}(x_1^2+x_2^2+x_3^2)}dx_1dx_2dx_3\Bigg].$$
My questions are:
(1) Are the upper and lower bounds correct?
(2) How can I evaluate this expression? I have tried applying the Leibniz's rule but end up getting something impossible to integrate. I have also looked at the Reynolds transport theorem, but I don't quite understand it, especially the region is unbounded in this case.
Any help is greatly appreciated.
Since no one answers it for 2 years and now I know where I went wrong, I'll post my answer here.
The short answer is: (1) The bounds are incorrect. (2) Once we fix that, the expression can be evaluated quite easily using Leibniz's rule and the Gaussian integral.
Let's consider the case with 2 random variables first, i.e. $Y = X_1 + X_2$. In this case, we do have $$\frac{d}{dy} F_Y(y)=\frac{d}{dy}P(X_1+X_2\leq y)=\frac{d}{dy}\Bigg[\int_{-\infty}^{\infty}\int_{-\infty}^{y-x_2}(2\pi)^{-1}e^{-\frac{1}{2}(x_1^2+x_2^2)} \: dx_1dx_2\Bigg].$$
We then use the Leibniz integral rule for differentiation under the integral sign. Noting that the upper bound of the inner integral involves $y$ but the integrand does not involve $y$, we get
\begin{align} \frac{1}{2\pi}\frac{d}{dy}\Bigg[\int_{-\infty}^{\infty}\int_{-\infty}^{y-x_2}e^{-\frac{1}{2}(x_1^2+x_2^2)} \: dx_1dx_2\Bigg] &= \frac{1}{2\pi}\int_{-\infty}^{\infty}\Bigg[e^{-\frac{1}{2}[(y-x_2)^2+x_2^2]}\cdot 1+\int_{-\infty}^{y-x_2}0 \: dx_1\Bigg]dx_2 \\ &= \frac{1}{2\pi} \int_{-\infty}^{\infty}e^{-\frac{1}{2}[(y-x_2)^2+x_2^2]} \: dx_2 \\ &= \frac{1}{2\pi} \int_{-\infty}^{\infty}e^{-(x_2-\frac{y}{2})^2-\frac{y^2}{4}} \: dx_2 \\ &= \frac{1}{2\pi} e^{-\frac{y^2}{4}} \int_{-\infty}^{\infty}e^{-(x_2-\frac{y}{2})^2} \: dx_2 \\ &= \frac{1}{2\pi} e^{-\frac{y^2}{4}} \sqrt{\pi} \\ &= \frac{1}{\sqrt{2} \cdot \sqrt{2\pi}} e^{-\frac{1}{2}(\frac{y}{\sqrt{2}})^2}, \end{align}
where the second last equality is obtained by a simple change of variables and using the well-known Gaussian integral. We see that the resulting function is indeed the density function of the normal distribution with mean 0 and variance 2.
Now consider the original 3 variables problem, i.e. $Y = X_1+X_2+X_3$.
To see why the bounds that the OP (that's me) gave are incorrect, consider the point $P := (x_1, x_2, x_3) := (-1, y, 1) \in \mathbb{R}^3$. We have $x_1+x_2+x_3 \leq y \:$. (In fact, $P$ lies on the plane$X_1+X_2+X_3=y$). However, it is false that $-\infty < x_2 \leq y-x_3$, as described by OP's bounds.
The reason why I made this mistake was that I was too used to the standard Calculus 3 exercises where we are asked to find the volume of a $\textbf{bounded}$ solid. In our case here, we really don't have any constraints on $x_2$ and $x_3$: given any $x_2$ and $x_3$, we only need to make sure that $-\infty < x_1 \leq y - x_2 - x_3$. Therefore, very similar to the 2 variables case, we get:
\begin{align} \frac{d}{dy} F_Y(y) &=\frac{d}{dy}P(X_1+X_2+X_3\leq y) \\ &=\frac{d}{dy}\Bigg[\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}\int_{-\infty}^{y-x_2-x_3}(2\pi)^{-\frac{3}{2}}e^{-\frac{1}{2}(x_1^2+x_2^2+x_3^2)} \: dx_1dx_2dx_3\Bigg] \\ &=(2\pi)^{-\frac{3}{2}} \int_{-\infty}^{\infty}\int_{-\infty}^{\infty} e^{-\frac{1}{2}[(y-x_2-x_3)^2+x_2^2+x_3^2]} \: dx_2dx_3\\ &=(2\pi)^{-\frac{3}{2}} e^{-\frac{1}{6}y^2} \int_{-\infty}^{\infty} e^{-\frac{3}{4}\big(x_3-\frac{y}{3}\big)^2} dx_3 \int_{-\infty}^{\infty} e^{-\big(x_2 + \frac{x_3-y}{2}\big)^2} dx_2 \\ &=(2\pi)^{-\frac{3}{2}} e^{-\frac{1}{6}y^2}\bigg(\frac{2}{\sqrt{3}}\sqrt{\pi} \bigg) \sqrt{\pi} \\ &=\frac{1}{\sqrt{3} \cdot \sqrt{2\pi}} e^{-\frac{1}{2}(\frac{y}{\sqrt{3}})^2} \end{align}
So there we have it, the density function of the normal distribution with mean 0 and variance 3.