  The Bessel-Neumann expansion (of integer order) of a function
$g:\mathbb{C}\rightarrow\mathbb{C}$ corresponds to representing $g$ as a linear
combination of basis functions $\phi_0,\phi_1,\ldots$, i.e., $g(z)=\sum_{\ell =
0}^\infty w_\ell \phi_\ell(s)$, where $\phi_i(z)=J_i(z)$, $i=0,\ldots$, are the
Bessel functions. In this work, we study an expansion for a more general class
of basis functions. More precisely, we assume that the basis functions satisfy
an infinite dimensional linear ordinary differential equation associated with a
Hessenberg matrix, motivated by the fact that these basis functions occur in
certain iterative methods. A procedure to compute the basis functions as well
as the coefficients is proposed. Theoretical properties of the expansion are
studied. We illustrate that non-standard basis functions can give faster
convergence than the Bessel functions.
