Chapter 3, Section 3.13, Figure 3.30

For all the starting values we are going to take the same initial y = 200 with initial \(TraditionalForm\`\(x\& . \)\) = 0.

y0=200.0;  xd0=0;

Now we solve the scaled form of Hill's equations (given in Eq.(3.226) and (3.227) for 81 starting values of x. For each we calculate the initial \(TraditionalForm\`\(y\& . \)\) and solve for 1000 time units selecting only those points of the trajectory within ±200 in x and ±10 in y.  A list of points is created for each trajectory.  
N.B. This can take a while!

Do[
x0=i*0.1;
yd0=-1.5*x0;
sol = NDSolve[
{x''[t]-2y'[t] ==  3x[t]*(1-1/(Sqrt[x[t]^2+y[t]^2])^3),
 y''[t]+2x'[t] == -3y[t]/(Sqrt[x[t]^2+y[t]^2])^3,
 x[0] == x0, y[0] == y0, x'[0] == xd0, y'[0] == yd0},
 {x,y}, {t,0,1000}, MaxSteps->10000];
p=
Map[Flatten,
 Table[Evaluate[{-y[t],x[t]} /. sol], {t,0,1000,0.1}]];
pt = Select[p, (Abs[#[[1]]]<=200 && Abs[#[[2]]]<=10)&];
plist[i]=pt,
{i,1,81,1}
]

Because a lot of points are devoted to turning orbits close to y = 0 we select only every tenth point of the first ten trajectories.

Do[
pt=plist[i];
len = Length[pt];
tenth = Table[j,{j,1,len,10}];
If[i<=10, sp=pt[[tenth]], sp=pt];
splist[i]=sp,
{i,1,81}
];

Because we would get identical trajectories starting from the diametrically opposite side (look at the equations of motion) we make use of the results we have to derive plots for two sets of trajectories.

Do[
sp=splist[i];
pp[i]=
ListPlot[sp,
PlotJoined->True,PlotRange->{{-200,200},{-9,9}},
DisplayFunction->Identity];
pp[i+1]=
ListPlot[-sp,
PlotJoined->True,PlotRange->{{-200,200},{-9,9}},
DisplayFunction->Identity],
{i,1,79,2}
];

Now we can display the plots together.

Show[
Table[pp[i],{i,1,80,1}],
 PlotRange->{{-200,200},{-9,9}},
 DisplayFunction->$DisplayFunction,
 AxesOrigin->{-200,-9},Frame->True]