Chapter 10, Section 10.3.4, Table 10.1

Use the definition of LaplaceCoefficient from the DisturbingFunction package.

LaplaceCoefficient[s_, jj_, a_, n_] :=
Module[{outside, j=Abs[jj], cterm, cbterm, csum, i, k},
outside = 2*Product[s+k,{k,0,j-1}]/j!*a^(j-n);
cterm = 1;
cbterm = cterm * Product[j-k, {k,0,n-1}];
For[i=1; csum = cbterm, 
    (cbterm > $MachineEpsilon) || (cbterm == 0),
    i++, 
    cterm = cterm * (s+i-1) * (s+j+i-1)/i/(j+i) * a^2;
    cbterm = cterm * Product[j+2*i-k, {k,0,n-1}];
    csum = csum + cbterm
   ];
Expand[outside*csum]
]

Define a function Resonance that calculates the semi-major axis of the resonance given by the argument arg.

Resonance[gm_, rp_, as_, ns_, cpds_, omds_,
          ms_, es_, is_, j2_, j4_, arg_List, inner_] :=
(*
Determines the resonant location of the given
argument.

 gm   = gravitational constant * mass of planet (km^3/s^2)
 rp   = radius of planet (km)
 as   = semimajor axis of perturbing satellite (km)
 ns   = mean motion of perturbing satellite (deg/d)
 cpds = rate of change of pericentre of satellite (deg/d)
 omds = rate of change of node of satellite (deg/d)
 ms   = mass of satellite in units of planet's mass
 es   = eccentricity of satellite
 is   = inclination of satellite (deg)
 j2   = J2 of planet
 j4   = J4 of planet
 arg  = resonant argument, i.e. the coefficients of
        lambda', lambda, curly pi', curly pi, Omega',
        Omega, given as a list
 inner= denotes if resonance is interior (inner=1)
        or exterior (inner=0)
*)
Module[{ratio=Abs[arg[[2]]/arg[[1]]], fact=60*60*24},
If[inner==1, n0 = ns*ratio, n0 = ns/ratio];
a0 = N[(gm*fact*fact/(n0 Degree)^2)^(1/3)];
n = N[(Sqrt[gm/a^3]/Degree)*fact*
     (1 + (3/4)*j2*(rp/a)^2
      - ((9/32)*j2^2+(15/16)*j4)*(rp/a)^4)];
If[inner==1,
seccpdot = N[ms*n*(a/as)*0.25*(
    2*(a/as)*LaplaceCoefficient[1/2, 0, a/as, 1]+
    (a/as)^2*LaplaceCoefficient[1/2, 0, a/as, 2])],
seccpdot = N[ms*n*(as/a)*0.25*(
    2*(as/a)*LaplaceCoefficient[1/2, 0, as/a, 1]+
    (as/a)^2*LaplaceCoefficient[1/2, 0, as/a, 2])]
];
cpdot =  N[(Sqrt[gm/a^3]/Degree)*fact*
          ( (3/2)*j2*(rp/a)^2
           -(15/4)*j4*(rp/a)^4)] + seccpdot;
If[inner==1,
secomdot = N[-ms*n*(a/as)*0.5*(
    (a/as)*LaplaceCoefficient[3/2, 1, a/as, 0])],
secomdot = N[-ms*n*(as/a)*0.5*(
    (as/a)*LaplaceCoefficient[3/2, 1, as/a, 0])]
];
omdot = -N[(Sqrt[gm/a^3]/Degree)*fact*
          ( (3/2)*j2*(rp/a)^2
        -((9/4)*j2*j2+(15/4)*j4)*(rp/a)^4)] + secomdot;
If[inner==1,
ares = a /. FindRoot[
         arg.{ns, n, cpds, cpdot, omds, omdot}==0,
         {a, a0}],
ares = a /. FindRoot[
         arg.{n, ns, cpdot, cpds, omdot, omds}==0,
         {a, a0}]];
nres = n /. a->ares;
{inner, arg, nres, ares}
]

Define the resonant arguments of interest.

\(resarg\  = \ 
{{3, \(-1\), \(-2\), 0, 0, 
          0}, {3, \(-1\), \(-1\), \(-1\), 0, 
          0}, 
{3, \(-1\), 0, \(-2\), 0, 0}, {3, \(-1\), 0, 
          0, \(-2\), 0}, {3, \(-1\), 0, 0, \(-1\), \(-1\)}, {3, \(-1\), 0, 0, 
          0, \(-2\)}, {2, \(-1\), \(-1\), 0, 0, 0}, {2, \(-1\), 0, \(-1\), 0, 
          0}, 
{4, \(-2\), \(-2\), 0, 0, 
          0}, {4, \(-2\), \(-1\), \(-1\), 0, 0}, {4, \(-2\), 0, \(-2\), 0, 
          0}, {4, \(-2\), 0, 0, \(-2\), 0}, {4, \(-2\), 0, 
          0, \(-1\), \(-1\)}, {4, \(-2\), 0, 0, 0, \(-2\)}, {5, \(-3\), 
          0, \(-2\), 0, 0}, \ {5, \(-3\), \(-1\), \(-1\), 0, 
          0}, {5, \(-3\), \(-2\), 0, 0, 0}, \ {5, \(-3\), 0, 0, 
          0, \(-2\)}, {5, \(-3\), 0, 0, \(-1\), \(-1\)}, \ {5, \(-3\), 0, 
          0, \(-2\), 0}, {8, \(-5\), \(-3\), 0, 0, 
          0}, {8, \(-5\), \(-2\), \(-1\), 0, 0}, {8, \(-5\), \(-1\), \(-2\), 
          0, 0}, {8, \(-5\), 0, \(-3\), 0, 0}, {8, \(-5\), \(-1\), 0, \(-2\), 
          0}, {8, \(-5\), \(-1\), 0, \(-1\), \(-1\)}, {8, \(-5\), \(-1\), 0, 
          0, \(-2\)}, {8, \(-5\), 0, \(-1\), \(-2\), 0}, {8, \(-5\), 
          0, \(-1\), \(-1\), \(-1\)}, {8, \(-5\), 0, \(-1\), 
          0, \(-2\)}}\)

Calculate the mean motion at the resonance and the corresponding resonant semi-major axis.  We are using Mimas as the perturber.

\(results\  = 

      Map[Resonance[3.79312*10^\((7)\), \ 60330, \ 185727.25, \ 
            381.9945087, \ 1.0008, \ \(-0.9995\), \ 6.46*10^\((\(-8\))\), \ 
            0.02014, \ 1.585, \ 
            1.6297*10^\((\(-2\))\), \ \(-9.10\)*10^\((\(-4\))\), \ \ #, 
            1] &, \ resarg]\)

Sort the results by semi-major axis.

Sort[results,  (#2[[4]] > #1[[4]]) &]