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

Ut-(k(U))xUx-k(U)Uxx=0.

for k(U)=U^sigma. Reference: Ibragimov, N. CRC Handbook Volume 1, page 110.

Clear all symbols in the current context.

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

This next command turns off spurious spelling error warnings.

Off[General::spell]

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

Needs["SymmetryAnalysis`IntroToSymmetry`"]

Define the diffusivity function.

k[u_[x_,t_]]:=u[x,t]^sigma

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

inputequation=
"D[u[x,t],t]-D[k[u[x,t]],x]*D[u[x,t],x]-k[u[x,t]]*D[u[x,t],x,x]";

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],t]->D[k[u[x,t]],x]*D[u[x,t],x]+k[u[x,t]]*D[u[x,t],x,x]"};

Enter the list of independent variables.

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

Enter the list of dependent variables.

dependentvariables={"u"};

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

frozennames={"k","sigma"};

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

porder=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.

rorder=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.

FindDeterminingEquations[
independentvariables,dependentvariables,frozennames,porder,rorder,xseon,inputequation,rulesarray,internalrules]

equationgenericvariables

invarconditiontable

The program expresses the determining equations in terms of zvariables. Here is the correspondence between z-variables and conventional variables.

ztableofrules

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

Remove redundancies to try to simplify the determining equations.

zdeterminingequationsreduced1=ToExpression[zdeterminingequations]/.
{D[xse2[z1,z2,z3],z1]->0,
D[xse2[z1,z2,z3],z1,z1]->0,
D[xse2[z1,z2,z3],z1,z2]->0,
D[xse2[z1,z2,z3],z1,z3]->0,
D[xse2[z1,z2,z3],z3]->0,
D[xse2[z1,z2,z3],z3,z3]->0,
D[xse2[z1,z2,z3],z2,z3]->0};

TraditionalForm[Column[zdeterminingequationsreduced1,{Left,Below}]]

zdeterminingequationsreduced2=zdeterminingequationsreduced1/.
{D[xse1[z1,z2,z3],z3]->0,D[xse1[z1,z2,z3],z2,z3]->0,
D[xse1[z1,z2,z3],z3,z3]->0,D[xse1[z1,z2,z3],z1,z3]->0};

TraditionalForm[Column[zdeterminingequationsreduced2,{Left,Below}]]

The second to last equation implies that the two derivatives of eta1 are each zero.

zdeterminingequationsreduced3=zdeterminingequationsreduced2/.
{D[eta1[z1,z2,z3],z2]->0,D[eta1[z1,z2,z3],z1,z1]->0,
D[eta1[z1,z2,z3],z1,z2]->0,D[eta1[z1,z2,z3],z2,z3]->0};

TraditionalForm[Column[zdeterminingequationsreduced3,{Left,Below}]]

The last equation implies that the derivative of xse1 with respect to z2 is zero.

zdeterminingequationsreduced4=zdeterminingequationsreduced3/.
{D[xse1[z1,z2,z3],z2]->0,D[xse1[z1,z2,z3],z1,z2]->0,
D[xse1[z1,z2,z3],z2,z3]->0,D[xse1[z1,z2,z3],z2,z2]->0};

TraditionalForm[Column[zdeterminingequationsreduced4,{Left,Below}]]

Remove the factor z3^sigma that multiplies every term.

imax=Length[zdeterminingequationsreduced4]

zdeterminingequationsreduced4strings=Table[ToString[InputForm[zdeterminingequationsreduced4[[i]]]],{i,1,imax}];

zdeterminingequationsnoequalszero=
StringReplace[zdeterminingequationsreduced4strings,
{"== 0"->" "}];

Now convert each determining equation to an expression and multiply through by z3^(2-sigma). The built-in function Together cancels terms in the numerator and denominator.

zdeterminingequationsnoequalszeroexpression=
Together[(z3^2)*(z3^(-sigma))*ToExpression[
zdeterminingequationsnoequalszero]];

The number of determining equations.

imax=Length[
zdeterminingequationsnoequalszeroexpression]

Convert each determining equation back to a string and reinstate the == 0.

zdeterminingequationsmod1=Table[
ToString[InputForm[
zdeterminingequationsnoequalszeroexpression[[i]]]]<>
"== 0",
{i,1,imax}]

TraditionalForm[Column[zdeterminingequationsmod1,{Left,Below}]]

Now let's put together the final set of determining equations with redundancies removed.

zdeterminingequationsfinal={"sigma*eta1[z1, z2, z3] - sigma^2*eta1[z1, z2, z3] - sigma*z3*Derivative[0, 0, 1][eta1][z1, z2, z3] - z3^2*Derivative[0, 0, 2][eta1][z1, z2, z3] - sigma*z3*Derivative[0, 1, 0][xse2][z1, z2, z3] + 2*sigma*z3*Derivative[1, 0, 0][xse1][z1, z2, z3]== 0","-(z3*(sigma*eta1[z1, z2, z3] + z3*Derivative[0, 1, 0][xse2][z1, z2, z3] - 2*z3*Derivative[1, 0, 0][xse1][z1, z2, z3]))== 0","-(z3*(2*sigma*Derivative[1, 0, 0][eta1][z1, z2, z3] + 2*z3*Derivative[1, 0, 1][eta1][z1, z2, z3] - z3*Derivative[2, 0, 0][xse1][z1, z2, z3]))== 0",
"D[xse2[z1,z2,z3],z1]== 0",
"D[xse2[z1,z2,z3],z1,z1]== 0",
"D[xse2[z1,z2,z3],z1,z2]== 0",
"D[xse2[z1,z2,z3],z1,z3]== 0",
"D[xse2[z1,z2,z3],z3]== 0",
"D[xse2[z1,z2,z3],z3,z3]== 0",
"D[xse2[z1,z2,z3],z2,z3]== 0",
"D[xse1[z1,z2,z3],z3]== 0",
"D[xse1[z1,z2,z3],z2,z3]== 0",
"D[xse1[z1,z2,z3],z3,z3]== 0",
"D[xse1[z1,z2,z3],z1,z3]== 0",
"D[eta1[z1,z2,z3],z2]== 0",
"D[eta1[z1,z2,z3],z1,z1]== 0",
"D[eta1[z1,z2,z3],z1,z2]== 0",
"D[eta1[z1,z2,z3],z2,z3]== 0",
"D[xse1[z1,z2,z3],z2]== 0",
"D[xse1[z1,z2,z3],z1,z2]== 0",
"D[xse1[z1,z2,z3],z2,z3]== 0",
"D[xse1[z1,z2,z3],z2,z2]== 0"};

Now solve the determining equations in terms of multivariable polynomials of a selected order. 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,rorder,xseon,zdeterminingequationsfinal,order=3]]

infinitesimalgroupsxy=infinitesimalgroups/.{z1->x,z2->t,z3->u};

Column[infinitesimalgroupsxy]

Check the solution in the determining equations.

Simplify[ToExpression[zdeterminingequations]]

MakeCommutatorTable[independentvariables,dependentvariables,infinitesimalgroupsxy]

MatrixForm[commutatortable]

MaxMemoryUsed[ ]