I'm having difficulty moving from the generating function, $h(x)$, expressed in the form of partial fractions (which I believe I've successfully found) to the closed form of the recurrence relation, $e_n$.
The recurrence relation is: $e_0=e_1=1, e_2=2$, and $e_n=3e_{n-1}-3e_{n-2}+e_{n-3}$.
I have calculated the generating function, $h(x)$, as: $h(x)=\frac{1-2x+2x^2}{1-3x+3x^2-x^3}$.
Applying partial fractions, obtain: $h(x)=\frac{2}{1-x}+\frac{2}{(1-x)^2}+\frac{1}{(1-x)^3}$.
From Wolfram, I know that $e_n=\frac{1}{2}(n^2-n+2)$, but I can't see how the coefficients of $x_n$ are extracted from the terms to arrive at this solution (with the exception of the first term of the partial fraction, which clearly becomes $2$).
The generating function to apply here is $$\sum_{n=0}^\infty \binom{m+n}{n}x^n = \frac{1}{(1-x)^{m+1}},$$ which has several recent proofs here.
The $m=0$ case is the one you already knew. Now use the $m=1$ and $m=2$ cases.