Clear["Global`*"]

kozai[m1_,m2_,m3_,x0_,y0_,z0_,vx0_,vy0_,vz0_,r_]:=NDSolve[{{x''[t],y''[t],z''[t]}==-m1 m3(({x[t],y[t],z[t]})/((Sqrt[x[t]^2+y[t]^2+z[t]^2])^3))-m2 m3 (({r Cos[t],r Sin[t],z[t]})/((Sqrt[r^2+z[t]^2])^3)),{x[0],y[0],z[0]}=={x0,y0,z0},{x'[0],y'[0],z'[0]}=={vx0,vy0,vz0}},{x,y,z},{t,0,2000}]

s1=kozai[1,1,0.1,1,0,2,0,0.06,-0.02,5];

plot[s_,n_]:=Table[Show[ParametricPlot3D[
Evaluate[{x[t],y[t],z[t]}/.s],{t,2000.(k-10)/n,2000. k/n},PlotRange->{{-3.,3.},{-3.,3.},{-1.,5.}},Boxed->False,Axes->False],Graphics3D[{PointSize[0.02],Brown,Point[{0,0,0}],Black,Arrow[{{-1.5,-1.5,-0.5},{-1.5,-1.5,1}}]}]],{k,10,n}]

ListAnimate[plot[s1,2000],AnimationRunning->False]

stabi[m1_,m2_,m3_,x0_,y0_,vx0_,vy0_,a_,e_]:=NDSolve[{{x''[t],y''[t]}==-m1 m3(({x[t],y[t]})/((Sqrt[x[t]^2+y[t]^2])^3))-m2 m3 (({a Cos[t/7.1],a (Sqrt[1-e^2]) Sin[t/7.1]})/((Sqrt[a^2-e^2 Sin[t/7.1]^2])^3)),{x[0],y[0]}=={x0,y0},{x'[0],y'[0]}=={vx0,vy0}},{x,y},{t,0,5000}]

s2=stabi[1,0,0.02,1,0,0,0.1415,3,0.6];
s3=stabi[1,0,0.01,1,0,0,0.1,3,0.6];

{Plot[{x[t]/.s2,Cos[t/7.1]},{t,0,100}],Plot[{x[t]/.s3,Cos[t/7.1]},{t,0,100}]}

s4=stabi[1,0.2,0.02,1,0,0,0.1415,3,0.6];
s5=stabi[1,0.2,0.01,1,0,0,0.1,3,0.6];

stabiplot[s1_,s2_,n_]:=Table[Show[ParametricPlot[{
Evaluate[{x[t],y[t]}/.s1],Evaluate[{x[t],y[t]}/.s2]},{t,500.(k-1)/n,500. k/n},PlotRange->{{-3.,3.},{-3.,3.}},Axes->False],Graphics[{PointSize[0.05],Brown,Point[{0,0}]}],ListPlot[{{3 Cos[500.(2 k-1)/(2 n 7.1)],3 (Sqrt[1-0.6^2]) Sin[500.(2 k-1)/(2 n 7.1)]}},Axes->False]],{k,1,n}]

ListAnimate[stabiplot[s4,s5,250]]