I am entirely self-taught when it comes to Laplace transforms, and I am seeking an independent opinion on my attempt to work out how to arrive at the below expression (note: I am interested specifically in the process, not the just the end result): $$ \begin{align} \mathcal{L}^{-1} \left\{\left(\frac{\lambda}{\lambda+s}\right)^P\frac{\lambda}{s^3} \right\} &= \frac{\lambda^{P+1}}{(P-1)!}\iiint{e^{-\lambda t}t^{P-1} (dt)^{3}}\\ &=\frac{\lambda t^{2}}{2}-Pt+\frac{P(P+1)}{2\lambda}-\frac{e^{\lambda t}}{2\lambda}\left[P(P+1)\sum_{j=0}^{P-1} \frac{(\lambda t)^j}{j!} -2P\sum_{j=0}^{P-2} \frac{(\lambda t)^{j+1}}{j!} +\sum_{j=0}^{P-3} \frac{(\lambda t)^{j+2}}{j!} \right] \end{align}$$ Please note that I disagree with the indefinite integral notation, but it is from the academic papers I've used to put this expression together, specifically: equation 2.13 in this paper and equations 7.4-7.5 in this other paper. These references provide little detail on how the final result is arrived at. Yet I thought I'd offer my best stab at replicating the above result from scratch. Grateful if anyone could advise how far off I am from guessing the correct process that brings about the end result.
Replication attempt strategy
A little bit of context (more in the references I've given at the start): the result I'm trying to verify arises in the context of production-inventory systems, and the inverse Laplace transform I'm referring to is meant to be a cumulative measure of expected "backlog", which is understood as the non-negative difference between cumulative demand and cumulative production over a time interval. The latter is just a given discrete quantity ($P$ in our notation), whereas the former is a Poisson random variable with parameter $\lambda$.
First two (inner) integrals
To start with, I exploit the result that $\mathcal{L}^{-1}\left\{ \left(\frac{\lambda}{s+\lambda}\right)^n \right\}=\lambda^n\frac{t^{n-1}}{(n-1)!}e^{-\lambda t}$ (see e.g. Muth (1977) p.198), and that $\int_{0}^{t}\lambda^n\frac{x^{n-1}}{(n-1)!}e^{-\lambda x}dx =1-\sum_{j=0}^{n-1}\frac{\lambda^j e^{-\lambda t}}{j!}t^{j}$ (see Wagner (1975) p.859) to show the following equivalence: $$ \begin{align} \left(\frac{\lambda}{\lambda+s} \right)^{P}\frac{\lambda}{s^{2}} &= \lambda\frac{1}{s^{2}} \mathcal{L}\left\{\lambda^{P}\frac{t^{P-1}}{(P-1)!}e^{-\lambda t}\right\}\\ &=\mathcal{L}\left\{\lambda\int_{0}^{t}\left[\int_{0}^{t}\lambda^{P}\frac{u^{P-1}}{(P-1)!}e^{-\lambda u} du \right]du\right\}\\ &=\mathcal{L}\left\{\lambda\int_{0}^{t}\left[ 1-\sum_{j=0}^{P-1}\frac{\lambda^j e^{-\lambda u}}{j!}u^{j} \right] du\right\}\\ &=\mathcal{L}\left\{ \lambda t - \lambda\left[\int_{0}^{t}\frac{\lambda^{0}e^{-\lambda u}}{0!}u^{0} du + \int_{0}^{t}\frac{\lambda e^{-\lambda u}}{1!}u du + \dots + \int_{0}^{t}\frac{\lambda^{P-1}e^{-\lambda u}}{(P-1)!}u^{P-1}du \right] \right\} \\ &=\mathcal{L}\left\{ \lambda t - P + e^{-\lambda t}\sum_{i=0}^{P-1}\sum_{j=0}^{i}\frac{(\lambda t)^j}{j!} \right\} \\ &=\mathcal{L}\left\{ \lambda t - P + e^{-\lambda t}\left[\sum_{i=0}^{P-1}P\frac{(\lambda t)^i}{i!}-\sum_{i=0}^{P-1}i\frac{(\lambda t)^i}{i!} \right] \right\} \\ \end{align} $$
The third (outer) integral
Based on the above equivalence (and by re-using some of the results mentioned in the previous section), I conclude that:
$$ \begin{align} \mathcal{L}^{-1} \left\{\left(\frac{\lambda}{\lambda+s}\right)^P\frac{\lambda}{s^3} \right\} &=\int_{0}^{t}\lambda x - P + e^{-\lambda x}\left[\sum_{i=0}^{P-1}P\frac{(\lambda x)^i}{i!}-\sum_{i=0}^{P-1}i\frac{(\lambda x)^i}{i!} \right] dx \\ &=\frac{\lambda t^2}{2}-Pt+P\sum_{i=0}^{P-1}\frac{\lambda^i}{i!}\int_{0}^{t} x^i e^{-\lambda x}dx - \sum_{i=0}^{n-1}i\frac{\lambda^i}{i!} \int_{0}^{t} x^i e^{-\lambda x}dx \\ &= \frac{\lambda t^2}{2}-Pt+\frac{P}{\lambda}\sum_{i=0}^{P-1} \left[-\frac{(\lambda t)^i}{i!}e^{-\lambda t} + \int_{0}^{t} \lambda^i \frac{x^{i-1}}{(i-1)!} e^{-\lambda x}dx \right] - \frac{1}{\lambda}\sum_{i=0}^{P-1} \left[-i\frac{(\lambda t)^i}{i!}e^{-\lambda t} + i\int_{0}^{t} \lambda^i \frac{x^{i-1}}{(i-1)!} e^{-\lambda x}dx \right]\\ &= \frac{\lambda t^2}{2}-Pt-\frac{P}{\lambda}\sum_{i=0}^{P-1}\frac{(\lambda t)^i}{i!}e^{-\lambda t} +\frac{P}{\lambda}\sum_{i=0}^{P-1} \left[1-\sum_{j=0}^{i-1}\frac{(\lambda t)^j}{j!}e^{-\lambda t} \right] + \frac{1}{\lambda}\sum_{i=0}^{P-1}i\frac{(\lambda t)^i}{i!}e^{-\lambda t} - \frac{1}{\lambda}\sum_{i=0}^{P-1} i \left[ 1-\sum_{j=0}^{i-1}\frac{(\lambda t)^j}{j!}e^{-\lambda t} \right] \\ &= \frac{\lambda t^2}{2}-Pt-\frac{P(P+1)}{2\lambda}-\frac{e^{-\lambda t}}{\lambda}\left[P \sum_{i=0}^{P-1}\frac{(\lambda t)^{i}}{i!} - \sum_{i=0}^{P-2}\frac{(\lambda t)^{i+1}}{i!} + P \sum_{i=0}^{P-2}(P-1+i)\frac{(\lambda t)^{i}}{i!}\right] + \frac{1}{\lambda}\sum_{i=0}^{P-1}i \left[\sum_{j=0}^{i-1} \frac{(\lambda t)^j}{j!}e^{-\lambda t} \right] \\ &=\frac{\lambda t^2}{2}-Pt-\frac{P(P+1)}{2\lambda}-\frac{e^{-\lambda t}}{\lambda}\left[P^2\sum_{i=0}^{P-2}\frac{(\lambda t)^{i}}{i!}-(P+1)\sum_{i=0}^{P-3}\frac{(\lambda t)^{i+1}}{i!} +\frac{(\lambda t)^{P-1}}{(P-1)!}\right]+\frac{1}{\lambda}\sum_{i=0}^{P-2}\left[\frac{P(P-1)}{2}-\frac{i(i+1)}{2}\right]\frac{(\lambda t)^i}{i!}e^{-\lambda t}\\ &=\frac{\lambda t^2}{2}-Pt-\frac{P(P+1)}{2\lambda}-\frac{e^{-\lambda t}}{\lambda}\left[\frac{P(P+1)}{2}\sum_{i=0}^{P-2}\frac{(\lambda t)^i}{i!}-\frac{2P}{2}\sum_{i=0}^{P-3}\frac{(\lambda t)^{i+1}}{i!}+\frac{1}{2}\sum_{i=0}^{P-4}\frac{(\lambda t)^{i+2}}{i!}+\frac{(\lambda t)^{P-1}}{(P-1)!} \right] \end{align} $$
This is the closer I've managed to get. Obviously it's an extremely tedious process and I've had to skip some steps. Also, the way I've gotten rid of double summations relied entirely on tabulating them - except for the most obvious ones. I have posted separately about the issue of double summations in this other post.
And yet, my process did not quite yield the same as, say equation 2.13 in the original paper. Namely, in the paper the upper bounds in the summations are shifted, and the last term in the square bracket is omitted (as you can see at the start of this post). Also I've followed the suggestion in the comments, and punched $\mathcal{L}^{-1} \left\{\left(\frac{\lambda}{\lambda+s}\right)^P\frac{\lambda}{s^3} \right\}$ into Wolfram Alpha. Once the solution (given for the most part in terms of Gamma and incomplete Gamma) is developed, I get a similar expression that differs - once again - for the upper bounds in the summations and the last term in the square brackets.
As a self-learner, I'd be grateful if I could get an independent view of what might be wrong (if anything) in the way I've tried to replicate this result. Thanks.