Sum of two gamma/Erlang random variables $\Gamma(m,\lambda)$ and $\Gamma(n, \mu)$ with integer numbers $m \neq n, \lambda \neq \mu$ The gamma distribution with parameters $m &gt; 0$ and $\lambda &gt; 0$ (denoted $\Gamma(m, \lambda)$) has density function
$$f(x) = \frac{\lambda e^{-\lambda x} (\lambda x)^{m - 1}}{\Gamma(m)}, x &gt; 0$$
Given two independent gamma random variables $X = \Gamma(m,\lambda)$ and $Y = \Gamma(n, \mu)$ with integer numbers $m \neq n, \lambda \neq \mu$, what is the density function of their sum $X + Y = \Gamma(m,\lambda) + \Gamma(n,\mu)$?
Notice that both $X$ and $Y$ are also Erlang distribution since $m,n$ are positive integers.
My attempt
First, I searched for well-known results about gamma distribution and I got:
(1) If $\lambda = \mu$, the sum random variable is a Gamma distribution $\sim \Gamma(m+n, \lambda)$ (See Math.SE).
(2) $\Gamma(m, \lambda)$ (or $\Gamma(n, \mu)$) is the sum of $m$ (or $n$) independent exponential random variables each having rate $\lambda$ (or $\mu$). The hypoexponential distribution is related to the sum of independent exponential random variables. However, it require all the rates distinct.
(3) This site is devoted to the problem of sums of gamma random variables. In section 3.1, it claims that if $m$ and $n$ are integer numbers (which is my case), the density  function can be expressed in terms of elementary functions (proved in section 3.4). The answer is likely buried under a haystack of formulas (however, I failed to find it; you are recommended to have a try).
Then, I try to calculate it:
$$f_{X+Y}(a) = \int_{0}^{a} f_{X}(a-y) f_{Y}(y) dy \\ = \int_{0}^{a} \frac{\lambda e^{-\lambda (a-y)} (\lambda (a-y))^{m-1}}{\Gamma(m)} \frac{\mu e^{-\mu y} (\mu y)^{n-1}}{\Gamma(n)} dy \\ = e^{-\lambda a} \frac{\lambda^m \mu^n}{\Gamma(m) \Gamma(n)} \int_{0}^{a} e^{(\lambda - \mu) y} (a-y)^{m-1} y^{n-1} dy$$
Here, I am stuck with the integral and gain nothing ... Therefore,


How to compute the density function of $\Gamma(m,\lambda) + \Gamma(n,\mu)$ with integer numbers $m \neq n, \lambda \neq \mu$?

Added: The answers assuming $m = n$ ($\lambda \neq \mu$) are also appreciated.
A closed form expression is provided in the following paper.

SV Amari, RB Misra, Closed-form expressions for distribution of sum of exponential random variables, IEEE Transactions on Reliability, 46 (4), 519-522.
see eqn 3.383(1) of book tables of integrals Gradshteyn and Ryzhik, it reads

$$\int\limits_0^u x^{\nu-1}(u-x)^{\mu-1}e^{\beta x}\,\mathrm{d}x=
B(\mu,\nu)\overline{u}^{\mu+\nu-1} {}_1F_1(\nu;\mu+\nu;\beta u)\\
[\operatorname{Re}\mu&gt;0,\,\operatorname{Re}\nu&gt;0].$$