In this example we use the package IntroToSymmetry.m to work out the determining equations of the nonlinear second-order ODE 

(Subscript[Y, xx])+x(Subscript[Y, x])+x(Subscript[Y, x]^3)-Y-YYx2=0.

Clear all symbols in the current context.

ClearAll[Evaluate[Context[]<>"*"]]

First read in the package which is located in User Home Folder/Library/Mathematica/Applications/SymmetryAnalysis.

Needs["SymmetryAnalysis`IntroToSymmetry`"]

Enter the input equation as a string. Don't include the ==0 at the end.

inputequation=
"D[y[x],x,x]+x*D[y[x],x]+
x*D[y[x],x]^3-y[x]-y[x]*D[y[x],x]^2";

The function y is a solution of the equation and this constraint must be applied in the form of a rule to the invariance condition. Be sure to check signs. Note that the rule is carefully set up to replace one isolated term. For example, the following rule would give incorrect results.

x*D[y[x],x]->-x*D[y[x],x]^3+y[x]+y[x]*D[y[x],x]^2-D[y[x],x,x]

The problem is that this would only replace the particular product x*D[y[x],x] wherever it appears. This is not as general as replacing the first derivative everywhere. As another example, the rule,

y[x]->x*D[y[x],x]+x*D[y[x],x]^3-y[x]*D[y[x],x]^2+D[y[x],x,x]

will cause the FindDeterminingEquations function to fail to reach completion. The problem here is that y[x] is one of the variables that the infinitesimals are assumed to depend on and the software will attempt to replace y[x] with derivatives that are dependent variables in this context. Generally speaking the idea is to isolate and replace one of the higher derivative terms in the equation. In this case we choose the second derivative term.

rulesarray=
{"D[y[x],x,x]->-x*D[y[x],x]-x*D[y[x],x]^3+y[x]+y[x]*D[y[x],x]^2"};

Enter the list of independent variables.

independentvariables={"x"};

Enter the list of dependent variables.

dependentvariables={"y"};

Enter the list of function names or constant names that must be preserved when the equation is rewritten in terms of generic (x1,y1) variables. In this case there are no constants or function names and so a null list is entered.

frozennames={""};

Enter the maximum derivative order of the equation.

p=2;

The maximum derivative order that the infinitesimals are assumed to depend on is specified by the input parameter r. This parameter is only nonzero when the user is looking for Lie contact groups or Lie-Backlund groups. For the usual case where one is searching for point groups set r=0.

r=0;

When searching for Lie-Backlund groups (r=1 or greater) one can without loss of generality, leave the independent variables untransformed. The corresponding infinitesimals (the xse's) are set to zero by setting xseon=0. If one is searching for point groups then set  xseon=1. The choice xseon=1 is also an option when looking for Lie-Backlund groups and this can be useful when looking for contact symmetries for example.

xseon=1;

When searching for Lie-Backlund groups it is necessary to differentiate the input equation r times producing derivatives of order p+r. These higher order differential consequences are appended to the set of rules applied to the invariance condition. This process is carried out automatically when internalrules=1. For point groups the equation or equation system is the only rule or set of rules needed and one puts internalrules=0.

internalrules=0;

Now work out the determining equations of the Lie point group that leaves the equation invariant. The output is available as a table of strings called zdeterminingequations.

Timing[FindDeterminingEquations[
independentvariables,dependentvariables,frozennames,p,r,xseon,
inputequation,rulesarray,internalrules]]

equationgenericvariables

invarconditiontable

ztableofrules

ztableofrulesxy

Here are the determining equations expressed in terms of z-variables. This output may be quite long. The equations in the table can be distinguished by the ==0 at the end of each item.

zdeterminingequations

How many determining equations are there?

Length[zdeterminingequations]

Lets look at the determining equations in a slightly more readable form.

zdeterminingequationscolumn=Column[ToExpression[zdeterminingequations],{Left,Below}]

Now solve the determining equations in terms of multivariable polynomials of a selected order. This function can take a long time. The Mathematica function Solve uses Gaussian elimination to solve a large number of linear equations for the polynomial coefficients. 
The time roughly follows

time/timeref=((number of equations)/(number of equationsref))^n

where the exponent is between 2.4 and 2.7. The Mathematica function Timing outputs the time required for the SolveDeterminingEquations function to execute.

Timing[SolveDeterminingEquations[
independentvariables,dependentvariables,r,xseon,
zdeterminingequations,order=3]]

Express the xse functions in terms of x,y variables.

ToExpression[xsefunctions]/.ztableofrules

Express the eta function in terms of these variables.

ToExpression[etafunctions]/.ztableofrules

Flatten[infinitesimalgroups]/.ztableofrules

Check the infinitesimals in the determining equations.

ToExpression[zdeterminingequations]

The equation is invariant under the rotation group in the plane. This is used in example 8.3 to reduce the equation to first order. The resulting first order ODE is then solved using separation of variables.

MaxMemoryUsed[ ]