Searching for one-parameter Lie groups of point transformations of Newton's equations of motion of 2D Kepler problem

(If the description is in red, a modification is required)

Clear["Global`*"]

system of PDEs and infinitesimals of transformation group

Give variables and PDEs

(* Independent variables*)
IndepVar={t};
(* Dependent variables*)
DepVar={x,y};
(* The highest order of derivatives *)
order = 2;
(* PDE system *)
r=Sqrt[x[t]^2+y[t]^2];
PDEs = {D[x[t],t,t]+M*x[t]/r^3,D[y[t],t,t]+M*y[t]/r^3}

Give expressions to substitute for, usually the highest derivatives

subs = {D[x[t],t,t],D[y[t],t,t]}

Give infinitesimals of the one-parameter Lie group of point transformations

(* general ansatz *)
{ξ[t] = Ξ[t,x[t],y[t]],η[x][{}] =α[t,x[t],y[t]],η[y][{}] =β[t,x[t],y[t]]}

(* second ansatz *)
{ξ[t] = δx[t]x[t]+δy[t]y[t]+δ0[t],η[x][{}] =αy[t,x[t]]y[t]+α0[t,x[t]],η[y][{}] =βx[t,y[t]]x[t]+β0[t,y[t]]}

(* third ansatz *)
αy[t,x[t]]=αy0[t]+(δy)'[t]x[t];
βx[t,y[t]]=βx0[t]+(δx)'[t]y[t];
α0[t,x[t]]=α00[t]+α01[t]x[t]+(δx)'[t]x[t]^2;
β0[t,y[t]]=β00[t]+β01[t]y[t]+(δy)'[t]y[t]^2;
{ξ[t] = δx[t]x[t]+δy[t]y[t]+δ0[t],η[x][{}] =αy[t,x[t]]y[t]+α0[t,x[t]],η[y][{}] =βx[t,y[t]]x[t]+β0[t,y[t]]}

(* forth ansatz *)
{ξ[t] = c0 t^2+c1 t+c2,η[x][{}] =c3 y[t]+c0 t x[t]+c4 x[t],η[y][{}] =-c3 x[t]+c0 t y[t]+c5 y[t]}

Determination of infinitesimals of the extension of the point group and replacing derivatives by auxiliary variables in both infinitesimals and PDEs

Store all combinations of independent variables with respect to which one has to differentiate to the array der, including zeroth derivative as the empty set

der=Union[Subsets[Sort[Flatten[Table[IndepVar,{i,1,order}]]],order]]

Calculate infinitesimals and store them to η[f][der] where f stands for dependent variables and der for elements of the array der defined above

ivar = Delete[IndepVar,0]; (* get rid of {} *)
Do[ (* loop over dependent variables *)
Do[ (* loop over all derivatives *)
subder=j[[1;;Length[j]-1]];
varder=Last[j];
η[f][j]=D[η[f][subder],varder]-
Sum[D[ξ[k],varder] D[f[ivar],Delete[subder,0],k],{k,IndepVar}];
(* Print["der = ",subder,",η[",f,"][",j,"] = ",η[f][j]] *),
{j,der[[2;;]]} (* skip zeroth derivative, element {} *)
],
{f,DepVar}
]
η[Last[DepVar]][Last[der]]

Substitute auxiliary arrays like u[i] for derivatives in both PDEs and infinitesimals

Do[
p = 0;
Do[
If[Length[j]==0,
(* substitution for dependent variables *)
PDEs = PDEs/.{f[ivar]-> f[p]};
subs = subs/.{f[ivar]-> f[p]};
Do[Do[η[ff][jj]=η[ff][jj]/.{f[ivar]-> f[p]},{jj,der}],{ff,DepVar}],
(* substitution for their derivatives *)
PDEs = PDEs/.{D[f[ivar],Delete[j,0]]-> f[p]};
subs = subs/.{D[f[ivar],Delete[j,0]]-> f[p]};
Do[Do[η[ff][jj]=η[ff][jj]/.{D[f[ivar],Delete[j,0]]-> f[p]},{jj,der}],{ff,DepVar}]
];
p++;,
{j,der}
],
{f,DepVar}
]
{PDEs, subs}

Apply the extension of the infinitesimal operator to the PDEs

Solve the PDEs to obtain an expression to substitute

sol = Solve[PDEs==0, subs]

Check if the application of the extension of the infinitesimal operator to the PDEs gives zero when the condition PDEs = 0 is fulfilled

Do[
p = 0;
zero[i]=Sum[ξ[k] D[PDEs[[i]],k],{k,IndepVar}];
Do[
zero[i] += Sum[η[f][j] D[PDEs[[i]],f[p]],{f,DepVar}];
p++,
{j,der}
];
zero[i]=zero[i]/.sol(* //FullSimplify *),
{i,Length[PDEs]}
]
Table[zero[i],{i,1,Length[PDEs]}]

If it is not zero then find determining equations for point symmetry (give expressions which can have arbitrary values, usually derivatives of y(x))

(* for the first and second ansatz *) (* eqs = Flatten[CoefficientList[Table[zero[i]*(((((x[0])^2)+((y[0])^2)))^(5/2))//Simplify,{i,1,Length[PDEs]}],{x[1],y[1]}]];*)
(* for the third ansatz *)  (* eqs = Flatten[CoefficientList[Table[zero[i]*(((((x[0])^2)+((y[0])^2)))^(5/2))//Simplify,{i,1,Length[PDEs]}],{x[1],y[1],x[0],y[0]}]];Column[Union[Table[eqs[[i]],{i,1,Length[eqs]}]]] *)
(* for the forth ansatz *)  eqs = Flatten[CoefficientList[Table[zero[i]*(((((x[0])^2)+((y[0])^2)))^(5/2))//Simplify,{i,1,Length[PDEs]}],{x[1],y[1],x[0],y[0],t}]];Column[Union[Table[eqs[[i]],{i,1,Length[eqs]}]]]

Try to solve the determining equations

DSolve[{},{},{}]