Solution 1:

The usual way to do this is to consider the moment generating function, noting that if $S = \sum_{i=1}^n X_i$ is the sum of IID random variables $X_i$, each with MGF $M_X(t)$, then the MGF of $S$ is $M_S(t) = (M_X(t))^n$. Applied to the exponential distribution, we can get the gamma distribution as a result.

If you don't go the MGF route, then you can prove it by induction, using the simple case of the sum of the sum of a gamma random variable and an exponential random variable with the same rate parameter. Let's actually do this. Suppose $Y \sim {\rm Gamma}(a,b)$ and $X \sim {\rm Exponential}(b)$ are independent, so that $$f_Y(y) = \frac{b^a y^{a-1} e^{-by}}{\Gamma(a)} \mathbb 1(y > 0), \quad f_X(x) = be^{-bx} \mathbb 1(x > 0), \quad a, b > 0.$$ Then, we notice that if $a = 1$, $Y$ would also be exponential (i.e., the exponential distribution is a special case of the Gamma with $a = 1$). Now consider $Z = X+Y$. The PDF is $$\begin{align*} f_Z(z) &= \int_{y=0}^z f_Y(y) f_X(z-y) \, dy \\ &= \int_{y=0}^z \frac{b^{a+1} y^{a-1} e^{-by} e^{-b(z-y)}}{\Gamma(a)} \, dy \\ &= \frac{b^{a+1} e^{-bz}}{\Gamma(a)} \int_{y=0}^z y^{a-1} \, dy \\ &= \frac{b^{a+1} e^{-bz}}{\Gamma(a)} \cdot \frac{z^a}{a} = \frac{b^{a+1} z^a e^{-bz}}{\Gamma(a+1)}. \end{align*}$$ But this is just a gamma PDF with new shape parameter $a^* = a+1$. So, it is easy to see by induction that the sum of $n$ IID exponential variables with common rate parameter $\lambda$ is gamma with shape parameter $a = n$, and rate parameter $b = \lambda$.