A survey of the Lie Symetries of ordinary differential equations.

Section 1.  Infinitesimal versus finite Lie transforms.

Given  X = ξ (x, y)(∂_x)  + η (x, y)(∂_y) ,  find the corresponding change of variables.

First example :

xnew = FindGeneratingFunction[Table[((ϵ^k)/(k!))Nest[a( (∂_x) #)+b((∂_y) #)&,x,k],{k,0,10}],u]/.u->1

ynew = FindGeneratingFunction[Table[((ϵ^k)/(k!))Nest[a( (∂_x) #)+b((∂_y) #)&,y,k],{k,0,10}],u]/.u->1

Another example :

xnew = FindGeneratingFunction[Table[((ϵ^k)/(k!))Nest[a(x (∂_x) #)+b(y (∂_y) #)&,x,k],{k,0,10}],u]/.u->1

ynew = FindGeneratingFunction[Table[((ϵ^k)/(k!))Nest[a(x (∂_x) #)+b(y(∂_y) #)&,y,k],{k,0,10}],u]/.u->1

Another example :

xnew = FindGeneratingFunction[Table[((ϵ^k)/(k!))Nest[a((x^2) (∂_x) #)+b((y^2) (∂_y) #)&,x,k],{k,0,10}],u]/.u->1

ynew = FindGeneratingFunction[Table[((ϵ^k)/(k!))Nest[a((x^2) (∂_x) #)+b((y^2)(∂_y) #)&,y,k],{k,0,10}],u]/.u->1

Another example :

xnew = FindGeneratingFunction[Table[((ϵ^k)/(k!))Nest[(-y (∂_x) #)+(x (∂_y) #)&,x,k],{k,0,10}],u]/.u->1

ynew = FindGeneratingFunction[Table[((ϵ^k)/(k!))Nest[(-y (∂_x) #)+(x(∂_y) #)&,y,k],{k,0,10}],u]/.u->1

Given the change of variables, find the corresponding generator (Same four examples reversed).

{ξ,η} = {(∂_ϵ) (x+a ϵ),(∂_ϵ) (y+b ϵ)}/.ϵ->0

{ξ,η} = {(∂_ϵ) ((E^(a ϵ)) x),(∂_ϵ) ((E^(b ϵ)) y)}/.ϵ->0

{ξ,η} = {(∂_ϵ) ((x/(1-a x ϵ))),(∂_ϵ) ((y/(1-b y ϵ)))}/.ϵ->0

{ξ,η} = {(∂_ϵ) (x Cos[ϵ]-y Sin[ϵ]),(∂_ϵ) (y Cos[ϵ]+x Sin[ϵ])}/.ϵ->0

Section 2.  Prolongating Lie generators.

Generic calculation of (η^((k))) (k=1,2) .

Clear[f,g,η,V];     (*General clearing needed when restarting Section 2 *)

Clear[f,g,η];f[0]:=ξ[x,y[x]];g[0]:=η[x,y[x]];g[k_]:=g[k]=(∂_x) g[k-1]-((∂_({x,k})) y[x])(∂_x) f[0];(Subscript[η, k_]):=g[k];

(Subscript[η, 1])

(Subscript[η, 2])//Simplify

Calculation of the third prolongation (V^((3))) of the generator, X=-y(∂_x) +x(∂_y) .

Clear[f,g,η,V];f[0]:=-y[x];g[0]:=x;g[k_]:=g[k]=(∂_x) g[k-1]-((∂_({x,k})) y[x])(∂_x) f[0];(Subscript[η, k_]):=g[k];(Subscript[V, k_]):=Flatten[Append[{f[0],g[0]},Table[g[i],{i,k}]]];

(* The asked prolongation is dispalyed in vector notation : *)

(Subscript[V, 3])

Section 3.  Symmetry condition for the differential equation : (Subscript[y, xx])=ω(x,y,y')=(((Subscript[y, x])^2)/y)-(y^2).

Clear[f,g,ξ,η,V,symcond,splitting,n,ω,equation,mod,α,β,γ,δ,A,X,λ];     (*General clearing needed when restarting Section 3 *)

Clear[equation,mod];equation=(∂_({x,n})) y[x]-ω[x,y[x],y'[x]];
mod=(∂_({x,n})) y[x]->ω[x,y[x],y'[x]];

Clear[f,g];f[0]:=ξ[x,y[x]];g[0]:=η[x,y[x]];g[k_]:=g[k]=(∂_x) g[k-1]-((∂_({x,k})) y[x])(∂_x) f[0];(Subscript[η, k_]):=g[k];(Subscript[V, k_]):=Flatten[Append[{f[0],g[0]},Table[g[i],{i,k}]]];

Input of the order and of the details of the differential equation :

Clear[n,ω];n=2;ω[x0_,y0_,y1_]:=(y1^2)/y0-(y0^2) ;

(* Note : it is not recommended to use notations like y' and y'' for the derivatives without caution; here we use y0 for y, y1 for y' and so on ... instead. In fact Mathematica forbids y' or y'' to be variables and it does not print any warning when it calculates something like this : 
(∂_y) y'
answer  : y''   !!!! 

y'[x] is of course allowed *)

Symmetry condition : (X^((N-1)))ω=(η^((N)))    mod((y^((N)))=ω)

Clear[symcond];symcond=((Subscript[V, n-1]).Prepend[Table[(∂_((∂_({x,j})) y[x])) ω[x,y[x],y'[x]],{j,0,n-1}],(∂_x) ω[x,y,y']/.{y->y[x],y'->y'[x],y''->y''[x],y'''->y'''[x],y''''->y''''[x]}]-g[n])/.mod/.y[x]->y//Simplify//Expand

Resolution of the symmetry condition :

Clear[splitting];splitting=CoefficientList[symcond,y'[x]]//Simplify//Expand              (* Observe the four monomials and the splitting of the condition in four parts *)

Length[splitting]

Clear[ξ,η];DSolve[splitting[[4]]==0,ξ[x,y],{x,y}]             (* Solve the fourth equation *)

Clear[α,β];ξ[x_,y_]:=Log[y] α[x]+β[x]             (* Take the solution for ξ into account *)

Clear[η];DSolve[splitting[[3]]==0,η[x,y],{x,y}]            (* Solve the third equation *)

Clear[γ,δ];η[x_,y_]:=y γ[x]+y Log[y] δ[x]+y ((Log[y])^2) α'[x]            (* Take the solution for η into account *)

splitting//Simplify//Expand            (* Detection of the distinct monomials in the two remaining equations : y and Log[y] *)

CoefficientList[splitting,{y,Log[y]}]//Flatten             (* Final extra conditions for α(x), β(x), γ(x), δ(x); easily solved at hand, see main text *)

(* Solve at hand the elementary remaining equations : *)

α[x_]:=0;β[x_]:=(Subscript[k, 1])+(Subscript[k, 2]) x;γ[x_]:=-2(Subscript[k, 2]) ;δ[x_]:=0;

{ξ[x,y],η[x,y]}      (* = Final solution : (Subscript[X, 1])= (∂_x) and (Subscript[X, 2])= x(∂_x)-2y(∂_y) *)

splitting              (* Final verification : OK *)

Checking the validity of the symmetry condition (same equation, second version : [(X^((N-1))),A ]=λA ):

Clear[A];A:=((∂_x0) #+y1 (∂_y0) #+(((y1^2)/y0)-(y0^2))(∂_y1) #)&

Clear[X];X:= (∂_x0) #&        (* We check the generator (Subscript[X, 1])=(∂_x) (prolonged once) *)

X[A[f[x0,y0,y1]]]-A[X[f[x0,y0,y1]]]//Simplify//Expand

Clear[λ];λ=Coefficient[X[A[f[x0,y0,y1]]]-A[X[f[x0,y0,y1]]],Derivative[1,0,0][f][x0,y0,y1]]

X[A[f[x0,y0,y1]]]-A[X[f[x0,y0,y1]]]-λ  A[f[x0,y0,y1]]//Simplify

Clear[X];X:=x0 (∂_x0) #-2y0  (∂_y0) #-3y1 (∂_y1)#&        (* We check the generator (Subscript[X, 2])(∂_x) - 2y(∂_y)  (prolonged once) *)

X[A[f[x0,y0,y1]]]-A[X[f[x0,y0,y1]]]//Simplify//Expand

Clear[λ];λ=Coefficient[X[A[f[x0,y0,y1]]]-A[X[f[x0,y0,y1]]],Derivative[1,0,0][f][x0,y0,y1]]

X[A[f[x0,y0,y1]]]-A[X[f[x0,y0,y1]]]-λ  A[f[x0,y0,y1]]//Simplify

Section 4.  Changing the variables in a differential equation of order N.

Clear[x,y];     (*General clearing needed when restarting Section 4 *)

Clear[x,y];x[0]=x[r,s[r]];
y[k_]:=Nest[((Dt[#])/(Dt[x[r,s[r]]]))&,y[r,s[r]],k];

x[r_,s_]:=r Exp[s];           (* The change of variables : s is the new dependant variable *)
y[r_,s_]:=Exp[-s];

y[1]-y[0]^2+y[0]/x[0]//Simplify           (* Input of the differential equation : (Subscript[y, x])-(y^2)+y/x = 0 *)

Special case : variables x and y are simply permuted (x=s,y=r).

Clear[x,y];x[0]=x[r,s[r]];
y[k_]:=Nest[((Dt[#])/(Dt[x[r,s[r]]]))&,y[r,s[r]],k];

x[r_,s_]:=s;
y[r_,s_]:=r;

Table[y[j],{j,4}]/.{s->x,r->y}//Simplify

Section 5.  Lie algebra and commutation table of the group (Subscript[R, 3]).

Clear[(Subscript[X, 1]),(Subscript[X, 2]),(Subscript[X, 3])];(Subscript[X, 1]):=y0 (∂_z0) #-z0  (∂_y0) #&;(Subscript[X, 2]):=z0 (∂_x0) #-x0  (∂_z0) #&;(Subscript[X, 3]):=x0 (∂_y0) #-y0  (∂_x0) #&;

Replace[Table[(Subscript[X, i])[(Subscript[X, j])[f[x0,y0,z0]]]-(Subscript[X, j])[(Subscript[X, i])[f[x0,y0,z0]]],{i,3},{j,3}]//Simplify//Expand,{(Subscript[X, 1])[f[x0,y0,z0]]->"(Subscript[X, 1])",(Subscript[X, 2])[f[x0,y0,z0]]->"(Subscript[X, 2])",(Subscript[X, 3])[f[x0,y0,z0]]->"(Subscript[X, 3])",-(Subscript[X, 1])[f[x0,y0,z0]]->"-(Subscript[X, 1])",-(Subscript[X, 2])[f[x0,y0,z0]]->"-(Subscript[X, 2])",-(Subscript[X, 3])[f[x0,y0,z0]]->"-(Subscript[X, 3])"},{2}]//Transpose//MatrixForm

Section 6.  Quest for first integrals.

Clear[splitting,splitting2,Λ,α,β,γ,δ,A,μ,ν,σ,c1,c2,c3,c4,c5];     (*General clearing needed when restarting Section 6 *)

First method : cocharacteristic equation.  Example : (Subscript[y, xxx])=3(Subscript[yy, x])

Clear[A];A:=((∂_x) #+y1 (∂_y0) #+y2 (∂_y1) #+3y0 y1(∂_y2)#)&

Clear[splitting,Λ];splitting=CoefficientList[A[A[A[Λ[x,y0,y1]]]]-3A[y0 Λ[x,y0,y1]]+3y1 Λ[x,y0,y1]//Simplify,y2]         (*Note that Λ does not depend on y2 *)

Length[splitting]

DSolve[splitting[[4]]==0,Λ[x,y0,y1],{x,y0,y1}]

Clear[α,β,γ];Λ[x_,y0_,y1_]:=α[x,y0]+y1 β[x,y0]+(y1^2) γ[x,y0]

Clear[splitting2];splitting2=CoefficientList[splitting[[3]],y1]//Simplify

Clear[σ];γ[x_,y0_]:=σ[x]

DSolve[Derivative[0,1][β][x,y0]+2 Derivative[1,0][γ][x,y0]==0,β[x,y0],{x,y0}]

Clear[δ];β[x_,y0_]:=δ[x]-2y0 σ'[x]

CoefficientList[splitting//Simplify,y1]

DSolve[18 y0 σ[x]-12 σ''[x]+3 Derivative[0,2][α][x,y0]==0,α[x,y0],{x,y0}]

Clear[μ,ν];α[x_,y0_]:=-(y0^3) σ[x]+μ[x]+y0  ν[x]+2 (y0^2) σ''[x]

CoefficientList[splitting,{y0,y1}]//Flatten//Simplify

(* Solve at hand *)

σ[x_]:=c1;ν[x_]:=c2;δ[x_]:=c3 x+c4; μ[x_]:=c5;c3=0;c4=0;

(* Three possible integrand factors are : 1, y and ((Subscript[y, x])^2)-(y^3) : *)

Λ[x,y0,y1]

splitting//Simplify

Second (purely algebraic) method.  Example : (Subscript[y, xxx])=(3/2)((Subscript[y, xx]^2)/(Subscript[y, x]))

(* Six possibilities are to be explored corresponding to (Subscript[f, i])=-1 (i= 1, 2, ..., 6).  
In each case the solution is underdetermined (containing two free parameters) : each coefficient is a constant or a first integral *)

Solve[{f1+x f2+(x^2) f3+g==0,f4+y f5+(y^2) f6+yx g==0,-yx f2-2 x yx f3+yx f5+2y yx f6+yxx g==0,-2yxx f2-(2yx+4x yxx)f3+yxx f5+(2yx^2+2y yxx)f6+3/2 yxx^2/yx  g==0}/.f1->-1,{f2,f3,f4,f5,f6},{g}]

Solve[{f1+x f2+(x^2) f3+g==0,f4+y f5+(y^2) f6+yx g==0,-yx f2-2 x yx f3+yx f5+2y yx f6+yxx g==0,-2yxx f2-(2yx+4x yxx)f3+yxx f5+(2yx^2+2y yxx)f6+3/2 yxx^2/yx  g==0}/.f2->-1,{f1,f3,f4,f5,f6},{g}]

Solve[{f1+x f2+(x^2) f3+g==0,f4+y f5+(y^2) f6+yx g==0,-yx f2-2 x yx f3+yx f5+2y yx f6+yxx g==0,-2yxx f2-(2yx+4x yxx)f3+yxx f5+(2yx^2+2y yxx)f6+3/2 yxx^2/yx  g==0}/.f3->-1,{f1,f2,f4,f5,f6},{g}]

Solve[{f1+x f2+(x^2) f3+g==0,f4+y f5+(y^2) f6+yx g==0,-yx f2-2 x yx f3+yx f5+2y yx f6+yxx g==0,-2yxx f2-(2yx+4x yxx)f3+yxx f5+(2yx^2+2y yxx)f6+3/2 yxx^2/yx  g==0}/.f4->-1,{f1,f2,f3,f5,f6},{g}]

Solve[{f1+x f2+(x^2) f3+g==0,f4+y f5+(y^2) f6+yx g==0,-yx f2-2 x yx f3+yx f5+2y yx f6+yxx g==0,-2yxx f2-(2yx+4x yxx)f3+yxx f5+(2yx^2+2y yxx)f6+3/2 yxx^2/yx  g==0}/.f5->-1,{f1,f2,f3,f4,f6},{g}]

Solve[{f1+x f2+(x^2) f3+g==0,f4+y f5+(y^2) f6+yx g==0,-yx f2-2 x yx f3+yx f5+2y yx f6+yxx g==0,-2yxx f2-(2yx+4x yxx)f3+yxx f5+(2yx^2+2y yxx)f6+3/2 yxx^2/yx  g==0}/.f6->-1,{f1,f2,f3,f4,f5},{g}]

Section 7.  Dynamical symmetries. Example : (Subscript[y, xxxx])=(((Subscript[y, xxx])+((Subscript[y, xxx])^2))/(Subscript[y, xx]))

Clear[f,g,ξ,η,V,symcond,splitting,n,ω,equation,mod,C1,C2,C3,C4,C5,C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17,A,X,λ];     (*General clearing needed when restarting Section 7 *)

Clear[mod];mod=(∂_({x,n})) y[x]->ω[x,y[x],y'[x],y''[x],y'''[x]];

Clear[f,g,η,V];f[0]:=0;g[0]:=η[x,y[x],y'[x],y''[x]];g[k_]:=g[k]=((∂_x) g[k-1]-((∂_({x,k})) y[x])(∂_x) f[0])/.mod;(Subscript[η, k_]):=g[k];(Subscript[V, k_]):=Flatten[Append[{f[0],g[0]},Table[g[i],{i,k}]]];

Input of the differential equation :

Clear[n,ω];n=4;ω[x0_,y0_,y1_,y2_,y3_]:=((y3+(y3^2) )/y2);

(* Note : it is not recommended to use notations like y' and y'' in place of y1 and y2. In fact Mathematica forbids y' or y'' to be variables and it does not print any warning when you calculate something like this : 
(∂_y) y'
answer  : y''   !!!! *)

Symmetry condition :    (*  Mathematica helps those who help themselves ! *)

Clear[symcond];symcond=((Subscript[V, n-1]).Prepend[Table[(∂_((∂_({x,j})) y[x])) ω[x,y[x],y'[x],y''[x],y'''[x]],{j,0,n-1}]/.mod,((∂_x) ω[x,y,y',y'',y''']/.mod)/.{y->y[x],y'->y'[x],y''->y''[x],y'''->y'''[x]}]-g[n])//Simplify//Expand;

Clear[splitting];splitting=CoefficientList[symcond,y'''[x]]/.{x->x0,y[x]->y0,y'[x]->y1,y''[x]->y2};

Length[splitting]

DSolve[splitting[[5]]==0,η[x0,y0,y1,y2],{x0,y0,y1,y2}]

Clear[C1,C2,C3,C4,C5];η[x0_,y0_,y1_,y2_]:=-y2 C1[x0,y0,y1]+y2 Log[y2] C1[x0,y0,y1]-Log[y2] C2[x0,y0,y1]+C3[x0,y0,y1]+y2 C4[x0,y0,y1]

splitting[[4]]//Simplify

CoefficientList[(y2^4)splitting[[4]],{y2,Log[y2]}]//Flatten

(* All these components must vanish.  This leads to the following : *)

C2[x0_,y0_,y1_]=0;C1[x0_,y0_,y1_]=0;C4[x0_,y0_,y1_]=C5[x0,y0];

splitting[[3]]//Simplify

CoefficientList[(y2^2)splitting[[3]]//Simplify,y2]//Flatten

(* All these components must vanish : start with the last one : *)

DSolve[4 Derivative[0,1][C5][x0,y0]+2 Derivative[0,0,2][C3][x0,y0,y1]==0,C3[x0,y0,y1],{x0,y0,y1}]

Clear[C6,C7,C8,C9,C10,C11,C12,C13,C14,C15,C16,C17];C3[x0_,y0_,y1_]:=C6[x0,y0]+y1 C7[x0,y0]-(y1^2) Derivative[0,1][C5][x0,y0]

(* Consider now the remaining components of test[[3]] : *)

CoefficientList[(y2^2)splitting[[3]]//Simplify,{y1,y2}]//Flatten

(* Solve at hand : *)

C5[x0_,y0_]:=a x0 +b y0+s;C7[x0_,y0_]:=C8[x0]+y0 C9[x0];C6[x0_,y0_]:=-((a (x0^2))/2)+C10[y0]+x0 C11[y0];C9[x0_]:=C12-(1/2) x0 C13''[y0]-(1/4) (x0^2) C11''[y0];C13[y0_]:=C14+y0 C15+C10[y0];C8[x0_]:=-((b (x0^2))/2)+C16+x0 C17-(x0^2) C11'[y0]+(1/4) (x0^2) y0 C11''[y0];

CoefficientList[y2 splitting//Simplify,{y1,y2}]//Flatten

(* Solve at hand : *)

C11[y0_]:=-b y0+c;C10[y0_]:=f y0+g;C12=0;b=0;C17=-f/3;

(* The remaining two first components are identically satisfied : *)

(y2^2)splitting//Simplify

(* Final result : six symmetries are found corresponding to the free parameters, g, c, a, f, C16, s : *)

η[x0,y0,y1,y2]

(* Example of the calculation of the third prolongation of a dynamical symmetry : *)

Clear[mod];mod=(∂_({x,4})) y[x]->((y'''[x]+((y'''[x])^2))/(y''[x]));   (* Each time (y^((N))) appears it must be replaced by ω *)
Clear[f,g,η,V];f[0]:=0;g[0]:=y'[x];g[k_]:=g[k]=(∂_x) g[k-1]-((∂_({x,k})) y[x])(∂_x) f[0]/.mod;(Subscript[η, k_]):=g[k];(Subscript[V, k_]):=Flatten[Append[{f[0],g[0]},Table[g[i],{i,k}]]];

(Subscript[V, 4])//Simplify

(* Verification of the symmetry condition : *)

Clear[equ];equ=(∂_({x,4})) y[x]-((y'''[x]+((y'''[x])^2))/(y''[x]))

(Subscript[V, 4]).{(∂_x) equ,(∂_(y[x])) equ,(∂_(y'[x])) equ,(∂_(y''[x])) equ,(∂_(y'''[x])) equ,(∂_(y''''[x])) equ}//Simplify

Section 8.  The Kepler problem.

Clear[f,g,ξ,η,A,V1,ξ2,η2,V2,ξ3,η3,V3,ξ4,η4,V4,symcond4,splitting,equation];     (*General clearing needed when restarting Section 8 *)

Clear[A];A={1,x'[t],y'[t],-((μ x[t])/(((((x[t])^2)+((y[t])^2)))^(3/2))),-((μ y[t])/(((((x[t])^2)+((y[t])^2)))^(3/2)))}

(* Four symmetries are (see the main text)  : 
X1 and X2 (point symmetries) and X3 and X4 (dynamical).  (Subscript[Vn, 1]) is the first prolongation of Xn in vector notation *)

Clear[V1];(Subscript[V1, 1])={1,0,0,0,0}

Clear[ξ2,η2,V2];ξ2[0]:=0;η2[0]:={-y[t],x[t]};η2[k_]:=η2[k]=(∂_t) η2[k-1]-((∂_({t,k})) y[t])(∂_t) ξ2[0];(Subscript[V2, k_]):=Flatten[Append[{ξ2[0],η2[0]},Table[η2[i],{i,k}]],1]//Simplify//Flatten;(Subscript[V2, 2])

Clear[ξ3,η3,V3];ξ3[0]:=0;η3[0]:={y[t]y'[t],x'[t]y[t]-2x[t]y'[t]};η3[k_]:=η3[k]=(∂_t) η3[k-1]-((∂_({t,k})) y[x])(∂_t) ξ3[0]/.{x''[t]->-μ x[t](((x[t]^2+y[t]^2))^(-3/2)),y''[t]->-μ y[t](((x[t]^2+y[t]^2))^(-3/2))};(Subscript[V3, k_]):=Flatten[Append[{ξ3[0],η3[0]},Table[η3[i],{i,k}]],1]//Simplify//Flatten;(Subscript[V3, 2])

Clear[ξ4,η4,V4];ξ4[0]:=0;η4[0]:={x[t]y'[t]-2x'[t]y[t],x[t]x'[t]};η4[k_]:=η4[k]=(∂_t) η4[k-1]-((∂_({t,k})) y[x])(∂_t) ξ4[0]/.{x''[t]->-μ x[t](((x[t]^2+y[t]^2))^(-3/2)),y''[t]->-μ y[t](((x[t]^2+y[t]^2))^(-3/2))};(Subscript[V4, k_]):=Flatten[Append[{ξ4[0],η4[0]},Table[η4[i],{i,k}]],1]//Simplify//Flatten;(Subscript[V4, 1])

(* Checking the validity of the symmetry condition in the case of X4 : (X4^((2)))[equation]=0 *)

Clear[symcond4];symcond4=Flatten[(Subscript[V4, 2])//Simplify]

Clear[equation];equation={x''[t]+μ x[t](((x[t]^2+y[t]^2))^(-3/2)),y''[t]+μ y[t](((x[t]^2+y[t]^2))^(-3/2))};symcond4.{(∂_t) equ,(∂_(x[t])) equ,(∂_(y[t])) equ,(∂_(x'[t])) equ,(∂_(y'[t])) equ,(∂_(x''[t])) equ,(∂_(y''[t])) equation}//Simplify

(* Search for a linear expression of V1 as a function of V2, V3, V4 and A : *)

Solve[- (Subscript[V1, 1])+f2 (Subscript[V2, 1])+f3 (Subscript[V3, 1])+f4 (Subscript[V4, 1])+g A==0,{f2,f3,f4,g}]//Simplify

(* Here is the linear combination exhibiting first integrals : *)

- (Subscript[V1, 1])+2(En/J) (Subscript[V2, 1])-(Rx/(J^2)) (Subscript[V3, 1])-(Ry/(J^2)) (Subscript[V4, 1])+ A//Simplify

(*  The first integrals : *)

Rx=y[t] x'[t] y'[t]+x[t] ((μ/(Sqrt[((x[t])^2)+((y[t])^2)]))-((y'[t])^2));Ry=y[t] ((μ/(Sqrt[((x[t])^2)+((y[t])^2)]))-((x'[t])^2))+x[t] x'[t] y'[t];J=x[t] y'[t]-y[t] x'[t];En=-(μ/(Sqrt[((x[t])^2)+((y[t])^2)]))+((((x'[t])^2)+((y'[t])^2))/2);

(* They are not completely independant : *)

Rx^2+Ry^2-μ^2-2J^2 En//Simplify

Section 9.  Quest for first integrals in the Lorenz system (β=0, σ=(1/3)).

Clear[A,splitting,splitting2,splitting3,C1,C2,C3,C4,C5,C6,C7,p0,p1,p2,α,β,γ,Λ];     (*General clearing needed when restarting Section 3 *)

(* The method of cocharacteristic is used *)

Clear[A];A:=((∂_t) #+x1 (∂_x0) #+x2 (∂_x1) #+((x1-4/3x0)x2/x0+4/3(x1^2)/x0-(x0^2) x1-1/3(x0^3))(∂_x2)#)&

(* Cocharacteristic equation : *)

Clear[splitting];splitting=CoefficientList[A[A[A[Λ[t,x0,x1]]]]+A[A[(x1/x0-4/3) Λ[t,x0,x1]]]-A[(x2/x0+8/3x1/x0-(x0^2))Λ[t,x0,x1]]-(x1 x2/(x0^2)+4/3(x1^2)/(x0^2)+2x0 x1+(x0^2)) Λ[t,x0,x1]//Simplify,x2]         (*Note that Λ does not depend on y2 *)

Length[splitting]

DSolve[splitting[[4]]==0,Λ[t,x0,x1],{t,x0,x1}]

Clear[Λ,α,β,γ];Λ[t_,x0_,x1_]:=α[t,x0]+x1 β[t,x0]+(x1^2) γ[t,x0]

Clear[splitting2];splitting2=CoefficientList[splitting//Simplify,x1]//Flatten

Length[splitting2]

DSolve[splitting2[[12]]==0,γ[t,x0],{t,x0}]

Clear[C1];γ[t_,x0_]:=((C1[t])/x0)

DSolve[splitting2[[11]]==0,β[t,x0],{t,x0}]

Clear[C2];β[t_,x0_]:=((C2[t])/(x0^(2/3)))+(1/3) (16 C1[t]-9 C1'[t])

{splitting2[[11]],splitting2[[10]]}//Simplify

Clear[splitting3];splitting3=CoefficientList[(x0^(8/3))splitting2[[9]],(x0^(2/3))]//Simplify

DSolve[-((32 C1[t])/3)+8 C1'[t]==0,C1[t],t]

Clear[C3];C2[t_]:=0;C1[t_]:=(E^(4 t/3)) C3;

CoefficientList[(x0^2) splitting//Simplify,x1]//Flatten

DSolve[-6 C3 (E^(4 t/3)) (x0^3)-3 α[t,x0]+3 x0 Derivative[0,1][α][t,x0]+3 (x0^2) Derivative[0,2][α][t,x0]==0,α[t,x0],{t,x0}]

Clear[C4,C5];α[t_,x0_]:=(1/4) C3 (E^(4 t/3)) (x0^3)+(((1+(x0^2)) C4[t])/(2 x0))+((I (-1+(x0^2)) C5[t])/(2 x0))

CoefficientList[(x0^3) splitting//Simplify,{x0,x1}]//Flatten       (* Solve at hand *)

C3=0;

Clear[C6];C4[t_]:=I C5[t]+C6

DSolve[2 C6+4 I C5[t]-3 I C5'[t]==0,C5[t],t]

Clear[C7];C5[t_]:=((I C6)/2)+(E^(4 t/3)) C7

C6=0;

splitting//Simplify

(* This achieves the resolution of the cocharacteristic equation *)

(*  Here is the only possible integrand factor : *)

Λ[t,x0,x1]//Simplify

Clear[p0,p1,p2];p2[t_,x0_,x1_]:=(E^(4 t/3)) x0;p1[t_,x0_,x1_]:=-A[p2[t,x0,x1]]-(x1/x0-4/3)(E^(4 t/3)) x0  ;p0[t_,x0_,x1_]:=-A[p1[t,x0,x1]] -(x2/x0+8/3x1/x0-(x0^2))(E^(4 t/3)) x0  ;

(* Combining with (Subscript[X, 1]) or (Subscript[X, 2])  results in the same first integral *)

-x1 p0[t,x0,x1]+ p1[t,x0,x1]A[-x1]+p2[t,x0,x1]A[A[-x1]]//Simplify

(E^(4 t/3))x1 (x0 x2-x1^2+1/4x0^4)p0[t,x0,x1]+ p1[t,x0,x1]A[(E^(4 t/3))x1 (x0 x2-x1^2+1/4x0^4)]+p2[t,x0,x1]A[A[(E^(4 t/3))x1 (x0 x2-x1^2+1/4x0^4)]]//Simplify