  An algorithm for constructing a $J$-orthogonal basis of the extended Krylov
subspace $\mathcal{K}_{r,s}=\operatorname{range}\{u,Hu, H^2u,$ $ \ldots, $
$H^{2r-1}u, H^{-1}u, H^{-2}u, \ldots, H^{-2s}u\},$ where $H \in \mathbb{R}^{2n
\times 2n}$ is a large (and sparse) Hamiltonian matrix is derived (for $r =
s+1$ or $r=s$). Surprisingly, this allows for short recurrences involving at
most five previously generated basis vectors. Projecting $H$ onto the subspace
$\mathcal{K}_{r,s}$ yields a small Hamiltonian matrix. The resulting HEKS
algorithm may be used in order to approximate $f(H)u$ where $f$ is a function
which maps the Hamiltonian matrix $H$ to, e.g., a (skew-)Hamiltonian or
symplectic matrix. Numerical experiments illustrate that approximating $f(H)u$
with the HEKS algorithm is competitive for some functions compared to the use
of other (structure-preserving) Krylov subspace methods.
