There are a few points to be addressed in your question.
First of all, exponential distributions are supported on the entirety of the positive real line, meaning that $X_1, X_2$ take values in $[0,\infty)$, rather than $[0,60]$ as you claim; moreover their sum $X = X_1 + X_2$ also takes values in $[0,\infty)$.
There are two immediate approaches to calculate the variance of $X$. The first one depends only on the fact that they are independent. A basic fact in probability theory asserts that if $U, V$ are independent random variables, then
\begin{align*}
\text{Var}(U+V) &= \mathbf E[ (U+ V)^2] - \mathbf E[U + V]^2 \\
& = \mathbf E[U^2] + \mathbf E[V^2] + 2\mathbf E[U]\mathbf E[V] - (\mathbf E[U]^2 + \mathbf E[V]^2 + 2\mathbf E[U]\mathbf E[V] )\\
& = \text{Var}[U] + \text{Var}[V]
\end{align*}
From this it follows from the fact that the variance of an $\text{Exp}(\lambda)$ variable is $\lambda^{-2}$, that
\begin{align*}
\text{Var}(X_1 + X_2) & = \lambda_1^{-2} + \lambda_2^{-2} \\
& = \frac{101}{4}.
\end{align*}
for $\lambda_1 = 1/5$, $\lambda_2 = 2$.
Note that in this approach we did not need any properties of the distributions, other than knowledge of their variances (i.e. if you gave me two distributions $U,V$, with $\text{Var}[U] = \lambda_1, \, \text{Var}[V] = \lambda_2$, the answer would not change).
A second approach would be to argue via the probability density functions (p.d.f.s), as you suggest in the question. This will in general be a harder approach, as we will have to derive the variance of the sum from the p.d.f., as well as computing the p.d.f.
In general for independent $U, V$, with p.d.f. $f_U, f_V$ respectively, the distribution of their sum is given as the convolution of their p.d.f.s
\begin{align*}
f_{U +V}(x) = \int_{-\infty}^\infty f_U(y)f_V(x-y)dy.
\end{align*}
In the case of a exponential distributions, the p.d.f. of a $\lambda$ exponential distribution is
\begin{align*}
f(x) = \lambda e^{-\lambda{x}}
\end{align*}
And so for two independent exponentials with parameters $ \lambda_1 \neq \lambda_2$
\begin{align*}
f_{X}(x) &= f_{X_1 + X_2}(x) \\
& = \int_{0}^x \lambda_1 e^{-\lambda_1 y} \lambda_2 e^{-\lambda_2(x-y)}dy \\
& = \lambda_1 \lambda_2 e^{-\lambda_2 x}\int_{0}^x e^{(\lambda_2-\lambda_1)y}dy \\
& = \frac{\lambda_1 \lambda_2}{\lambda_2-\lambda_1}( e^{-\lambda_1 x}-e^{-\lambda_2 x}).
\end{align*}
Note that we only integrate over the region supported by the p.d.f.s. From this we can then proceed to calculate the variance from the p.d.f. by the standard formulae
\begin{align*}
\mathbf E[V] &= \int_{-\infty}^\infty x f_V(x) d x,\\
\mathbf E[V^2] &= \int_{-\infty}^\infty x^2 f_V(x) dx
\end{align*}
In all, I hope that my post has made it clear that this second approach requires a lot more work, and is more prone to calculation error!
Suppose that $\left(X_{i}\right)_{1\leq i\leq n}$
is independant and identically distributed according to an exponential law with a parameter $\lambda>0$
. For all $x\in\mathbb{R}$
, we have$$\mathbb{P}\left[X_{i}\leq x\right]=1-e^{-\lambda x}.$$
Now, it is straightforward to use the characteristic functions : if $S_{n}=\sum_{i=1}^{n}X_{i}$
, then we have for all $t\in\mathbb{R}$
$$\chi_{_{S_{n}}}\left(t\right)=\prod_{i=1}^{n}\chi_{_{X_{i}}}\left(t\right)=\left(\chi_{_{X_{1}}}\left(t\right)\right)^{n}=\left(\left(1-i\lambda t\right)^{-1}\right)^{n}=\left(1-i\lambda t\right)^{-n}$$
which is the characteristic function of a random variable that follow a gamma law of parameters $\left(n,\lambda\right)$ (I used independance and identical distribution here). Thus, the density of the law of $S_{n}$
is the density of the gamma law (see for example http://en.wikipedia.org/wiki/Gamma_distribution).
You can also show by induction that the density of the sum of INDEPENDANT random variables is the convolution of the densities.
Best Answer
The mistake is in the upper bound in the very last integral. It should be $z$ because $x+y\leq z$ and since $Y$ is exponentially distributed it is nonnegative. Then we have
$$f_Z(z)=\lambda^2\int_0^ze^{-\lambda z}dx=\lambda^2 ze^{-\lambda z}$$
Note that this is the answer we expect since the sum of two exponential with parameter $\lambda$ is distributed as $\Gamma(2,\lambda)$ which has density
$$f_Z(z)=\frac{\lambda^2z^{2-1}e^{-\lambda z}}{\Gamma(2)}$$