One way to understand the calculation is to recall that for a gamma distribution with shape $\alpha$ and scale $\beta$, $$f_X(x) = \frac{x^{\alpha-1} e^{-x/\beta}}{\beta^\alpha \Gamma(\alpha)}, \quad x > 0.$$ The denominator, being independent of $x$, suggests that $1/(\beta^\alpha \Gamma(\alpha))$ is the required multiplicative factor for the density such that $$\int_{x=0}^\infty f_X(x) \, dx = 1.$$ In other words, $$\beta^\alpha \Gamma(\alpha) = \int_{x=0}^\infty x^{\alpha - 1} e^{-x/\beta} \, dx.$$ This holds true for any $\alpha, \beta > 0$. Now, with this in mind, $$\operatorname{E}[X^\nu] = \int_{x=0}^\infty x^\nu f_X(x) \, dx = \int_{x=0}^\infty \frac{x^{\nu + \alpha - 1} e^{-x/\beta}}{\beta^\alpha \Gamma(\alpha)} \, dx = \frac{\beta^{\nu + \alpha} \Gamma(\nu + \alpha)}{\beta^\alpha \Gamma(\alpha)}\int_{x=0}^\infty \frac{x^{\nu + \alpha - 1} e^{-x/\beta}}{\beta^{\nu + \alpha}\Gamma(\nu + \alpha)} \, dx.$$ This is what the provided solution does, but the motivation for doing so should now be clear, because now the integrand represents a gamma density with shape parameter $\nu + \alpha$, and rate $\beta$. Therefore, its integral over its support is also $1$, provided $\nu + \alpha > 0$. It follows that $$\operatorname{E}[X^\nu] = \frac{\beta^\nu \Gamma(\nu + \alpha)}{\Gamma(\alpha)}, \quad \nu > -\alpha,$$ as claimed.
Don't make this so hard for yourself. Simply compute the kernels rather than explicitly integrating. I will change your notation because $\bar X$ is typically used for the sample mean $$\bar X = \frac{1}{n} \sum_{i=1}^n X_i$$ rather than the sample total. Then
$$p(\lambda \mid X, \alpha, \beta) \propto \lambda^n e^{-\lambda n \bar X} \lambda^{\alpha - 1} e^{-\beta \lambda} = \lambda^{n + \alpha - 1} e^{-(n\bar X + \beta)\lambda}$$ which is the kernel of a gamma density with posterior shape hyperparameter $\alpha^* = n + \alpha$ and rate hyperparameter $\beta^* = n \bar X + \beta$, which agrees with your computation (keeping in mind my $\bar X$ differs from yours by a factor of $1/n$). In performing the computation, we discarded any factors that were not functions of $\lambda$. You inadvertently did this by ignoring the fact that your computation resulted in a posterior likelihood which does not integrate to unity; i.e., your expression for $p(\lambda \mid X, \alpha, \beta)$ is not a proper density since the constant terms with respect to $\lambda$ are not the required normalizing factors for a gamma density with shape $\alpha^*$ and rate $\beta^*$.
Next, if we know that the original gamma prior has a normalizing factor of $\beta^\alpha/\Gamma(\alpha)$ since $$p(\lambda, \mid \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} K(\lambda \mid \alpha, \beta)$$ where $K$ is the kernel, and the posterior is gamma with hyperparameters $\alpha^*$, $\beta^*$, it immediately follows that the marginal likelihood is the ratio of the prior normalizing factor divided by the posterior normalizing factor; i.e., $$p(X \mid \alpha, \beta) = \frac{\beta^\alpha/\Gamma(\alpha)}{(\beta^*)^{\alpha^*}/\Gamma(\alpha^*)} = \frac{\beta^\alpha \Gamma(n+\alpha)}{(n\bar X + \beta)^{n + \alpha} \Gamma(\alpha)},$$ because of Bayes' rule: $$p(\lambda \mid X, \alpha, \beta) = \frac{p(X \mid \lambda)p(\lambda \mid \alpha,\beta)}{p(X \mid \alpha, \beta)}.$$ No integration is required. It is important to note that $p(X \mid \alpha, \beta)$ is multivariate with respect to the sample $X = (X_1, \ldots, X_n)$, thus is not itself gamma distributed.
For the posterior predictive distribution, we apply the same principles as described above. First, by Bayes' rule, $$p(x \mid X , \alpha, \beta) = \frac{p(x \mid \lambda) p(\lambda \mid X, \alpha, \beta)}{p(\lambda \mid X, x, \alpha, \beta)}$$ where the denominator is the posterior given the sample $X$ and the new observation $x$, and the numerator is the likelihood; thus the RHS is again the ratio of normalizing factors, but this time the normalizing factors correspond to the posterior density in the numerator, and the posterior plus new observation in the denominator: $$p(x \mid X, \alpha, \beta) = \frac{(\beta^*)^{\alpha^*}/\Gamma(\alpha^*)}{(\beta')^{\alpha'}/\Gamma(\alpha')},$$ where $$\alpha' = \alpha^* + 1 = n+\alpha + 1,$$ and $$\beta' = \beta^* + x = n\bar X + \beta + x.$$ Note that the posterior predictive is a univariate density in the new observation $x$, hence it is instructive to consider its kernel: $$p(x \mid X, \alpha, \beta) \propto \frac{1}{(\beta')^{\alpha'}} = (n \bar X + \beta + x)^{-\alpha'}$$ which is proportional to a Pareto (Type II) density with minimum value parameter $0$, scale parameter $\beta^* = n \bar X + \beta$ and shape parameter $\alpha^* = n + \alpha$.
Best Answer
It is very very easy...just to be solved in a couple of passages.
Observe that
$$P(a\leq X_{\alpha, \beta} \leq b)= \int_{a}^{b}\underbrace{\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}}_{=f}\times \underbrace{e^{-x \beta} }_{=g'}dx=f\times g\Bigg]_a^b -\int_a^b f' \times g=$$
$$=\frac{\beta^{\alpha}}{\Gamma(\alpha)}x^{\alpha-1}\cdot \frac{e^{-x\beta}}{-\beta}\Bigg]_a^b+\int_a^b \frac{\beta^{\alpha}(\alpha-1)}{\Gamma(\alpha)}x^{\alpha-2}\cdot \frac{e^{-x\beta}}{\beta}dx $$
Simplify the expression and get immediately your proof
EDIT: it results to me
$$\frac{1}{\beta}[f_{X_{\alpha,\beta}}(a)-f_{X_{\alpha,\beta}}(b) ]+P(a<X_{\alpha-1,\beta}<b)$$
Can it be a typo in the posted solution which shows $\beta$ instead of $\frac{1}{\beta}$?