In this example we use the package IntroToSymmetry.m to work out, first the point group of the Huxley equation, 

Uxx-Ut+2*U^2(1-U)=0.

then a nonclassical symmetry. 
Reference: Hydon, P. 2000 Symmetry Methods for Differential Equations Cambridge Press (see Section 9.3).

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.

huxleyequation=
"D[U[x,t],x,x]-D[U[x,t],t]+2*U[x,t]^2-2*U[x,t]^3";

The function u[x,t] is a solution of the equation and this constraint must be applied in the form of a rule to the invariance condition. Be careful to check signs.

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

Enter the list of independent variables.

independentvariables={"x","t"};

Enter the list of dependent variables.

dependentvariables={"U"};

Enter the list of function and/or constant names that must be preserved when the equation is converted to generic y1[x1,x2] variables.

frozennames={};

Enter the maximum derivative order of the input equation(s).

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.

xseon=1;

When searching for Lie-Backlund groups it is necessary to differentiate the input equation with respect to each of the independent variables 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 sets 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,huxleyequation,rulesarray,internalrules]]

equationgenericvariables

invarconditiontable

ztableofrules

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]

Now solve the determining equations in terms of multivariable polynomials. 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]]

Here are the infinitesimals for the independent variables expressed in terms of z-variables.

xsefunctions

and the infinitesimals for the dependent variables.

etafunctions

infinitesimalgroups1=infinitesimalgroups/.{z1->x,z2->t,z3->U}

List the groups.

Column[infinitesimalgroups1]

Now generate the commutator table.

MakeCommutatorTable[
independentvariables,dependentvariables,infinitesimalgroups1]

Put the commutator table in readable form.

MatrixForm[Simplify[commutatortable]]

The equation admits a two parameter group of translations in space and time. Now let's look for a nonclassical symmetry. Let G=U-f[x,t]. The function G[x,t,U] is an invariant solution of the Huxley equation if and only if it satisfies the invariance condition 

xse1(∂G/∂x)+xse2(∂G/∂t)+eta1(∂G/∂U)=0.

We will introduce this condition to the problem by adding it to rulesarray but with the ansatz xse2=1 and the infinitesimals xse1 and eta1 treated as arbitrary functions of (x,t,U) to be identified with infinitesimals later.

rulesarraynew=
{"D[U[x,t],x,x]->D[U[x,t],t]-2*U[x,t]^2+2*U[x,t]^3","D[U[x,t],t]->eta1prime[x,t,U[x,t]]-xse1prime[x,t,U[x,t]]*D[U[x,t],x]"};

frozennames={"xse1prime","eta1prime"};

Timing[FindDeterminingEquations[independentvariables,dependentvariables,frozennames,p,r,xseon,huxleyequation,rulesarraynew,internalrules]]

zdeterminingequations

Length[zdeterminingequations]

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

Now change eta1prime and xse1prime to the infinitesimals eta1 and xse1. In addition set xse2=1.

eta1prime[z1_,z2_,z3_]:=eta1[z1,z2,z3]

xse1prime[z1_,z2_,z3_]:=xse1[z1,z2,z3]

xse2[z1_,z2_,z3_]:=1

zdeterminingequationssimplified=ToExpression[zdeterminingequations];

Column[zdeterminingequationssimplified,{Left,Below}]

Note that the determining equations are nonlinear in this approach. The results above should be compared with Example 9.7 in the book "Symmetry Methods for Differential Equations", Cambridge Press by Peter Hydon. Now (knowing the answer) let's work to simplify these determining equations. Assume the infinitesimals xse1 and eta1 do not depend on x and t

zdeterminingequationssimplified1=zdeterminingequationssimplified/.
{D[xse1[z1,z2,z3],z1]->0,D[xse1[z1,z2,z3],z1,z1]->0,D[xse1[z1,z2,z3],z2]->0,D[xse1[z1,z2,z3],z1,z3]->0,D[eta1[z1,z2,z3],z2]->0,D[eta1[z1,z2,z3],z1,z1]->0}

Column[zdeterminingequationssimplified1,{Left,Below}]

Use Mathematica to solve the first of these relations.

DSolve[2*eta'[z3]*z3^3-6*eta[z3]*z3^2-2*eta'[z3]*z3^2+4*eta[z3]*z3==0,eta[z3],z3]

eta1[z1_,z2_,z3_]:=a*z3^2-a*z3^3

zdeterminingequationssimplified2=Simplify[zdeterminingequationssimplified1];

Column[zdeterminingequationssimplified2,{Left,Below}]

zdeterminingequationssimplified3=zdeterminingequationssimplified2/.{a->3}

Column[zdeterminingequationssimplified3,{Left,Below}]

DSolve[-9*z3+xse[z3]*xse'[z3]+3==0,xse[z3],z3]

The only value of the constant of integration that provides a solution that satisfies the last equation just above is C=1. Finally

xse1[z1_,z2_,z3_]=±(3*z3-1)

Check this result.

Expand[zdeterminingequationssimplified3]

This confirms the nonclassical symmetry,

xse1[z1, z2, z3_]:=±(3*z3-1)
xse2[z1, z2, z3]:=1
eta1[z1, z2, z3]:=3*z3^2-3*z3^3

with group operator

X=±(3*U-1)(∂/∂x)+(∂/∂t)+(3*U^2-3*U^3)(∂/∂U)=0.

Let's check this result in the original invariance condition.

n=2;

m=1;

r=0;

xse1[z1_,z2_,z3_]=(3*z3-1)

xse2[z1_,z2_,z3_]=1

eta1[z1_,z2_,z3_]:=3*z3^2-3*z3^3

eta1noderiv=3*y1[x1,x2]^2-3*y1[x1,x2]^3

eta1t=PthInfinitesimal[
n,m,yindex=1,plocal=1,r,xseon,j1=2]

eta1xx=PthInfinitesimal[
n,m,yindex=1,plocal=2,r,xseon,j1=1,j2=1]

invariancecondcoefficients={
0,0,eta1noderiv,0,eta1t,eta1xx,0,0}

invarianceconditionsum=
Expand[
invarconditiontable.invariancecondcoefficients]/.
{D[y1[x1,x2],x1,x1]->D[y1[x1,x2],x2]-2*y1[x1,x2]^2+2*y1[x1,x2]^3}

Simplify[%]

%/.{D[y1[x1,x2],x2]->3*y1[x1,x2]^2-3*y1[x1,x2]^3-(3*y1[x1,x2]-1)*D[y1[x1,x2],x1]}

Simplify[%]

The equation is indeed invariant under the nonclassical symmetry

xse1[z1, z2, z3_]:=±(3*z3-1)
xse2[z1, z2, z3]:=1
eta1[z1, z2, z3]:=3*z3^2-3*z3^3

with group operator

X=±(3*U-1)(∂/∂x)+(∂/∂t)+(3*U^2-3*U^3)(∂/∂U)=0.
Note that this symmetry holds only over the subset of solutions that satisfy the condition XU=0.

MaxMemoryUsed[ ]