Let $X,Y$ be independent exponential RV's with respective pdf's $f(x) = \lambda e^{-\lambda x}$ and $f(y) = \mu e^{-\mu y}$. We want to find the pdf of $Z=X-Y$.
I originally tried the convolution method (but with wrong limits). $$ P[X-Y \le z]= \int_0^{\infty} \Bigg( \int_0^{z+y} \lambda e^{-\lambda x} dx \Bigg) \mu e^{-\mu y} dy$$
which is wrong.
I was told that I'm supposed to split the inner integral into two, depending on $X-Y$ being positive or negative. I don't understand why.
I'm thinking that since we've got $X \in (0,\infty)$ and $Y \in (0, \infty)$ we should have $Z \in (-\infty, \infty)$.
Any hints?
I just thought of another way: the moment generating function. But there's a catch, we only know that product of MGF's is an MGF of a SUM (not difference). I'm wondering if a ratio of MGF's might be an MGF of a difference?

This is a perfect example to show that indicator functions should not be omitted from the expressions of densities--or rather that, if they are, absurd conclusions may follow.
Here the densities are NOT what you wrote but the functions $f$ and $g$ defined on $\mathbb R$ by $$ f(x)=\lambda\mathrm e^{-\lambda x}\mathbf 1_{x\geqslant0},\qquad g(y)=\mu\mathrm e^{-\mu y}\mathbf 1_{y\geqslant0}. $$ Thus, for every $z$, $$ P[X-Y\leqslant z]=\iint f(x)g(y)\mathbf 1_{x-y\leqslant z}\mathrm dx\mathrm dy=\int_0^{+\infty}\int_0^{+\infty}\lambda\mathrm e^{-\lambda x}\mu\mathrm e^{-\mu y}\mathbf 1_{x-y\leqslant z}\mathrm dx\mathrm dy, $$ that is, $$ P[X-Y\leqslant z]=\int_0^{+\infty}\mu\mathrm e^{-\mu y}\int_0^{(y+z)^+}\lambda\mathrm e^{-\lambda x}\mathrm dx\mathrm dy=\int_0^{+\infty}\mu\mathrm e^{-\mu y}(1-\mathrm e^{-\lambda (y+z)^+})\mathrm dy, $$ and from this point, one should treat separately the cases $z\leqslant0$ and $z\geqslant0$.
A more direct approach is to note that the density $h$ of $X-Y$ is defined by $$ h(z)=\int f(y+z)g(y)\mathrm dy=\int \lambda\mathrm e^{-\lambda (y+z)}\mu\mathrm e^{-\mu y}\mathbf 1_{y+z\geqslant0}\mathbf 1_{y\geqslant0}\mathrm dy, $$ that is, $$ h(z)=\int_{z^-}^{+\infty}\lambda\mathrm e^{-\lambda (y+z)}\mu\mathrm e^{-\mu y}\mathrm dy=\frac{\lambda\mu\mathrm e^{-\lambda z}}{\lambda+\mu}\mathrm e^{-(\lambda+\mu)z^-}, $$ which can be rewritten as the more symmetrical $$ h(z)=\frac{\lambda\mu}{\lambda+\mu}(\mathrm e^{-\lambda z}\mathbf 1_{z\geqslant0}+\mathrm e^{\mu z}\mathbf 1_{z\leqslant0}). $$