I want to propose a (recursive) method to obtain $F_k(x)$ by exploiting a formula well known from the theory of Borel re-summation. In the present context it looks like this:
\begin{equation} \begin{array}{ll} \displaystyle \sum_{n \geq 0} \, \langle E_k \rangle_n \, z^n &\displaystyle = \ \int_0^\infty \, dt \, e^{-t} \sum_{n \geq 0} \, {\langle E_k \rangle_n \over {n!}} \, (zt)^n \\ &\displaystyle = \ \int_0^\infty \, dt \, e^{-t} \, F_k(zt) \\ &\displaystyle = \ \Big( {1 \over z} \Big) \cdot \,\text{$\frak{L}$} \big\{ F_k \big\} \Big({1 \over z} \Big) \quad \text{Laplace transform} \end{array} \end{equation}
I'll illustrate the approach for $k=1$. Let $E_1(x) := \sum_{n \geq 0} \, \langle E_1 \rangle_n \, x^n$. In this case the recurrence relation, as indicated the initial post, is
\begin{equation} \underbrace{n \langle E_1 \rangle_n}_{xE_1'(x) - x} \ = \ \left\{ \begin{array}{l} \displaystyle \underbrace{\langle E_1 \rangle_{n-1}}_{xE_1(x)} \ + \ \underbrace{(n-1) \langle E_1 \rangle_{n-2}}_{x^3E_1'(x) \ + \ x^2E_1(x)} \ + \ \\ \displaystyle \underbrace{\ n^2 \langle E_0 \rangle_{n-1} \ }_{\sum_{n \geq 2} n^2 x^n} \ + \ \underbrace{(n-1)(n^2+1)\langle E_0 \rangle_{n-2}}_{\sum_{n \geq 2} (n-1)(n^2+1) \, x^n} \end{array} \right. \end{equation}
for $n \geq 2$ with the initial values $\langle E_1\rangle_0 = 0$ and $\langle E_1 \rangle_1 = 1$. I've "underbraced" the terms of the recurrence by their respective contributions to the $1$-st order non-homogeneous ODE for the generating function $E_1(x)$ which is, upon simplification,
\begin{equation} E'_1(x) \ + \ {1 \over {x-1}} \, E_1(x) \ = \ {x^3 -x^2 + 5x +1 \over {(1+x)(1-x)^5}} \end{equation}
Mathematica tells me that the general solution is of the form
\begin{equation} \ {\text{const} \over {1-x}} \ + \ {1 \over 4} {x^2+x+2 \over {(1-x)^4}} \ + \ {3 \over 8} {1 \over {(1-x)}} \log \Big( {1-x \over {1+x}} \Big) \end{equation}
The value of the constant term is forced to be $\text{const} = - {1 \over 2}$ since $\langle E_1\rangle_0 = 0$. After performing the requisite inverse transforms we should (excuse me for mimicking Wikipedia) we get
\begin{equation} \begin{array}{ll} \displaystyle F_k(zt) &\displaystyle = \ \sum_{n \geq 0} \, {\langle E_1 \rangle_n \over {n!}} \, (zt)^n \\ &\displaystyle = \ {1 \over {2\pi}} \, \int_{-\pi}^{\pi} \, d\theta \, E_1 \big(zt e^{-i\theta} \big) \exp \big( e^{i\theta} \big) \\ &\displaystyle \ \ \ \ \ \text{(Not sure what this is yet.)} \end{array} \end{equation}
ines.