Know that the solution to Kepler's equation can in fact be expressed as a series of Bessel functions of the first kind:
$$E=M+2\sum_{k=1}^\infty \frac{\sin(k M)J_k(k\epsilon)}{k}$$
but since the Bessel function is a bit more complicated to use, one is better off using Newton-Raphson or Halley to solve Kepler's equation (particularly convenient due to the properties of the derivatives of the sine and cosine, so you can recycle common subexpressions in numerical evaluation). That leaves you with the problem of starting up the iteration; I personally prefer the starting approximation due to G.R. Smith:
$$E\approx M+\frac{\epsilon\;\sin\;M}{1+\sin\;M-\sin(\epsilon+M)}$$
which is easily polished off subsequently with Halley, Newton-Raphson, or the simple fixed-point iteration mentioned by Christian. the paper shows how one can arrive at this via linear interpolation.
FWIW, the first thing you should have done before asking was to search the archives of the Celestial Mechanics and Dynamical Astronomy journal at Springer's web page or The SAO/NASA Astrophysics Data System; there are quite a number of survey articles on the efficient routes for evaluating the solutions to Kepler's (elliptic) equation, as well as methods for the parabolic ($\epsilon=1$) and hyperbolic ($\epsilon > 1$) cases.
Lastly, as a meta-note since I can't comment, this should be tagged [astronomy]
and/or [celestial-mechanics]