Here is another method:
So, you have the product of two infinite series. You want to multiply out:
Those are all of the terms with the power of less than or equal to 3.
Another way to write this:
The substitution method doesn't work for the same reason as with what I first tried. I have used the method given by SlipEternal but the combined sigma notation at the end I don't know how to obtain. I haven't seen subfactorial before so I won't be looking at the last one.
So, applying that to the product of infinite series I gave you, you get the solution I provided (the double summation at the end of post #6).
Btw, Wolframalpha can compute the subfactorial, so if you want to check your answer: link. In the link I offered for subfactorial, it actually states that the exponential generating function is , so plugging in , you get the same expansion I offered, so this is a correct solution.