  A class of bivariate infinite series solutions of the elliptic and hyperbolic
Kepler equations is described, adding to the handful of 1-D series that have
been found throughout the centuries. This result is based on an iterative
procedure for the analytical computation of all the higher-order partial
derivatives of the eccentric anomaly with respect to the eccentricity $e$ and
mean anomaly $M$ in a given base point $(e_c,M_c)$ of the $(e,M)$ plane.
Explicit examples of such bivariate infinite series are provided, corresponding
to different choices of $(e_c,M_c)$, and their convergence is studied
numerically. In particular, the polynomials that are obtained by truncating the
infinite series up to the fifth degree reach high levels of accuracy in
significantly large regions of the parameter space $(e,M)$. Besides their
theoretical interest, these series can be used for designing 2-D spline
numerical algorithms for efficiently solving Kepler's equations for all values
of the eccentricity and mean anomaly.
