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 get (excuse me for mimicking Wikipedia)

\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.