Finding the probability density function of $Z = X + Y$ The problem
Let there be two independent random variables given. Let's name them $X$ and $Y$. Suppose the variables have both normal distribution $N(\mu, \lambda^2).$
Let's define third random variable $Z = X + Y$. The task is to find the density function of $Z$.
My attempt
I know that the density function can be easily calculated using the distribution function (I would need to differentiate it). That's why I started with finding the distribution.
$$F_Z(t) = Pr(Z \le t) = Pr(Y + X  \le t) = Pr(Y \le t - X) = \int \limits_{-\infty}^{\infty}Pr(Y \le t-x)f_X(x) dx$$
The equation above leads us to the following one:
$$ \int \limits_{-\infty}^{\infty} \int \limits_{-\infty}^{t-x} f_Y(y)f_X(x) \mbox{d} y\mbox{d}x \tag{1}.$$
Using the knowledge about $f_X$ and $f_Y$ we can write $(1)$ in this form:
$$\theta\int \limits_{-\infty}^{\infty} \int \limits_{-\infty}^{t-x}e^{-\frac{(y- \mu)^2}{\lambda^2}} e^{-\frac{(x- \mu)^2}{\lambda^2}} \mbox{d} y\mbox{d}x,$$
where $\theta = \frac{1}{\sqrt{2 \lambda^2}}.$
I wonder if my attempt was correct because the integral is quite horrible. Is it possible to calculate it in a "nicer" way?

Is $(1)$ equivalent to this expression $(f_X*f_Y)(x)$?
If yes I would get the following integral:
$$\int \limits_{\mathbb{R}} e^{-\frac{(y- \mu)^2}{\lambda^2}} e^{-\frac{(x- y - \mu)^2}{\lambda^2}} \mbox{d}y.$$
Is there a method of calculating such things?
 A: From your expression, $F_Z(t) = \int_\mathbb{R}P(Y \leq t-x)f_X(x)dx$, differentiate to obtain the density (you may need to justify why you can differentiate under the integral), 
$$ f_Z(t) = \int_\mathbb{R}f_Y(t-x)f_X(x)dx = \int_\mathbb{R}\frac{1}{\sqrt{2\pi}}e^{-(t-x-\mu)^2/2\lambda^2}\frac{1}{\sqrt{2\pi}}e^{-(x-\mu)^2/2\lambda^2}dx$$
which we recognize as $(f_Y*f_X) (x)$. To actually evaluate the integral, expand the exponent terms inside the integrand,
$$-\frac{1}{2\lambda^2}\left[(t-x-\mu)^2 + (x-\mu)^2\right] = -\frac{1}{2\lambda^2}(t^2-2\mu t + 2\mu^2 + \underbrace{2x^2-2tx}_{(*)})$$
Complete the square on the $(*)$ term, $2x^2 - 2tx = 2\left(x-\frac{t}{2}\right)^2-\frac{t^2}{2}$. Putting things together, 
$$-\frac{1}{2\lambda^2}(t^2-2\mu t + 2\mu^2 + 2x^2-2tx) = -\frac{1}{4\lambda^2}(t^2-4\mu t + 4\mu^2)-\frac{1}{\lambda^2}(x-t/2)^2$$
So, we can pull out the $t$ terms from the integrand to be left with,
$$ f_Z(t) = \frac{1}{\sqrt{4\pi\lambda^2}}e^{-(t-2\mu)^2/4\lambda^2}\underbrace{\int_\mathbb{R}\frac{1}{\sqrt{\pi\lambda^2}}e^{-(x-t/2)^2/\lambda^2}dx}_{=1} = \frac{1}{\sqrt{4\pi\lambda^2}}e^{-(t-2\mu)^2/4\lambda^2}$$
The integral above equals $1$ because the integrand is the density function of a $N(t/2,\lambda^2/2)$ distribution.
A: The problem can also be solved by using the following property:
If $X\sim N(\mu_1,\sigma_1^2) $ and $Y \sim N(\mu_2,\sigma_2^2)$ and are independent, then $aX + bY \sim N(a\mu_1 + b\mu_2,a^2\sigma_1^2 +b^2\sigma_2^2)$ and thus in your case $Z \sim N(2\mu,2\lambda^2)$.
Edit: I just noticed the comment by @saulspatz above who already pointed this out. 
