Calculate probability density function from moment generating function

Inverting generating functions is not always easy. It does help to state any constraints on the parameters (which are left unstated here by the OP). Anyway, the approach is: First, convert your mgf into a characteristic function (i.e. replace t -> I t):

Next, invert the characteristic function to yield the pdf (using Mathematica here), using an inverse fourier transform. To get a 'neat' result, it appears we need the parameter constraints {a > 0, $\lambda$ > 0}, so I have added these as assumptions:

HeavisideTheta[x] is a function that is 0 for x<0, and 1 for x > 0, so our pdf is:

Finally, a quick check to make sure all is OK: generate the mgf from the pdf to see it is the same as that we started with (here using the Expect function in mathStatica/Mathematica):

... and all is OK :)

Finally, here is a quick plot of the derived pdf, as parameter $\lambda$ varies:

PlotDensity[f /. {a -> 2, \[Lambda] -> {2,3,4}}]

(This explains in detail a suggestion by @Dilip then expands briefly on another one by @bryanis2010, both leading to excellent general approaches. Thus, the aim of this post is to provide a mathematical answer to a question which has none, so far.)

As every rational function, $M_Z$ can be decomposed into partial fractions, here, for every $a\ne1$, $$ M_Z(t)=\frac{u}{\lambda-at}+\frac{v}{\lambda-t}, $$ for some parameters $u$ and $v$ one can readily identify as $u=-a\lambda/(1-a)$ and $v=\lambda/(1-a)$.

Now, one of the easiest to remember generating functions is that, if $X$ is standard exponential, then $$ M_X(t)=E[\mathrm e^{tX}]=\int_0^{+\infty}\mathrm e^{tx}\mathrm e^{-x}\mathrm dx=\frac1{1-t}. $$ Hence, for every positive $\mu$, $$ \frac1{1-\mu t}=E[\mathrm e^{t\mu X}]. $$ Applying this to $\mu=1/\lambda$ and to $\mu=a/\lambda$ yields $$ M_Z(t)=\frac{u}{\lambda}E[\mathrm e^{taX/\lambda}]+\frac{v}{\lambda}E[\mathrm e^{tX/\lambda}]. $$ This indicates that $$ M_Z(t)=\frac{u}{\lambda}\int_0^{+\infty}\mathrm e^{tax/\lambda}\mathrm e^{-x}\mathrm dx+\frac{v}{\lambda}\int_0^{+\infty}\mathrm e^{tx/\lambda}\mathrm e^{-x}\mathrm dx, $$ that is, $$ \int_0^{+\infty}\mathrm e^{tx}f_Z(x)\mathrm dx=-\frac{a}{1-a}\int_0^{+\infty}\mathrm e^{tx}\mathrm e^{-x\lambda/a}\frac{\lambda}a\mathrm dx+\frac1{1-a}\int_0^{+\infty}\mathrm e^{tx}\mathrm e^{-x\lambda}\lambda\mathrm dx, $$ hence, finally, $$ f_Z(x)=\frac{\lambda}{1-a}(\mathrm e^{-x\lambda}-\mathrm e^{-x\lambda/a}). $$ The simplest method to solve the case $a=1$ might be to consider the limit when $a\to1$ of this expression, yielding the correct gamma density $$ f_Z^{(a=1)}(x)=\lambda^2x\mathrm e^{-x\lambda}. $$

Another method, suggested by @bryanis2010 in comments is to note that $M_Z$ is the product of two generating functions, each easily identifiable, namely $$ M_Z(t)=\frac{\lambda}{\lambda-at}\cdot\frac{\lambda}{\lambda-t}=E[\mathrm e^{atX/\lambda}]\cdot E[\mathrm e^{tX/\lambda}], $$ hence $Z$ is distributed as $(aX+Y)/\lambda$ where $Y$ is also standard exponential, independent from $X$. This yields $f_Z$ as the convolution of two well known densities, for the same end result as above, the $a=1$ case (gamma random variable) being particularly direct.