(*parametersofmysystem*)ms=1.989*10^(30);m=0.00000003003*{1,5,10,50,100,500,1000,5000,10000};mA=1;au=1;mB=0.0003003;mub=mA*mB/(mA+mB);g=4*Pi^2;ab=10;Lb=mub*Sqrt[g*(mA+mB)*ab];eb=Table[j/500,{i,1},{j,0,500}];a1=60;alpha1=ab/a1;A1=g*mub*m*alpha1^2/(a1*8);L1=m*Sqrt[g*(mA+mB)*a1];psi1=g*mub*m*15*(mA-mB)*alpha1^3/(a1*64*(mA+mB));final=0;leng=Length[m];lbin=Sqrt[1-eb^2];l=Table[0,{i,leng},{j,Length[eb[[1]]]}];For[i=1,i≤leng,i++,l[[i]]=(lbin+(L1[[i]]/Lb))[[1]];];(*Hamiltonian*)Hor[x_,L_,o_,t_]=(-A1[[6]]*(5-3(L-(L1[[6]]*x)/Lb)^2)*(1/x)^3+psi1[[6]]*(Sqrt[(1-(L-(L1[[6]]*x)/Lb)^2)*(1-x^2)]*(1/x)^5)*(7-3*(L-(L1[[6]]*x)/Lb)^2)*Cos[o]);(*Equationsofmotion*)k[x_,L_,o_,t_]=D[Hor[x,L,o,t],x]*(1/L1[[6]]);w[x_,L_,o_,t_]=-D[Hor[x,L,o,t],o]*(1/L1[[6]]);q[x_,lb_,L_,t_]=lb+(L1[[6]]/Lb)*x-L;k[x_,L_,o_,t_]=k[x,L,o,t]/.xx[t];k[x_,L_,o_,t_]=k[x,L,o,t]/.oo[t];w[x_,L_,o_,t_]=w[x,L,o,t]/.xx[t];w[x_,L_,o_,t_]=w[x,L,o,t]/.oo[t];q[x_,lb_,L_,t_]=q[x,lb,L,t]/.xx[t];q[x_,lb_,L_,t_]=q[x,lb,L,t]/.lblb[t];num=301;n1[x_,o_,t_]=With[{L=l[[6,num]]},k[x,L,o,t]];p1[x_,o_,t_]=With[{L=l[[6,num]]},w[x,L,o,t]];ql[x_,lb_,t_]=With[{L=l[[6,num]]},q[x,lb,L,t]];(*Conditionofallowedinitialconditionsoftheangularmomentum*)min1=1-alpha1*(eb[[1,num]]+1);min=Sqrt[1-min1^2];(*Solvingthesystemformultipleinitialconditionsofangularmomentumbutsameinitialconditionofperiapsis*)time=10^8;number=500;initial=Subdivide[min+0.001,0.99999,30];data=Table[0,{i,1},{j,Length[initial]}];tval=Subdivide[0,time,number];xval=Table[0,{i,1},{j,number}];oval=Table[0,{i,1},{j,number}];For[p=1,p≤Length[initial],p++,{n}=NDSolve[{x'[t]==p1[x,o,t],o'[t]==n1[x,o,t],x[0]==initial[[p]],o[0]==Pi},{x,o},{t,0,time}];For[i=2,i<number,i++,xval[[1,i]]=x[tval[[i]]]/.n;oval[[1,i]]=Mod[o[tval[[i]]]/.n,2*Pi];];data[[1,p]]=Transpose@{oval[[1,1;;]],Sqrt[1-xval[[1,1;;]]^2]};];