Let's try using the moment-generating functions of $X$ and $Y$
Let $\lambda_1=\frac{1}{\theta_1}$ and $\lambda_2=\frac{1}{\theta_2}$
Then $M_{X+Y}(t)=M_X(t)M_Y(t)=\frac{\lambda_1}{\lambda_1-t}\cdot\frac{\lambda_2}{\lambda_2-t}=\frac{\lambda_1\lambda_2}{\lambda_1\lambda_2-\lambda_1t-\lambda_2t+t^2}=\dots$
I'm not sure what to do next. Is it possible to manipulate this expression into the moment-generating function of some well-known probability distribution, or is there a better approach?
The sum of independent but not identically distributed exponential random variables is not exponential. Let's do the calculation via the CDFs: define $Z = X + Y$, and suppose $z \ge 0$, so that $$F_Z(z) = \Pr[Z \le z] = \Pr[X + Y \le z] = \int_{y=0}^\infty \Pr[X \le z - y]f_Y(y) \, dy = \int_{y=0}^z F_X(z-y) f_Y(y) \, dy.$$ Continuing, we obtain $$\begin{align*} F_Z(z) &= \frac{1}{\theta_2} \int_{y=0}^z (1 - e^{-(z-y)/\theta_1}) e^{-y/\theta_2} \, dy \\ &= \frac{1}{\theta_2} \int_{y=0}^z e^{-y/\theta_2} \, dy - \frac{1}{\theta_2} \int_{y=0}^z e^{-z/\theta_1} e^{-y(1/\theta_2 - 1/\theta_1)} \, dy \\ &= F_Y(z) - \frac{e^{-z/\theta_1}}{\theta_2} \frac{1}{\frac{1}{\theta_2} - \frac{1}{\theta_1}} \left(1 - e^{-z(1/\theta_2 - 1/\theta_1)}\right) \\ &= 1 - e^{-z/\theta_2} - \frac{\theta_1}{\theta_1 - \theta_2} \left( e^{-z/\theta_1} - e^{-z/\theta_2} \right) \\ &= 1 - \frac{\theta_1 e^{-z/\theta_1} - \theta_2 e^{-z/\theta_2}}{\theta_1 - \theta_2}. \end{align*}$$ Differentiating with respect to $z$ yields the density $$f_Z(z) = \frac{e^{-z/\theta_1} - e^{-z/\theta_2}}{\theta_1 - \theta_2}.$$ When $\theta_1 = \theta_2$, then we get a gamma distribution with shape $2$ and scale $\theta_1 = \theta_2$.
