For the textbook explanation of this problem (https://i.stack.imgur.com/DYFdT.jpg), in step 2, they use the gamma function in order to integrate the integral. But this gives them a result of 1/lambda?
In this, I substituted lambda for 5 (random number i picked) to check the integration, which I did by parts twice instead of the gamma integration.
And the answer would end up being 2/lambda^2 ? This difference is huge and I'm really confused and convinced their integration isn't correct. Thanks!
Making the problem more general, the link you give seems to address the problem of $$I_n=\int x^n\, e^{-\lambda x}\, dx$$ (in which $\lambda >0$). You can compute the antiderivative by a series of integration by parts of use the definition of the gamma function to get $$I_n=-x^{n+1} (\lambda x)^{-(n+1)}\,\Gamma (n+1,x \lambda )$$ where appears the incomplete gamma function.
Concerning $$J_n=\int_0^\infty x^n\, e^{-\lambda x}\, dx=\lambda ^{-(n+1)} \,\Gamma (n+1)$$ provided that $\Re(\lambda )>0\land \Re(n)>-1$. If $n$ is an integer, it reduces to $$J_n=\lambda ^{-(n+1)} n!$$
So, in step $2$, if I am not wrong, it should be $$E(x^2)=\lambda \int_0^\infty x^2\, e^{-\lambda x}\, dx=\frac 2 {\lambda^2}$$