Integrating the pdf of a normal distribution I need to find the distribution of $Y=X_1+X_2$ where both $X_1$ and $X_2$ are normally distributed with $(\mu,\sigma^2)$. 
So I'm looking for $f_y(y)=P(Y=y)$$=P(X_1+X_2=Y)$$=P(X_1=x_1)P(X_2=Y-x_1)$$$=\int_{x=-\infty}^{y} \frac1{\sqrt{2\pi}\sigma} e^{\frac{(x_1-\mu)^2}{2\sigma^2}} \frac1{\sqrt{2\pi}\sigma} e^{\frac{(y-x_1-\mu)^2}{2\sigma^2}} dx$$
$$=\int_{x=-\infty}^{y} \frac1{2\pi\sigma^2} e^{\frac{(x_1-\mu)^2+(y-x_1-\mu)^2}{2\sigma^2}}  dx$$$$= \frac1{2\pi\sigma^2} \int_{x=-\infty}^{y} e^{\frac{(x_1-\mu)^2+(y-x_1-\mu)^2}{2\sigma^2}} dx $$
which is about where I get stuck. I tried expanding the brackets and removing any terms from the intergral that don't depend on x and what I get is:
$$= \frac1{2\pi\sigma^2} e^{\frac{y^2+2\mu^2-2y\mu}{2\sigma^2}}\int_{x=-\infty}^{y} e^{\frac{x_1(x_1-y)}{\sigma^2}} dx $$
I've been looking for a way to use u substitution on the exponent term $\frac{x_1(x_1-y)}{\sigma^2}$ but I can't come up with anything that works. Am I missing something obvious? Any nudges in the right direction are welcomed. 
 A: Have you learned about the characteristic function of a random variable? The characteristic function of a sum of independent random variables, in general, is the product of the characteristic functions of the components of the sum. So if we call the characteristic function of a variable $X$ as $\Phi_X$, then $\Phi_{X+Y} = \Phi_X \Phi_Y$.
The characteristic function of a Normal random variable $X$ with mean $\mu_X$ and variance $\sigma_X$ is $$\Phi_X(\nu) = exp\left({j \nu \mu_X - \frac{\sigma_X^2 \nu^2}{2}}\right) $$ where $j = \sqrt{-1}$. So the sum of two Normal random variables, $X$ as above and $Y$ with mean $\mu_Y$ and variance $\sigma_Y$ is $$\Phi_{X+Y}(\nu) = exp\left({j \nu \mu_X - \frac{\sigma_X^2 \nu^2}{2}}\right)exp\left({j \nu \mu_Y - \frac{\sigma_Y^2 \nu^2}{2}}\right).$$ From there you can find your answer.
A: The sum of independent (a necessary assumption you left out) normal random variables has a normal distribution.  Well, maybe the purpose of the exercise is to prove that.
It makes things simpler if you first standardize your random variables: $X_i = \mu + \sigma Z_i$ where $Z_i$ have standard normal distribution (mean $0$ and variance $1$).
Your line $f_y(y)=P(Y=y)=P(X_1+X_2=Y)=P(X_1=x_1)P(X_2=Y−x_1)$ is all wrong.  These are continuous random variables, and you're looking for a density, not a probability.
The correct statement is 
$$f_Y(y) = \int_{-\infty}^\infty f_{X_1}(x) f_{X_2}(y-x)\ dx$$
A: Following up on Robert Israel's suggestion to use $X_i = \mu + \sigma Z_i$ where the $Z_i$ are independent standard normal random variables, note that $Y = X_1+X_2 = 2\mu + \sigma (Z_1+Z_2)$. So you could just find the density of $Z_1+Z_2$ and then transform it into the density of $X_1+X_2$ using the
standard result for linear transformations, viz.
$$f_{a+bW}(w) = \frac{1}{|b|}f_W\left(\frac{w - a}{b}\right)$$
which gives
$$f_{X_1+X_2}(y) = f_{Y}(y) = \frac{1}{\sigma}f_{Z_1+Z_2}\left(\frac{y-2\mu}{\sigma}\right).$$
This will save you a lot of carrying around $\mu$ and $\sigma$ which merely
clutter up the formulas.
A calculation of the density of $\alpha Z_1 + \beta Z_2$ for independent standard normal random variables can be found in the answers to 
this earlier question. Here,
of course, you have $\alpha = \beta = 1$. Note
that convolution is not explicitly required to obtain the answer in this
simple case; an appeal to the circular symmetry of the joint density of
$Z_1$ and $Z_2$ allows one to deduce the answer easily.
