How to compute $\sum\limits_{a=1}^{\infty}\int_0^b\lambda\left(\int_0^{\lambda}e^{-t}t^{a-1}dt\right)d\lambda$ Please suggest an efficient method to compute the following integral
\begin{equation}
I = \sum\limits_{a=1}^{\infty}\int_0^b\lambda\left(\int_0^{\lambda}e^{-t}t^{a-1}dt\right)d\lambda
\end{equation}
One can see that the quantity in the parenthesis is incomplete gamma function. So, one of the methods is to write this gamma function as an infinite series and perform the integration term by term. But since it will create one infinite series inside another, it will be computationally expensive to follow this method.
Can anyone suggest any other faster way?
What would happen if the integral is modified as follows?
\begin{equation}
I = \sum\limits_{a=1}^{\infty}\int_0^b\lambda\left(\int_0^{\lambda}e^{-t}t^{a-1}dt\right) \left(\int_0^{\alpha\lambda}e^{-t_1}t_1^{a-1}dt_1\right) d\lambda
\end{equation}

You can reverse the order of integration. The set you are integrating is
$$\{(t,\lambda): \lambda\in [0,b], t\in [0,\lambda]\}$$
This set can also be expressed as
$$ \{ (t,\lambda):t\in [0,b],\ \lambda\in [t,b] \} $$
Therefore you get
$$I = \sum_{a=1}^\infty \int_0^b\left(\int_t^b\lambda e^{-t}t^{a-1}d\lambda \right)dt = \sum_{a=1}^\infty \int_0^be^{-t}t^{a-1} \left(\int_t^b\lambda d\lambda \right)dt$$
$$ = \sum_{a=1}^\infty \int_0^be^{-t}t^{a-1} \frac{b^2-t^2}{2}dt $$
The rest of the calculation is pretty straight forward.
