Chapter 3, Section 3.9, Figure 3.16

Set the mass ratio, μ.

μ = 0.001;
μ1 = 1-μ;
μ2 = μ;

Create a graphics object to represent the locations of both masses and the L4 point.

masses = 
Graphics[
{{PointSize[0.02], Point[{ μ1, 0}]},
 {PointSize[0.04], Point[{-μ2, 0}]},
 {PointSize[0.01], Point[{0.5-μ2,N[Sqrt[3]/2]}]}}];

Define a function ThreeBodyOrbit that solves for the trajectory of the particle in the planar, circular restricted problem.

ThreeBodyOrbit[{x0_, y0_, xd0_, yd0_}, tmax_, 
                ppoints_:100, maxsteps_:1000] :=
Module[{sol},
  sol = NDSolve[
  {x''[t]-2*y'[t] ==  
    x[t] - μ1*(x[t]+μ2)/(Sqrt[(x[t]+μ2)^2+y[t]^2]^3) -
           μ2*(x[t]-μ1)/(Sqrt[(x[t]-μ1)^2+y[t]^2]^3),
   y''[t]+2*x'[t] == 
    y[t] - (μ1/(Sqrt[(x[t]+μ2)^2+y[t]^2]^3) +
            μ2/(Sqrt[(x[t]-μ1)^2+y[t]^2]^3))*y[t],
   x[0]==x0, y[0]==y0, x'[0]==xd0, y'[0]==yd0},
   {x,y}, {t,0,tmax}, MaxSteps->maxsteps];
  ParametricPlot[Evaluate[{x[t],y[t]} /. sol],
                 {t,0,tmax}, 
                 PlotRange->All, 
                 AspectRatio->Automatic,
                 PlotPoints->ppoints]
]

Set up the two separate starting conditions for the tadpole orbits.

tadstart1 = 
{0.0065 + 0.5 - μ2, 0.0065 + N[Sqrt[3]/2], 0., 0.};
tadstart2 = 
{0.008  + 0.5 - μ2, 0.008  + N[Sqrt[3]/2], 0., 0.};

Integrate the first trajectory.

tadpole1 =
ThreeBodyOrbit[
{0.0065 + 0.5 - μ2, 0.0065 + N[Sqrt[3]/2], 0., 0.}, 30Pi]

Integrate the second trajectory.

tadpole2 =
ThreeBodyOrbit[
{0.008 + 0.5 - μ2, 0.008 + N[Sqrt[3]/2], 0., 0.}, 31Pi]

Add the masses and the L4 point to the first trajectory plot.

p1 = Show[tadpole1, masses, Axes->None,
          PlotRange->{{-0.8,1.1},{-0.1,1.1}},
          AspectRatio->Automatic]

Add the masses and the L4 point to the second trajectory plot.

p2 = Show[tadpole2, masses, Axes->None,
          PlotRange->{{-0.8,1.1},{-0.1,1.1}},
          AspectRatio->Automatic]

Display both plots together.

Show[GraphicsArray[{p1,p2}]]