When e< 0.5, M = E - e*sin(E) is approximately solved by: E = M + e/(1-e)*sin(M)-(2*e^2)/(1-e^2)/Pi^2*M*(2*Pi-M)*sin(M), This seems to be the most simple three-term solution, which matches E(M) and E'(M) at integer multiples of Pi. The error increases with e, so the approximation becomes unusable after e = 0.5, where the error is still less than 1% of 2*Pi. All planets are well within this bound, but some bound comets are not. Check with the wiki-tabulation at: https://en.wikipedia.org/wiki/Orbital_eccentricity For worst-case planet Mercury with e=0.21, accuracy is better than 99.7%. The approximate E(M) is so simple, I think it would have been discovered hundreds of years ago, but this is the first time I've seen it. Anyone with a Reference? Or how about a simple 99% approximation for e > 0.5? --Brad (* Mma Accuracy Check *) Manipulate[Plot[Times[Subtract[ M + e/(1 - e)*Sin[M] - (2*e^2)/(1 - e^2)/Pi^2*M*(2*Pi - M)* Sin[M], InverseFunction[# - e Sin[#] &][M]], 100/2/Pi], {M, 0, 2 Pi}, PlotRange -> {-1, 1}], {e, 0.001, 0.5}]