Make a first approximation, , and substitute into the rhs of the Kepler equation to get a better approximation. Keep doing this. One finds the following succession of approximations:

phi[0]=M;
for j to N do
phi[j] = subs(phi=rhs(%),rhs(Keq));
subs( f=phi, expansion( subs(phi=f,%), e, j ) );
collect(combine(%,trig),e);
print(%)
od: