Chapter 2, Section 2.5

Derive series for the first five Bessel functions of the first kind using Mathematica's BesselJ function. (cf. Eq.(2.79))

Table[
 {s,Series[BesselJ[s,x],{x,0,5}]},{s,5}] // TableForm

Use the definition given in the text to derive a series for the eccentric anomaly, E, in terms of the mean anomaly, M, and up to the fifth power of the eccentricity, e. (cf. Eq.(2.80))

ea=Collect[
    M+Expand[
      2*Normal[
       Series[
        Sum[(1/s)BesselJ[s, s e] Sin[s M],{s,1,5}],
       {e,0,5}]]], e]

Use the definition given in the text to derive a series for r/a in terms of the mean anomaly, M, and up to the fifth power of the eccentricity, e.  Because of the derivatives with respect to e we need to ensure we take enough terms in the series so first we create a temporary series to sixth order ...

troa = 
Expand[1+e^2/2-2e*
Series[
Sum[(1/s^2)*
D[BesselJ[s, s e],e]Cos[s*M],{s,1,6}],{e,0,6}]]

... and then we truncate this to fourth order.   (cf. Eq.(2.80))

roa = Collect[Expand[Series[Normal[troa],{e,0,4}]], e]

Now we create a similar temporary series for Cos[E] using the formula in the text ...

tcosea = 
Normal[Series[-e/2 + 2*
Sum[(1/s^2)*D[BesselJ[s, s*e],e]*Cos[s*M],{s,1,6}],
         {e,0,6}]]

... and truncate it to fourth order to get the series.  (cf. Eq.(2.82))

cosea = Collect[Normal[Series[tcosea,{e,0,4}]],e]

Use the series for r/a to derive a series for (((a/r))^3).  (cf. Eq.(2.83))

aor3 = 
Collect[
TrigReduce[Normal[
Series[1/(Expand[1+e^2/2-2e*Series[
Sum[(1/s^2)D[BesselJ[s, s e],e]Cos[s M],{s,1,6}],
{e,0,6}]])^3,{e,0,4}]]],e]

Use the formula in the text to derive a series for Sin[f] where f is the true anomaly.  
(cf. Eq.(2.84))

sf = Normal[Series[ExpandAll[Series[2 Sqrt[1-e^2]*
       Sum[(1/s)D[BesselJ[s, s e],e]Sin[s M],{s,1,6}],
       {e,0,6}]],{e,0,4}]]

Now derive a series for Cos[f].  (cf. Eq.(2.85))

cf = Collect[ExpandAll[
       Series[-e+Series[2 (1-e^2)/e *
       Sum[BesselJ[s, s e]*Cos[s M],{s,1,6}],
       {e,0,6}],{e,0,4}]],e]

We can derive a series for f by first calculating d f from derivatives of the series for E derived above.

df=
Collect[TrigReduce[
Normal[Series[Sqrt[1-e^2]*(D[ea, M])^2, {e,0,5}]]],e]

Integrate the resulting series term by term to get the series for f.  (cf. Eq.(2.88))

f = Collect[Integrate[df, M], e]

The equation of centre is just the difference f–M so just take the series for each and subtract.  (cf. Eq.(2.88))

eqc=f-M