(***********************************************************************
This file was generated automatically by the Mathematica front end.
It contains Initialization cells from a Notebook file, which typically
will have the same name as this file except ending in ".nb" instead of
".m".

This file is intended to be loaded into the Mathematica kernel using
the package loading commands Get or Needs.  Doing so is equivalent to
using the Evaluate Initialization Cells menu command in the front end.

DO NOT EDIT THIS FILE.  This entire file is regenerated automatically 
each time the parent Notebook file is saved in the Mathematica front end.
Any changes you make to this file will be overwritten.
***********************************************************************)

(* :Title:  IntroToSymmetry.m *)

(* :Author: Brian J. Cantwell (February, 1998),
Department of Aeronautics and Astronautics
Stanford University, Stanford, CA 94305. Questions,
comments and corrections should be directed to 
cantwell@stanford.edu. Version 1.0 of the package 
was originally developed in Mathematica version 2.2.2
in Jan-Feb 1998. Several improvements were made in 
September of 1998. In August 2000 additional features 
were added and the package was adapted to 
Mathematica 4.0. In May of 2001 an improved user interface 
was added to the package.This is Version 2.0 of the package.*)

(* :Copyright \[Copyright] 1998 Brian J. Cantwell.

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License version 2 June 1991 as 
    published by the Free Software Foundation.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA.*)


 (* :Summary:
This application program is used to solve for the 
determining equations of the Lie-Backlund group which 
leaves the equation 
             \[Psi]=\[CapitalPsi][x,y,y{1},y{2},...,y{p}] 
invariant. The equation is a differential function which 
depends on the independent variables,
                 x={x1,x2,....,xn},
the dependent variables,
             y={y1[x],y2[x],...,ym[x]},
first derivatives, 
            y{1}={\[PartialD]y1/\[PartialD]x1,...,\[PartialD]ym/\[PartialD]xn},
second derivatives,      
       y{2}={\[PartialD](\[PartialD]y1/\[PartialD]x1)/\[PartialD]x1,...,\[PartialD](\[PartialD]ym/\[PartialD]xn)/\[PartialD]xn},
and so on up to derivatives of order p. Note that the
equation could be one member of a system of equations.
The infinitesimal form of the transforming group is
 xj'=xj+s\[Xi]j[x,y,y{1},y{2}...,y{r}]
 yi'=yi+s\[Eta]i[x,y,y{1},y{2}...,y{r}]
 \[PartialD]yi'/\[PartialD]xj'=\[PartialD]yi/\[PartialD]xj+s\[Eta]ij[x,y,y{1}...,y{r+1}]
 \[PartialD](\[PartialD]yi'/\[PartialD]xj1')/\[PartialD]xj2'=\[PartialD](\[PartialD]yi/\[PartialD]xj1)/\[PartialD]xj2+s\[Eta]ij1j2[x,y,...]
 ......................................................
 \[PartialD](...\[PartialD](\[PartialD](pyi'/\[PartialD]xj1')/\[PartialD]xj2'...)/\[PartialD]xjp'=
 \[PartialD](...\[PartialD](\[PartialD](pyi/\[PartialD]xj1)/\[PartialD]xj2...)/\[PartialD]xjp+
             s\[Eta]ij1j2..jp[x,y,y{1},y{2},..,y{r+p}]
where s is the group parameter (assumed small) and r is 
the maximum derivative order upon which the infinitesimals 
are assumed to depend. The infinitesimals
of the point transformation (x,y)->(x',y') are \[Xi]j and \[Eta]i.
The infinitesimals of the extended transformation are
given by the formulae,
        \[Eta]ij=D\[Eta]i/Dxj-(\[PartialD]yi/\[PartialD]xk)D\[Xi]k/Dxj
     \[Eta]ij1j2=D\[Eta]ij1/Dxj-(\[PartialD](\[PartialD]yi/\[PartialD]xj1)/\[PartialD]xk)D\[Xi]k/Dxj2
\[Eta]ij1j2...jp=D\[Eta]ij1...j(p-1)/Dxjp-
        (\[PartialD]...(\[PartialD]yi/\[PartialD]xj1)/\[PartialD]xj2...\[PartialD]xj(p-1)\[PartialD]xjk)D\[Xi]k/Dxjp      
where D( )/Dxj is the operator of total differentiation 
with respect to xj. The invariance condition is
                      X\[CapitalPsi]=0
where the p-th order extended operator X is given by
  X=\[Xi]j\[PartialD]/\[PartialD]xj+\[Eta]i\[PartialD]/\[PartialD]yi+\[Eta]ij1\[PartialD]/\[PartialD]yij1+...+\[Eta]ij1j2..jp\[PartialD]/\[PartialD]yij1...jp.
The invariance condition, X\[CapitalPsi]=0, is satisfied subject to the 
condition that the y's and their derivatives that appear
in the invariance condition must satisfy the original
equation and its higher order differential consequences.
This leads to a series of transformation rules which must
be applied to the y-derivatives appearing in the invariance 
condition.

If the sought after group is of Lie-Backlund type with the
infinitesimals depending on x, y and at least y{1} then 
the \[Xi]'s can be set to zero without loss of generality. For
point groups the \[Xi]'s must be included in the computation.

Also included is a function called solvdeterminingequations
which attempts to solve the determining equations for the 
unknown infinitesimals in the form of multivariate polynomials.

  \[Xi]j[z1,z2,z3,...,zq]=
      aj1+aj2*z1+aj3*z2+aj4*z3+...+aj(q+1)*zq+
          aj(q+2)*z1^2+aj(q+3)*z1*z2+aj(q+4)*z1*z3+...

  \[Eta]i[z1,z2,z3,...,zq]=
      bi1+bi2*z1+bi3*z2+bi4*z3+...+bi(q+1)*zq+
          bi(q+2)*z1^2+bi(q+3)*z1*z2+bi(q+4)*z1*z3+...
      

where [x,y,y{1},y{2},...,y{p}]->[z1,z2,z3,...,zq] and the
aj's and bi's are constants to be solved for using the
determining equations. The variable count is 

     q=n+m*Sum[(n+k-1)!/(k!*(n-1)!), {k,0,r}]

Note that I have chosen to use 'xse' as the spelling of the
greek symbol \[Xi] rather than the usual 'xi'. The reason is
that 'xi' is too close to 'xj' which refers to the jth 
independent variable. To avoid any possible conflict, I
decided to play it safe and use 'xse'instead. The rationale 
for using a multinomial expansion is that the infinitesimals 
are the solutions of a, usually, highly overdetermined 
system of linear partial differential equations called 
the determining equations of the group. As a result, the 
infinitesimals usually turn out to be relatively simple
power multinomials. However the infinitesimals can contain 
arbitrary functions. The program can be used to identify 
such functions by searching over several orders of the trial 
multinomial. The ability to identify those groups which
are expressed as truncated multinomials or arbitrary functions
is extremely useful and a major first step toward the complete 
solution of the determining equations. But the user must be 
aware that quite often it is just a first step and there are 
many examples where the infinitesimals contain transcendental 
functions. These must be determined interactively by using 
known groups to simplify the determining equations and 
Mathematica built-in functions or adhoc methods to solve 
the remaining set.

HOW TO USE THE PACKAGE
The equation is entered as a string called 'inputequation' with the
derivatives written as Mathematica input expressions. The user is free
to pick whatever symbols are desired for the various independent and
dependent variables that appear in the equation. In addition, the 
equation or system of equations are entered by the user as a list of 
rules in string format called 'rulesarray'. The rules are constraints 
on the invariance condition and are used to replace selected derivatives 
of the dependent variable(s) that appear in the invariance condition.

The function FindDeterminingEquations first takes the input equation and 
re-expresses it in terms of generic variable names; (x1,x2,x3,...) for 
the independent variables and (y1,y2,y3,...) for the dependentvariables. 
Then the long list of terms in the invariance condition is generated. 
The rules in 'rulesarray are then applied to the various relevant derivatives of 
the dependent variables appearing in the invariance condition.The terms 
in the invariance condition are then parsed to form the determining 
equations of the group. These are available as a list of strings called
zdeterminingequations. The unknown infinitesimals are designated 'xsej' 
and 'etai'and are expressed as functions of z variables where

(z1,z2,z3,z4,....z(n+m))=(x1,x2,x3,...,xn,y1,y2,y3,...,ym).

For example the Blasius equation would be entered as:
  inputequation="D[f[x],x,x,x]+f[x]*D[f[x],x,x]"
The heat equation would be entered as
     inputequation="D[u[x,t],t]-kappa*D[u[x,t],x,x]"
In addition to the equation, there are six other inputs.

       independentvariables - the list of independent variable names, for 
       example, independentvariables={"x","y","t"}
         
       dependentvariables - the list of dependent variable names, for
       example, dependentvariables={"u","v"}
       
       frozenstrings - the list of strings that must remain unchanged
       when the input equation is changed to generic (xj,yi) variables. 
       These are usually function names and constant names that might appear 
       in the equation. For the heat equation this would be entered as
          frozennames={"kappa"}
       If no such names appear a null list is entered, frozennames={}. The
       reason this is required is that when the equation is converted to
       generic variables where, say, x is replaced by x1, a name like the
       Mathematica built-in function Exp will be mistakenly changed to Ex1p.
       The list frozennames prevents this.
         
       porder - the order of the highest derivative in the input equation
       or in the system of equations being considered.
             
       rorder - the order of the highest y-derivative that the infinitesimals 
       'xse' and 'eta' can depend on.
         
       xseon - This parameter is zero if the xse's are set to zero, and 
       one if the xse's are not zero. If the search is for Lie point groups 
       then 'rorder=0' and 'xseon' must be set to one. In the case of 
       Lie-Backlund groups with 'rorder=1'or higher, then one can, without 
       loss of generality, transform only the independent variables.In this 
       case 'xseon' can be set to zero or one. If zero is used then only 
       the dependent variables are transformed.
         
       rulesarray - This is an array of replacement rules in string format 
       used to apply the equation to the invariance condition as a constraint.
       If this is not done then the invariance condition is too restrictive 
       and symmetries may be lost.
       
       internalrules - This parameter is set to zero if no further rules 
       need to be constructed. This is the case if one is searching for point
       groups. But if the package is being used to search for Lie-Backlund
       groups with rorder=1 or more then internalrules needs to be set to 
       one so that the differential consequences of the equation in question are generated 
       internally and appended to 'rulesarray'. The enlarged set of rules is
       contained in a table called 'rulesarrayexpanded'.
         
Once the data has been entered, the user calls the function
   FindDeterminingEquations[
       independentvariables,dependentvarables,frozenstrings,porder,rorder,
       xseon,inputequation,rulesarray,internalrules].

Several tables are generated as output. The first is called

                 invarconditiontable
                 
and contains the various derivatives of the input equation.
The entries in this table are expressed in terms of generic
x1,z2,x3,... independent variables and generic y1,y2,y3,...
dependent variables.The second table is called

            termsoftheInvarianceConditionNoRules
            
which, as the name implies, contains all the terms in the 
invariance condition before the application of transformation 
rules. These rules are contained in the table 'rulesarray' 
input by the user. The program can also generate a table of 
rules corresponding to the higher order differential
consequences of the original equation, ie, derivatives of the 
equation with respect to all possible combinations of 
independent variables up to the order 'rorder'selected by the user. 
The expanded set of rules, which includes the original set
input by the user, is generated by the routine 
'MakeRulesArray'. The expanded set of rules is contained 
in the table 'rulesarrayexpanded'. If the parameter 'internalrules=1' 
the internally generated rules are used. If 'internalrules=0' then 
only the rules in 'rulesarray' are used. Thus the user has 
maximum flexibility in choosing the rules used to restrict 
the invariance condition. The invariance condition 
after the application of the rules is available in the table

         termsoftheInvarianceConditionRulesApplied
         
The y-derivative terms are separated from the infinitesimal
derivative terms to produce two final tables. The y-derivative 
products are contained in the table

                   yderivfactortable
                   
and the determining equations are contained in the table.

                  determiningequations
                  
For convenience, the variables upon which the infinitesimals
depend are converted to a set of dummy variables, ie,
x1->z1, x2->z2, y1->z3, etc. The result is contained in

                  zdeterminingequations
                  
and the correspondence between variables is contained in

                     ztableofrules
                     
These tables are the main output of the package.

OTHER USES
The package can also be used simply to generate a 
table of variables including derivatives up to order porder. 
This is done simply by calling the function

              GenerateVariableTable[n_,m_,porder_].
              
where n is the number of independent variables and m is the
number of dependentvariables. Similarly, one can generate 
the terms of any infinitesimal by calling the function

     PthInfinitesimal[n_,m_,yindex_,p_,r_,xseon_,
     j1_,j2_:0,j3_:0,j4_:0,j5_:0,j6_:0,j7_:0,j8_:0,
     j9_:0,j10_:0,j11_:0,j12_:0,j13_:0,j14_:0]
     
The input yindex is the index of the particular
dependent variable in question. The indices 'j1' through
'j14' are the lower indices of the infinitesimal 
corresponding to the derivative being transformed. These
indices are optional and any missing ones are simply set to 
zero. The limit of 14 seems to be a limitation of the 
application Mathematica which only allows up to 13 optional arguements 
in a defined function. The fact is, almost any computer will 
be brought to its knees if the user tries to compute Lie-Backlund 
transformations that depend on derivatives above 4th or 5th order.

SOLVING THE DETERMINING EQUATIONS
A multivariable polynomial solution for the infinitesimals can 
be attempted by calling the function

     SolveDeterminingEquations[
          independentvariables,dependentvariables,rorder,xseon,
          zdeterminingequations,polyorder]
          
The parameter 'polyorder' is used to select the order of the 
polynomial used in the expansion. One approach is to call the 
function succesively with increasing orders, stopping when the 
computed infinitesimals stop changing. The results are presented 
as two tables of strings, xsefunctions and etafunctions ordered 
in the usual way. These strings can be converted to expressions 
using the built-in function 'ToExpression'. 

SAMPLE RUN
A typical session for calculating the point group of the heat 
equation would be as follows.
Inputs
independentvariables={"x","t"}
dependentvariables={"u"}
frozenstrings={"kappa"}
porder=2
rorder=0
xseon=1
internalrules=0
inputequation="D[u[x,t],t]-kappa*D[u[x,t],x,x]"
rulesarray={"D[u[x,t],x,x]->(1/kappa)*D[u[x,t],t]"}

Call the function
FindDeterminingEquations[
independentvariables,dependentvariables,frozenstrings,porder,rorder,xseon,
inputequation,rulesarray,internalrules]

Main Output
determiningequations
zdeterminingequations
ztableofrules

Try to solve for xse1, xse2 and eta1
polyorder=3
SolveDeterminingEquations[
independentvariables,dependentvariables,rorder,
xseon,zdeterminingequations,polyorder]

Output
xsefunctions
etafunctions
infinitesimalgroups

Express the table infinitesimal groups in terms of user variables
infinitesimalgroups1=infinitesimalgroups/.{z1->x,z2->y,z3->u}

Construct the commutator table of the group
MakeCommutatorTable[
independentvariables,dependentvariables,infinitesimalgroups1]

Output
commutatortable
*)

(*Prepend the context SymmetryAnalysis` to $ContextPath*)

BeginPackage["SymmetryAnalysis`IntroToSymmetry`"]

GenerateVariableTable::usage=
"GenerateVariableTable[numberindependentvars,numberdependentvars,porderoftheequation] This function generates a 
table of all the various variables in x,y,y{1},y{2},...,y{porderoftheequation}.
The output is a table of strings contained in 
variablestringtable.";

SecondTerm::usage=
"SecondTerm[numberindependentvars,numberdependentvars,yindex,plocal,rorderofinfinitesimals,
j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,j13,j14] This function 
generates the second term in the p-th order infinitesimal.
Note that derivatives up to 14th order are permitted by
program.";

PthInfinitesimal::usage=
"PthInfinitesimal[numberindependentvars,numberdependentvars,yindex,porderoftheequation,rorderofinfinitesimals,xseon,
j1,j2,j3,j4,j5,j6,j7,j8,j9,j10,j11,j12,j13,j14] This nested
function generates the first term in the p-th order 
infinitesimal. In the process it calls the function
'SecondTerm'.";

GenerateInfinitesimalTable::usage=
"GenerateInfinitesimalTable[numberindependentvars,numberdependentvars,porderoftheequation,rorderofinfinitesimals,xseon] 
This function generates a table of strings with the same 
indices as 'GenerateVariableTable'. Each element in  the 
table is a call to 'PthInfinitesimal' which will be 
evaluated once it is determined which infinitesimals will be 
needed.";

InvarianceConditionNoRules::usage=
"InvarianceConditionNoRules[numberindependentvars,numberdependentvars,porderoftheequation,rorderofinfinitesimals,xseon,inputequation] 
This function generates a table of expressions for the terms 
in the invariance condition before the transformation rules 
have been applied. The result is contained in 
'termsoftheInvarianceConditionNoRules'.";

InvarianceConditionRulesApplied::usage=
"InvarianceConditionRulesApplied[
numberindependentvars,numberdependentvars,porderoftheequation,rorderofinfinitesimals,xseon,inputequation,rulesarray,internalrules]
This function generates a table of expressions for the terms in 
the invariance condition after the transformation rules have 
been applied. The result is contained in 
'termsoftheInvarianceConditionRulesApplied1'
and 'termsoftheInvarianceConditionRulesApplied'.In the first table
each factor is multiplied by two to insure that every term has at
least two factors to insure the proper operation of the function 
'Cases' which is used in the next step. In the second table the
factor of two has been removed.";

MakeRulesArray::usage=
"MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98,rorderofinfinitesimals,rulesarray] This function
generates the array of rules which are applied to the invariance
condition in the case of a Lie-Backlund transformation. The rules
are standard in the sense they are simply all the possible
derivatives of the input equation up to the order of the
Lie-Backlund variables. The results are contained in the
table rulesarrayout.";

FindDeterminingEquations::usage=
"FindDeterminingEquations[
independentvariables,dependentvariables,frozenstrings,porderoftheequation,rorderofinfinitesimals,xseon,
inputequationuservariables,rulesarrayuservariables,internalrules]
This is the main function of the program. The results are 
presented as strings in the table 'determiningequations' 
which contains the determining equations of the group. The table 
'yderivfactortable' contains expressions for the corresponding 
products of free variables. For convenience, the variables 
upon which the infinitesimals depend are converted to a set of 
dummy variables,ie, x1->z1, x2->z2, y1->z3, etc. The result is 
a table of strings contained in 'zdeterminingequations'. 
The correspondence between variables is contained in 
'ztableofrules'. These tables are the main output of the 
package.";

SolveDeterminingEquations::usage=
"SolveDeterminingEquations[
independentvariables,dependentvariables,rorderofinfinitesimals,xseon,
zdeterminingequationstable,order] This
function can be used to try to solve the determining equations
contained in the array 'zdeterminingequations'. The solution
is sought in the form of a polynomial expansion in powers of the 
z-variables. This function should be thought of as a taking a 
first cut at the general solution of the determining equations 
which can admit functions that are beyond the reach of a truncated 
power series expansion. The results are contained in the tables
'xsefunctions' and 'etafunctions'. The expressions in 
these tables are ordered according to the usual ordering of
independent and dependent variables.";

MakeCommutatorTable::usage=
"MakeCommutatorTable[
independentvariables,dependentvariables,infinitesimalgroupsxy] This
function is used to generate the commutator table of the Lie 
algebra of the groups contained in 'infinitesimalgroupsxy'.
The symbols used for 'independentvariables' and
'dependentvariables' must correspond to the symbols in
'infinitesimalgroupsxy'. The output is contained in the table
'commutatortable'. This is the last function in the package.";



(*Many of the names used are similar to one another.
The function 'Off' prevents the application from 
generating a bunch of annoying warnings.*)
Off[General::spell1];

(******BEGINNING*OF*FUNCTION*GenerateVariableTable*******)
(*This is the name of the function that generates the
table of variables beginning with independent variables,
followed by dependent variables, then by derivatives up
to order p. The result is generated in the table
"variabletable".*)
GenerateVariableTable[numberindependentvars_,numberdependentvars_,porderoftheequation_]:=

(*The function "Module" hides the variable names in 
brackets from the rest of the program*) 
(Module[{indepvariablestring,indepvariablestring1,
dependvariablestring,makederivtablestring,derivindtable,
derivindsumtable,postbl,postbllength,
yderivtableofstrings1,
yderivtableofstrings2,yderivtableofstrings3,
yderivtableofstrings4,yderivtableofstrings,i,j,jindex,
ppndx,indtbl,indtbl2,indtbl3,indexm},

(*Generate a table of independent variables.*)
indepvariablestring=
Table["x"<>ToString[j],{j,1,numberindependentvars}];

(*Generate an unbracketed string of independent variables.*)
indepvariablestring1=
StringReplace[
StringReplace[
StringReplace[
ToString[Table["x"<>ToString[j],{j,1,numberindependentvars}]],
"{"->""],"}"->""]," "->""];

(*Generate a table of dependent variables.*)
dependvariablestring=
Table["y"<>ToString[i],{i,1,numberdependentvars}];

(*Generate a table of dependent variables each with the 
bracketed string of independent variables attached to 
indicate the dependence on the xj's.*)
dependvariablestring1=
Table["y"<>ToString[i]<>"["<>
indepvariablestring1<>"]",{i,1,numberdependentvars}];

(*This step prepares a string which in the next line is
converted to a function.*)
makederivtablestring="Table["<>
ToString[Table["j"<>ToString[
jindex],{jindex,1,numberindependentvars}]]<>", "<>
StringReplace[
StringReplace[
ToString[
Table[{"j"<>ToString[numberindependentvars+1-jindex]<>",0,"<>ToString[porderoftheequation]},
{jindex,1,numberindependentvars}]],"{{"->"{"],"}}"->"}"]<>"]";

(*This step generates a table of lists of derivative
indices.*)
derivindtable=
Flatten[
ToExpression[makederivtablestring],numberindependentvars-1];

(*It is necessary to identify those combinations of
indices which add up to the derivative order, p.
First add together the indices in each term. The
list is flattened to remove extra brackets and
create a one-dimensional list.*)
derivindsumtable=
Flatten[
ToExpression[
"Table["<>
ToString[
Apply[Plus,Table["j"<>ToString[jindex],{jindex,1,numberindependentvars}]]]<>
", "<>
StringReplace[
StringReplace[
ToString[
Table[{"j"<>ToString[jindex]<>",0,"<>ToString[porderoftheequation]},
{jindex,1,numberindependentvars}]],"{{"->"{"],"}}"->"}"]<>"]"],numberindependentvars-1];

(*This step identifies the positions of those terms
where the indices add up to 1,2,3...,p. These are
the terms corresponding to the various order
derivatives. This list is also flattened to make the
list one-dimensional.*)
postbl=Flatten[
Table[
Flatten[
Position[derivindsumtable,ppndx],1],{ppndx,1,porderoftheequation}],1];

(*The number of terms with the right combinations
of derivatives.*)
postbllength=Length[postbl];

(*This picks off the needed terms from the rather long
table produced above. The table 'derivativeindextable' 
is available as an output of this function, ie, it is
not protected by 'Module'.*)
derivativeindextable=
Table[derivindtable[[postbl[[i]]]],
{i,1,postbllength}];

(*At this point the name structure of the variables
is created by inserting yi before each list.*)
yderivtableofstrings1=Flatten[Table[
Table[
ToString[derivativeindextable[[indtbl]]]<>
"[y"<>ToString[indexm]<>"]",
{indtbl,1,postbllength}],{indexm,1,numberdependentvars}],1];

(*Now exchange the braces surrounding the lists of indices
for brackets and remove blanks between numbers.*)
yderivtableofstrings2=Table[
StringReplace[
StringReplace[
StringReplace[
yderivtableofstrings1[[indtbl2]],
"{"->"["],
"}["->"]["],
" "->""],
{indtbl2,1,Length[yderivtableofstrings1]}];

(*Each variable needs to be put into the form of an 
input expression which the application can recognize.
To accomplish this, add the word 'Derivative' in
front of the brackets and the string of bracketed
independent variables after the name of the 
independent variable.*)
yderivtableofstrings3=Table[
StringReplace[
yderivtableofstrings2[[indtbl3]],
"]"->"]["<>indepvariablestring1<>"]"],
{indtbl3,1,Length[yderivtableofstrings2]}];

(*Get rid of the extra list of independent variables.*)
yderivtableofstrings4=Table[
StringReplace[
yderivtableofstrings3[[indtbl3]],
"]["<>indepvariablestring1<>"]["->"]["],
{indtbl3,1,Length[yderivtableofstrings2]}];

(*Now add the word 'Derivative' at the front of each 
element of the table.*)
yderivtableofstrings=Table[
"Derivative"<>yderivtableofstrings4[[indtbl3]],
{indtbl3,1,Length[yderivtableofstrings2]}];

(*Join the tables of independent variables, dependent
variables and derivatives together.*)
variablestringtable=Join[
indepvariablestring,
dependvariablestring1,
yderivtableofstrings];

(*Produce a version of the table without the independent
variables.*)
withoutindvarstringtable=Join[
dependvariablestring1,
yderivtableofstrings];])
(**********END*OF*FUNCTION*GenerateVariableTable*******)

(************BEGINNING*OF*FUNCTION*SecondTerm**********)
(*This function creates the plocal-th order second 
term in the extended infinitesimal*)
SecondTerm[numberindependentvars_,numberdependentvars_,yindex_,plocal_,rorderofinfinitesimals_,j1_,j2_:0,j3_:0,j4_:0,j5_:0,j6_:0,j7_:0,j8_:0,j9_:0,j10_:0,j11_:0,j12_:0,j13_:0,j14_:0]:=

(*The function 'Module' hides the variable names in 
brackets from the rest of the program.*)
(Module[{indepvariables,dependvariables,derivaar,
beta,jindex,ind,j,i},


(*This generates a string corresponding to the 
independent variables.*)
indepvariables=StringReplace[
StringReplace[
ToString[
Table["x"<>ToString[j],{j,1,numberindependentvars}]],"{"->""],"}"->""];

(*Now use 'GenerateVariableTable' to create a 
table of dependent variables with brackets containing
the independent variables to indicate the 
functional dependence of the dependent variables.*)
GenerateVariableTable[numberindependentvars,numberdependentvars,rorderofinfinitesimals];

(*This is the table.*)
withoutindvarstringtable;

(*Now convert the table to a string.*)
dependvariables=" "<>
StringReplace[
StringReplace[
ToString[withoutindvarstringtable],
"{"->""],"}"->""];

(*This generates a table of indices corresponding
to the independent variables with respect to which
the p-th extension is generated.*)
derivaar=Table[
ToExpression["j"<>ToString[ind]],{ind,1,plocal}];

(*Now generate the sum of y's times xse's in the
second term of the plocal-th infinitesimal.*)
Apply[
Plus,
Table[
D[ToExpression["xse"<>ToString[beta]<>
"["<>indepvariables<>","<>dependvariables<>"]"],
ToExpression["x"<>ToString[derivaar[[plocal]]]]]*
D[ToExpression[
"D[y"<>ToString[yindex]<>"["<>indepvariables<>"]"<>
StringJoin[
Table[
", x"<>ToString[derivaar[[jindex]]],
{jindex,1,plocal-1}]]<>"]"],
ToExpression["x"<>ToString[beta]]],
{beta,1,numberindependentvars}]]])
(*****************END*OF*FUNCTION*SecondTerm*********)

(*******BEGINNING*OF*FUNCTION*PthInfinitesimal*******)
(*This function creates the p-th order extended 
infinitesimal corresponding to the transformation
of some derivative of the dependent variable 'yindex'.*)
PthInfinitesimal[numberindependentvars_,numberdependentvars_,yindex_,porderoftheequation_,rorderofinfinitesimals_,xseon_,j1_,j2_:0,j3_:0,j4_:0,j5_:0,j6_:0,j7_:0,j8_:0,j9_:0,j10_:0,j11_:0,j12_:0,j13_:0,j14_:0]:=

(*The function 'Module' hides the variable names in 
brackets from the rest of the program.*)
(Module[{indepvariables,dependvariables,
firstpthorderterm,pthextendedinfinitesimal,
j,i,pind,plocalind,indepvarindex},

(*This generates a string corresponding to the 
independent variables.*)
indepvariables=StringReplace[
StringReplace[
ToString[
Table["x"<>ToString[j],{j,1,numberindependentvars}]],"{"->""],"}"->""];

(*Now use 'GenerateVariableTable' to create a 
table of dependent variables with brackets containing
the independent variables to indicate the 
functional dependence of the dependent variables.*)
GenerateVariableTable[numberindependentvars,numberdependentvars,rorderofinfinitesimals];

(*This is the table.*)
withoutindvarstringtable;

(*Now convert the table to a string.*)
dependvariables=" "<>
StringReplace[
StringReplace[
ToString[withoutindvarstringtable],
"{"->""],"}"->""];

(*This nested function is generated by concatenating
a series of strings. The function generates the first 
term in the p-th order infinitesimal. The 'If' statement,
depending on the value of the parameter 'xseon', determines 
whether the infinitesimal transformations of the 
independent variables (the xse's) are to be included in 
the computation.*)
If[xseon>0,
firstpthorderterm=
ToExpression[Table["D[",{pind,1,porderoftheequation}]<>
"eta"<>ToString[yindex]<>
"["<>indepvariables<>","<>
dependvariables<>"]"<>", x"<>
ToString[j1]<>"]"<>
StringJoin[
Table[
"-SecondTerm["<>ToString[numberindependentvars]<>","<>ToString[numberdependentvars]<>","<>
ToString[yindex]<>","<>ToString[plocalind]<>","<>
ToString[rorderofinfinitesimals]<>","<>
StringReplace[StringReplace[
ToString[Table["j"<>ToString[indepvarindex]<>"="<>
ToString[ToExpression["j"<>ToString[indepvarindex]]],
{indepvarindex,1,plocalind}]],
"{"->" "],"}"->""]<>"], x"<>
ToString[ToExpression["j"<>ToString[plocalind+1]]]<>"]",
{plocalind,1,porderoftheequation-1}]]],
firstpthorderterm=
ToExpression[Table["D[",{pind,1,porderoftheequation}]<>
"eta"<>ToString[yindex]<>
"["<>indepvariables<>","<>
dependvariables<>"]"<>", x"<>
ToString[j1]<>"]"<>
StringJoin[
Table[", x"<>ToString[ToExpression[
"j"<>ToString[plocalind+1]]]<>"]",
{plocalind,1,porderoftheequation-1}]]]];

(*The expression in the previous step is now added to
the second expression in the p-th order infinitesimal.
If xseon is zero then there is no second term.*)
If[xseon>0,
pthextendedinfinitesimal=
ExpandAll[firstpthorderterm-
SecondTerm[numberindependentvars,numberdependentvars,yindex,porderoftheequation,rorderofinfinitesimals,
j1,j2,j3,j4,j5,j6,j7,j8,
j9,j10,j11,j12,j13,j14]],
pthextendedinfinitesimal=
ExpandAll[firstpthorderterm]]])
(**********END*OF*FUNCTION*PthInfinitesimal************)

(**BEGINNING*OF*FUNCTION*GenerateInfinitesimalTable****)
(*This function generates a table of strings 
corresponding to unevaluated calls to the function
'PthInfinitesimal'. The indices of derivatives are the
same as the indices produced by 'GenerateVariableTable'.
The result is contained in the table 
'infinitesimaltableofstringsally'.*)
GenerateInfinitesimalTable[numberindependentvars_,numberdependentvars_,porderoftheequation_,rorderofinfinitesimals_,xseon_]:=

(*The function "Module" hides the variable names in 
brackets from the rest of the program.*) 
(Module[{indepvariablestring,indepvariablestring1,
dependvariablestring,dependvariablestring1,
variabletableofstrings,stringofvariables,
tableofxses,tableofetas,makederivtablestring,
derivindtable,derivindsumtable,postbl,postbllength,
derivativeindextable,pindextable1,pindextable,
pindextableoflengths,pindextableofstrings,
infinitesimaltableofstrings,
infinitesimaltableofstringsally,
i,j,k,jindex,ppndx,pconstruct},

(*Generate a table of independent variables.*)
indepvariablestring=
Table["x"<>ToString[j],{j,1,numberindependentvars}];

(*Generate a string of independent variables.*)
indepvariablestring1=
StringReplace[
StringReplace[
ToString[Table["x"<>ToString[j],{j,1,numberindependentvars}]],
"{"->""],"}"->""];

(*Now use 'GenerateVariableTable' to create a 
table of dependent variables with brackets containing
the independent variables to indicate the 
functional dependence of the dependent variables.*)
GenerateVariableTable[numberindependentvars,numberdependentvars,rorderofinfinitesimals];

(*This is the table.*)
withoutindvarstringtable;

(*Equate this table to the table with the name 
used in the next command.*)
dependvariablestring1=withoutindvarstringtable;

(*Generate a string with independent variables
and dependent variables with brackets enclosing 
the independent variables.*)
variabletableofstrings=Join[indepvariablestring,
dependvariablestring1];

(*Generate a string of bracketed independent and
dependent variables.*)
stringofvariables=StringReplace[StringReplace[
ToString[variabletableofstrings],
"{"->"["],"}"->"]"];

(*Generate a table of xse's. If 'xseon' is zero 
then the xse's are all set to zero.*)
If[xseon>0,
tableofxses=Table["xse"<>ToString[j]<>stringofvariables,
{j,1,numberindependentvars}],
tableofxses=Table["0",{j,1,numberindependentvars}]];

(*Generate a table of eta's.*)
tableofetas=Table["eta"<>ToString[i]<>stringofvariables,
{i,1,numberdependentvars}];

(*This step prepares a string which in the next line is
converted to a function.*)
makederivtablestring="Table["<>
ToString[Table["j"<>ToString[
jindex],{jindex,1,numberindependentvars}]]<>", "<>
StringReplace[
StringReplace[
ToString[
Table[{"j"<>ToString[numberindependentvars+1-jindex]<>",0,"<>ToString[porderoftheequation]},
{jindex,1,numberindependentvars}]],"{{"->"{"],"}}"->"}"]<>"]";

(*This step generates a table of lists of derivative
indices*)
derivindtable=
Flatten[
ToExpression[makederivtablestring],numberindependentvars-1];

(*It is necessary to identify those combinations of
indices which add up to the derivative order, p.
First add together the indices in each term. The
list is flattened to remove extra brackets and
create a one-dimensional list.*)
derivindsumtable=
Flatten[
ToExpression[
"Table["<>
ToString[
Apply[Plus,Table["j"<>ToString[jindex],{jindex,1,numberindependentvars}]]]<>
", "<>
StringReplace[
StringReplace[
ToString[
Table[{"j"<>ToString[jindex]<>",0,"<>ToString[porderoftheequation]},
{jindex,1,numberindependentvars}]],"{{"->"{"],"}}"->"}"]<>"]"],numberindependentvars-1];

(*This step identifies the positions of those terms
where the indices add up to 1,2,3...,p. These are
the terms corresponding to the various order
derivatives. This list is also flattened to make the
list one-dimensional.*)
postbl=Flatten[
Table[
Flatten[
Position[derivindsumtable,ppndx],1],{ppndx,1,porderoftheequation}],1];

(*The number of terms with the right combinations
of derivatives.*)
postbllength=Length[postbl];

(*This picks off the needed terms from the rather long
table produced above.*)
derivativeindextable=
Table[derivindtable[[postbl[[i]]]],
{i,1,postbllength}];

(*Each element of the table just generated is itself a
table with 'n' positions with the numbers in each element
adding up to 'porderoftheequation'. Now we need to rearrange each element
into a table with 'porderoftheequation' positions corresponding to the
indices of the independent variables with respect to 
which the derivative is being taken.*)
pindextable1=Table[Table[Table[j,
{i,1,derivativeindextable[[k,j]]}],{j,1,numberindependentvars}],
{k,1,postbllength}];

(*The next step gets rid of empty brackets and 
consolidates terms of the same order in 'porderoftheequation'. This
is the sought after set of derivative indices.*)
pindextable=Table[Flatten[pindextable1[[k]],1],
{k,1,postbllength}];

(*Make a table with all the lengths at each order 'porderoftheequation'.*)
pindextableoflengths=Table[
Length[pindextable[[k]]],
{k,1,postbllength}];

(*Make a table of strings containing all the derivative
index assignments.*)
pindextableofstrings=Table[","<>"j"<>ToString[
jindex]<>
"="<>ToString[pindextable[[k,jindex]]],
{k,1,postbllength},
{jindex,1,pindextableoflengths[[k]]}];

(*Now make up a table of calls to the function
'PthInfinitesimal' with the various parameters set
except for yindex. The use of the variable pconstruct prevents
the parameter porderoftheequation from being changed when the calls to PthInfinitesimal
are executed later.*)
infinitesimaltableofstrings=Table[
"PthInfinitesimal"<>
"["<>"numberindependentvars="<>ToString[numberindependentvars]<>","<>"numberdependentvars="<>ToString[numberdependentvars]<>","<>
"yindex="<>","<>"pconstruct="<>ToString[
pindextableoflengths[[k]]]<>","<>
"rorderofinfinitesimals="<>ToString[rorderofinfinitesimals]<>","<>
"xseon="<>ToString[xseon]<>
pindextableofstrings[[k]]<>"]",
{k,1,postbllength}];

(*Now repeat the table for various 'yi'.*)
infinitesimaltableofstringsally=Flatten[
Table[
StringReplace[
infinitesimaltableofstrings[[k]],
"yindex="->"yindex="<>ToString[i]],
{i,1,numberdependentvars},{k,1,postbllength}],1];

(*Join the tables of xse's, eta's and 
porderoftheequation-th order infinitesimals together.*)
allpinfinitesimaltableofstrings=Join[
tableofxses,
tableofetas,
infinitesimaltableofstringsally];])
(****END*OF*FUNCTION*GenerateInfinitesimalTable****)

(**BEGINNING*OF*FUNCTION*InvarianceConditionNoRules**)
(*This function generates a table of expressions
corresponding to the various terms in the raw 
invariance condition, ie, before rules are applied which
impose the original equation and its higher order 
consequences on the y-derivatives which appear in 
the invariance condition.*)
InvarianceConditionNoRules[
numberindependentvars_,numberdependentvars_,porderoftheequation_,rorderofinfinitesimals_,xseon_,inputequation_]:=

(*The function "Module" hides the variable names in 
brackets from the rest of the program.*)
(Module[{inputequationexpression,allvariablestable,
allvariablestablelength,replacementvariabletable,
tableofrulesstrings,tableofrules,
replacementinputequation,replacementderivativestable,
reversetableofrulesstrings,reversetableofrules,idx,
invarconditiontablelength,masktable,
finalinfinitesimaltableofstrings1,
finalinfinitesimaltableofstrings,
finalinfinitesimaltable,
finalinvarianceconditiontable,numberofterms,
inx,ifinally},

(*Convert the string 'inputequation' put in by the user
to an expression.*)
inputequationexpression=ToExpression[inputequation];

(*Call the function 'GenerateVariableTable' to produce
the output table 'variablestringtable' which contains all 
of the possible variables up to derivatives of order 'porderoftheequation' 
which could appear in the input equation.*)
GenerateVariableTable[numberindependentvars,numberdependentvars,porderoftheequation];

(*Convert the terms in this table to expressions.*)
allvariablestable=ToExpression[variablestringtable];

(*Find the number of terms in the table.*)
allvariablestablelength=Length[allvariablestable];

(*In order to produce the invariance condition it is
necessary to differentiate the input equation with
respect to the various variables in the problem
including derivatives. But a problem comes up when 
differentiating the equation with respect to the
independent variables which might appear explicitely
in the equation. The problem is that the operator 'D'
is a total differentiation operator and therefore 
differentiates with respect to 'x' wherever it might 
appear including in the expressions indicating the 
dependence of the 'y' variables. To prevent this a 
correspondence is made between the variables in the 
problem and a set of replacement variables which 
can readily be treated as independent. Here a table of 
replacement variables of the required length is
generated.*)
replacementvariabletable=Table[
ToExpression["replacementvariable"<>ToString[idx]],
{idx,1,allvariablestablelength}];

(*Using string functions, a table of replacement rules 
for the variables is created.*)
tableofrulesstrings=Table[
ToString[InputForm[allvariablestable[[idx]]]]<>
"->"<>ToString[replacementvariabletable[[idx]]],
{idx,1,allvariablestablelength}];

(*Convert the replacement rules to expressions.*)
tableofrules=ToExpression[tableofrulesstrings];

(*Use the rules just generated to express the input
equation in terms of the replacement variables.*)
replacementinputequation=inputequationexpression/.tableofrules;

(*Create a table of derivatives of the input equation.*)
replacementderivativestable=Table[D[replacementinputequation,
ToExpression["replacementvariable"<>ToString[idx]]],
{idx,1,allvariablestablelength}];

(*Now generate a table of rules that reverse the 
assignment of variables so that the table of 
derivatives can be expressed in terms of conventional
variables.*)
reversetableofrulesstrings=Table[
ToString[replacementvariabletable[[idx]]]<>
"->"<>ToString[InputForm[allvariablestable[[idx]]]],
{idx,1,allvariablestablelength}];

(*Convert these reverse replacement rules to expressions.*)
reversetableofrules=
ToExpression[reversetableofrulesstrings];

(*Now express the table of derivatives of the input 
equation in terms of conventional variables.*)
invarconditiontable=
replacementderivativestable/.reversetableofrules;

(*Determine the length of the 'invarconditiontable'.*)
invarconditiontablelength=Length[invarconditiontable];

(*Produce a table of strings corresponding to the
infinitesimals and extensions. This has the same 
length as the 'invarconditiontable'.*)
GenerateInfinitesimalTable[numberindependentvars,numberdependentvars,porderoftheequation,rorderofinfinitesimals,xseon];

(*This is the table of infinitesimals and extensions.*)
allpinfinitesimaltableofstrings;

(*So far the infinitesimals and extensions are in 
the form of a table of strings which need to be converted
to evaluated functions. But we dont want to evaluate all
the expressions in the table, only those which correspond
to nonzero terms in the invariance condition. To begin
the process we create a mask in the form of a table of
zeros and one's with the one's corresponding to the nonzero
terms in the invariance condition. Note the use of
the triple equals sign to test equality symbolically.*)
masktable=Table[If[invarconditiontable[[inx]]===0,
0,1],{inx,1,invarconditiontablelength}];

(*Now produce a reduced list of infinitesimals 
corresponding to the nonzero terms in the invariance 
condition.*)
finalinfinitesimaltableofstrings1=masktable*
allpinfinitesimaltableofstrings;

(*Convert all terms including the zeros to strings.*)
finalinfinitesimaltableofstrings=Table[
ToString[finalinfinitesimaltableofstrings1[[inx]]],
{inx,1,invarconditiontablelength}];

(*Now evaluate the terms in this table.*)
finalinfinitesimaltable=Table[ToExpression[
finalinfinitesimaltableofstrings[[inx]]],
{inx,1,invarconditiontablelength}];

(*Here at last is the fully worked out invariance condition
in the form of a table where each term is a fully worked 
out term in the invariance condition.*)
finalinvarianceconditiontable=finalinfinitesimaltable*
invarconditiontable;

(*Now create a single sum of all these terms.*)
finalinvarianceconditionsumnorules=ExpandAll[Apply[Plus,
finalinvarianceconditiontable]];

(*Determine the number of terms in the fully expanded
invariance condition sum.*)
numberofterms=Length[finalinvarianceconditionsumnorules];

(*Finally the sum is sorted into a table of individual
terms.*)
termsoftheInvarianceConditionNoRules=Table[
finalinvarianceconditionsumnorules[[ifinally]],
{ifinally,1,numberofterms}];])
(******END*OF*FUNCTION*InvarianceConditionNoRules*****)

(*********BEGINNING*OF*FUNCTION*MakeRulesArray********)
(*This function takes one of the equations in 'rulesarray'
input by the user and generates all the various
higher order derivatives of the input equation up to
order 'rorder'. These derivatives are used to 
form an expanded set of rules which are appended to the
set put in by the user. The output of this function
is a table of rules called 'rulesarrayexpanded'.*)
MakeRulesArray[
numberindependentvars_,numberdependentvars_,mequindexsep98_,rorderofinfinitesimals_,rulesarray_]:=

(*The function 'Module' protects certain names from
the outside world.*)
(Module[{rulesarraystring,rulesarraystringlength,
arrowposition,righthandterm,lefthandterm,
derivativeindicessofar,derivativeindsstrings,
lefthandtermderivativetable,
righthandtermderivativetable1,
righthandtermderivativetable,
rulesarrayexpandedtableofstrings,
rulesarrayexpandedtable1,
mequationindexsep98},

(*Generate higher order differential consequences
for the mth equation*)
mequationindexsep98=mequindexsep98;

(*Turn the mequindexsep98 rule in 'rulesarray' into a string.*)
rulesarraystring=ToString[
InputForm[rulesarray[[mequationindexsep98]]]];

(*Find the number of characters in this string.*)
rulesarraystringlength=StringLength[rulesarraystring];

(*Locate the positions of the two characters which make up 
the arrow in the rule.*)
arrowposition=Flatten[StringPosition[
rulesarraystring,"->"]];

(*Pick off the right hand term and convert it to an 
expression.*)
righthandterm=ToExpression[StringTake[
rulesarraystring,
{arrowposition[[2]]+1,rulesarraystringlength}]];

(*Pick off the left hand term and convert it to an
expression.*)
lefthandterm=ToExpression[StringTake[
rulesarraystring,{1,arrowposition[[1]]-1}]];

(*Call this function to generate a table of
derivative indices up to order 'rorderofinfinitesimals'.*)
GenerateVariableTable[numberindependentvars,numberdependentvars,rorderofinfinitesimals];

(*This is the table of derivative indices.*)
derivativeindextable;

(*This is the number of derivatives of the input 
equation required for the selected order, 'rorder', 
of the Lie-Backlund transformation being considered.*)
lengthofderivativeindextable=
Length[derivativeindextable];

(*Make a table of the independent variables
corresponding to the derivative indices. In a 
moment these will be used to construct the
derivative expressions used to differentiate the
left and right hand sides of the original rule.*)
derivativeindicessofar=Table["x"<>ToString[jnd],
{ind,1,lengthofderivativeindextable},
{jnd,1,numberindependentvars},
{k,1,derivativeindextable[[ind,jnd]]}];

(*Convert each outer item in the table to a string.*)
derivativeindsstrings=Table[ToString[Flatten[
derivativeindicessofar[[ind]]]],
{ind,1,lengthofderivativeindextable}];

(*Now form the derivative expressions operating
on the left hand term in the original rule. Note
that the derivatives are executed.*)
lefthandtermderivativetable=Table[ToExpression[
"D["<>ToString[InputForm[lefthandterm]]<>","<>
StringReplace[
StringReplace[derivativeindsstrings[[ind]],"{"->""],
"}"->""]<>"]"],
{ind,1,lengthofderivativeindextable}];

(*Now form the derivative expressions operating
on the right hand term in the original rule. Note
that the derivatives are executed.*)
righthandtermderivativetable1=Table[ToExpression[
"D["<>ToString[InputForm[righthandterm]]<>","<>
StringReplace[
StringReplace[derivativeindsstrings[[ind]],"{"->""],
"}"->""]<>"]"],
{ind,1,lengthofderivativeindextable}];

(*Make sure the right hand sides are fully expanded.*)
righthandtermderivativetable=ExpandAll[
righthandtermderivativetable1];

(*Now convert the left and right sides to strings
and combine using -> to form the table of rules
in the form of strings.*)
rulesarrayexpandedtableofstrings=Table[
ToString[InputForm[
lefthandtermderivativetable[[ind]]]]<>"->"<>
ToString[InputForm[righthandtermderivativetable[[ind]]]],
{ind,1,lengthofderivativeindextable}];

(*Convert the strings to expressions.*)
rulesarrayexpandedtable1=ToExpression[
rulesarrayexpandedtableofstrings];

(*Rules for the differential consequences of the ith equation.
The user can view the additional rules by looking at 
rulesarrayexpandedtable1sep98.*)
rulesarrayexpandedtable1sep98=rulesarrayexpandedtable1;

(*Now combine with the original set of rules to make
the expanded set of rules.*)
rulesarrayexpanded=
Join[rulesarray,rulesarrayexpandedtable1];])
(************END*OF*FUNCTION*MakeRulesArray***********)

(**BEGINNING*OF*FUNCTION*InvarianceConditionRulesApplied**)
(*This function generates a table of expressions
corresponding to the terms in the invariance condition
with rules applied, ie, after the original equation
and its differential consequences have been applied
to the invariance condition.*)
InvarianceConditionRulesApplied[
numberindependentvars_,numberdependentvars_,porderoftheequation_,
rorderofinfinitesimals_,xseon_,inputequation_,rulesarray_,internalrules_]:=

(*The function "Module" hides the variable names in 
brackets from the rest of the program.*)
(Module[{termsoftheinvarianceconditiontimestwo,
termsoftheinvarianceconditiontimestwolength,
termsoftheinvarianceconditionsplittable,
termsoftheinvarianceconditionsplittablelength,
termsoftheinvarianceconditionsplittablelengths,
termsoftheinvarianceconditionsplittablestrings,
termsoftheinvarianceconditionsplittablestringetaxsegone,
termsetaxsegoneexpressions,
termsetaxsegoneexpressionswithrules,
termsetaxsegoneexpressionswithrulesstrings,
termsoftheinvarianceconditionsplittablestringetaxsebackin,
termsoftheinvarianceconditionsplittableexpressetaxsebackin,
invarianceconditiontablerulesapplied,
invarianceconditiontablerulesappliedsum1,
numberoftermsinvarsum},

(*Call the function 'InvarianceConditionNoRules' to
generate the raw invariance condition.*) 
InvarianceConditionNoRules[
numberindependentvars,numberdependentvars,porderoftheequation,rorderofinfinitesimals,xseon,inputequation];

(*This is the raw invariance condition in the form
of a table of terms each of which is a product of
y-derivatives and a derivative of one of the unknown
infinitesimals. In a few cases the infinitesimal
derivative appears by itself.*)
termsoftheInvarianceConditionNoRules;

(*Each term in the invariance condition is a product of 
derivatives of the dependent variables multiplied by
a derivative of one of the unknown infinitesimals. 
Eventually we want to isolate the infinitesimals to 
produce the set of determining conditions. This separation
is accomplished using string matching functions.*)
termsoftheinvarianceconditiontimestwo=
2*termsoftheInvarianceConditionNoRules;

(*Find the length of this table.*)
termsoftheinvarianceconditiontimestwolength=
Length[termsoftheinvarianceconditiontimestwo];

(*Separate each factor in each expression.*)
termsoftheinvarianceconditionsplittable=
Table[Cases[termsoftheinvarianceconditiontimestwo[[
indx]],x_],
{indx,1,Length[termsoftheinvarianceconditiontimestwo]}];

(*This table should be the same length as
'termsoftheinvarianceconditiontimestwolength'.*)
termsoftheinvarianceconditionsplittablelength=Length[
termsoftheinvarianceconditionsplittable];

(*Now make a table of lengths of each individual
term in the split up invariance condition, ie, the
number of factors in each term.*)
termsoftheinvarianceconditionsplittablelengths=Table[
Length[termsoftheinvarianceconditionsplittable[[ind]]],
{ind,1,termsoftheinvarianceconditionsplittablelength}];

(*Create a table of strings corresponding to the terms 
in the split table.*)
termsoftheinvarianceconditionsplittablestrings=Table[ToString[
InputForm[termsoftheinvarianceconditionsplittable[[ind,jnd]]]],
{ind,1,termsoftheinvarianceconditionsplittablelength},
{jnd,1,termsoftheinvarianceconditionsplittablelengths[[
ind]]}];

(*Create a two-dimensional table where, if an item in the 
split table matches the strings 'xse' or 'eta' set that 
position to one.*)
termsoftheinvarianceconditionsplittablestringetaxsegone=
Table[If[StringMatchQ[
termsoftheinvarianceconditionsplittablestrings[[ind,jnd]],
"*eta*"]||
StringMatchQ[
termsoftheinvarianceconditionsplittablestrings[[ind,jnd]],
"*xse*"],"1",
termsoftheinvarianceconditionsplittablestrings[[ind,jnd]]],
{ind,1,termsoftheinvarianceconditionsplittablelength},
{jnd,1,
termsoftheinvarianceconditionsplittablelengths[[ind]]}];

(*Convert the previous table to a table of expressions.*)
termsetaxsegoneexpressions=ToExpression[
termsoftheinvarianceconditionsplittablestringetaxsegone];

rulesarraylengthmmax=Length[rulesarray];

If[rulesarraylengthmmax>=1,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=1,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=1,
rulesarrayvector1=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=1,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1]];


If[rulesarraylengthmmax>=2,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=2,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=2,
rulesarrayvector2=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=2,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2]];


If[rulesarraylengthmmax>=3,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=3,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=3,
rulesarrayvector3=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=3,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3]];


If[rulesarraylengthmmax>=4,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=4,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=4,
rulesarrayvector4=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=4,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4]];


If[rulesarraylengthmmax>=5,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=5,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=5,
rulesarrayvector5=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=5,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5]];


If[rulesarraylengthmmax>=6,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=6,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=6,
rulesarrayvector6=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=6,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6]];


If[rulesarraylengthmmax>=7,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=7,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=7,
rulesarrayvector7=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=7,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7]];


If[rulesarraylengthmmax>=8,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=8,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=8,
rulesarrayvector8=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=8,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7,rulesarrayvector8]];


If[rulesarraylengthmmax>=9,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=9,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=9,
rulesarrayvector9=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=9,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7,rulesarrayvector8,
rulesarrayvector9]];


If[rulesarraylengthmmax>=10,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=10,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=10,
rulesarrayvector10=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=10,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7,rulesarrayvector8,
rulesarrayvector9,rulesarrayvector10]];


If[rulesarraylengthmmax>=11,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=11,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=11,
rulesarrayvector11=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=11,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7,rulesarrayvector8,
rulesarrayvector9,rulesarrayvector10,
rulesarrayvector11]];


If[rulesarraylengthmmax>=12,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=12,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=12,
rulesarrayvector12=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=12,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7,rulesarrayvector8,
rulesarrayvector9,rulesarrayvector10,
rulesarrayvector11,rulesarrayvector12]];

If[rulesarraylengthmmax>=13,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=13,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=13,
rulesarrayvector13=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=13,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7,rulesarrayvector8,
rulesarrayvector9,rulesarrayvector10,
rulesarrayvector11,rulesarrayvector12,
rulesarrayvector13]];

If[rulesarraylengthmmax>=14,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=14,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=14,
rulesarrayvector14=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=14,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7,rulesarrayvector8,
rulesarrayvector9,rulesarrayvector10,
rulesarrayvector11,rulesarrayvector12,
rulesarrayvector13,rulesarrayvector14]];

If[rulesarraylengthmmax>=15,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=15,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=15,
rulesarrayvector15=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=15,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7,rulesarrayvector8,
rulesarrayvector9,rulesarrayvector10,
rulesarrayvector11,rulesarrayvector12,
rulesarrayvector13,rulesarrayvector14,
rulesarrayvector15]];

If[rulesarraylengthmmax>=16,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=16,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=16,
rulesarrayvector16=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=16,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7,rulesarrayvector8,
rulesarrayvector9,rulesarrayvector10,
rulesarrayvector11,rulesarrayvector12,
rulesarrayvector13,rulesarrayvector14,
rulesarrayvector15,rulesarrayvector16]];

If[rulesarraylengthmmax>=17,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=17,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=17,
rulesarrayvector17=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=17,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7,rulesarrayvector8,
rulesarrayvector9,rulesarrayvector10,
rulesarrayvector11,rulesarrayvector12,
rulesarrayvector13,rulesarrayvector14,
rulesarrayvector15,rulesarrayvector16,
rulesarrayvector17]];

If[rulesarraylengthmmax>=18,
MakeRulesArray[numberindependentvars,numberdependentvars,mequindexsep98=18,rorderofinfinitesimals,rulesarray]];

If[rulesarraylengthmmax>=18,
rulesarrayvector18=
rulesarrayexpandedtable1sep98];

If[rulesarraylengthmmax>=18,
rulesarraysep98expanded1=Join[
rulesarray,
rulesarrayvector1,rulesarrayvector2,
rulesarrayvector3,rulesarrayvector4,
rulesarrayvector5,rulesarrayvector6,
rulesarrayvector7,rulesarrayvector8,
rulesarrayvector9,rulesarrayvector10,
rulesarrayvector11,rulesarrayvector12,
rulesarrayvector13,rulesarrayvector14,
rulesarrayvector15,rulesarrayvector16,
rulesarrayvector17,rulesarrayvector18]];

rulesarrayexpanded=rulesarraysep98expanded1;

(*Apply the rules in 'rulesarrayout' to each expression in
the modified invariance condition. The condition has been 
modified to avoid applying the rules to the independent 
variables which the infinitesimals depend on.*)
If[internalrules==1,
termsetaxsegoneexpressionswithrules=
termsetaxsegoneexpressions//.rulesarrayexpanded,
termsetaxsegoneexpressionswithrules=
termsetaxsegoneexpressions//.rulesarray];

(*Produce a corresponding table of strings.*)
termsetaxsegoneexpressionswithrulesstrings=Table[
ToString[InputForm[termsetaxsegoneexpressionswithrules[[
ind,jnd]]]],
{ind,1,termsoftheinvarianceconditionsplittablelength},
{jnd,1,
termsoftheinvarianceconditionsplittablelengths[[ind]]}];

(*Put xse's and eta's back in.*)
termsoftheinvarianceconditionsplittablestringetaxsebackin=
Table[If[StringMatchQ[
termsoftheinvarianceconditionsplittablestrings[[ind,jnd]],
"*eta*"]||
StringMatchQ[
termsoftheinvarianceconditionsplittablestrings[[ind,jnd]],
"*xse*"],
termsoftheinvarianceconditionsplittablestrings[[ind,jnd]],
termsetaxsegoneexpressionswithrulesstrings[[ind,jnd]]],
{ind,1,termsoftheinvarianceconditionsplittablelength},
{jnd,1,
termsoftheinvarianceconditionsplittablelengths[[ind]]}];

(*Now recreate the invariance condition with rules applied.*)
termsoftheinvarianceconditionsplittableexpressetaxsebackin=
ToExpression[
termsoftheinvarianceconditionsplittablestringetaxsebackin];

(*Multiply the factors in each term back together.*)
invarianceconditiontablerulesapplied=Table[
Apply[Times,
termsoftheinvarianceconditionsplittableexpressetaxsebackin[[
ind]]],
{ind,1,termsoftheinvarianceconditionsplittablelength}];

(*Here is the invariance condition with rules applied in 
the form of a fully expanded sum.Note the factor of 2
multiplying each term.*)
invarianceconditiontablerulesappliedsum1=
ExpandAll[Apply[Plus,
invarianceconditiontablerulesapplied]];

(*The actual invariance condition sum 
without the factor of two.*)
invarianceconditiontablerulesappliedsum=
(1/2)*invarianceconditiontablerulesappliedsum1;

(*The number of terms in the sum.*)
numberoftermsinvarsum=Length[
invarianceconditiontablerulesappliedsum1];

(*Here is the invariance condition with rules applied 
sorted into a table of terms.*)
termsoftheInvarianceConditionRulesApplied1=Table[
invarianceconditiontablerulesappliedsum1[[ind]],
{ind,1,numberoftermsinvarsum}];

(*The actual invariance condition table, 
without the factor of two.*)
termsoftheInvarianceConditionRulesApplied=
(1/2)*termsoftheInvarianceConditionRulesApplied1;])
(*****END*OF*FUNCTION*InvarianceConditionRulesApplied********)

(****BEGINNING*OF*FUNCTION*FindDeterminingEquations*******)
(*This function generates a table of expressions
corresponding to the determining equations of the group.*)
FindDeterminingEquations[
independentvariablesstrings_,dependentvariablesstrings_,frozenstrings_,porderoftheequation_,
rorderofinfinitesimals_,xseon_,inputequationuservariables_,rulesarrayuservariables_,internalrules_]:=

(*The function "Module" hides the variable names in 
brackets from the rest of the program.*)
(Module[{termsoftheInvarianceConditionRulesApplied1length,
termsoftheInvarianceConditionRulesApplied1splittable,
variablestringtable1,infinitesimalindependentvariables,
variablestringtable2,allvariables,freeyderivativevariables,
freeyderivativevariableslength,freeyderivativevariablesstrings,
lengthterms1table,ind,
termsoftheInvarianceConditionRulesApplied1splittablestrings,
truthtablefreevariables,truthtablenumbers,
freevariablelocations,freeyderivativetable,
infinitesimaltermtable1,infinitesimalproducttable,
ytermstablngth,ytermsmatchtable,ytermsmatchtableoflengths,
infinitesimaltablesums,infinitesimaltablesumsmatchtbl,
infinitesimaltablesumsmatchtblfirstindices,
finalinfinitesimalsumindices,
determiningequationstable1,
infinitesimalindepvariabletable,
infinitesimalindepvariabletablelength,
dummyvariabletable,
dummyvariabletableofrulesstrings,
xseclearstring1,xseclearstring,
etaclearstring1,etaclearstring,determiningequations,
zdeterminingequationstable,
zdeterminingequationstablelengthof,
determiningequationstable,
determiningequationstablelengthof,
zdeterminingequationstableofstrings,
determiningequationstableofstrings,numberindependentvars,numberdependentvars,
frozenstringslength,inputequationindependentvariablerules,
inputequationdependentvariablerules,inputequationvariablerules,
frozenstringsnewsymbols,frozenstringrules,
inputequationbeforefrozenstrings,inputequation,
rulesarraybeforefrozenstrings,
rulesarraystrings,rulesarray,ztableofrulesxytoz,
independentvariables,dependentvariables,
dummyvariabletableofrulesstringsztoxy,
dummyvariabletableofrulesstringsxytoinput,
inputequationdependentvariablerulesytoinput,
inputequationindependentvariablerulesxtoinput,
inputequationvariablerulesxytoinput,
inputequationvariablerulesxytoinputreverse},

(*Convert the list of independent variables to an expression*)
independentvariables=ToExpression[independentvariablesstrings];

(*Convert the list of dependent variables to an expression*)
dependentvariables=ToExpression[dependentvariablesstrings];

(*How many independentvariables are there?*)
numberindependentvars=Length[independentvariables];

(*How many dependentvariables are there?*)
numberdependentvars=Length[dependentvariables];

(*How many frozenstrings (function and/or constant names) appear in the equations?*)
frozenstringslength=Length[frozenstrings];

(*Construct a set of rules for converting the input independent variables
to generic variables, x1, x2, x3,...,xn.*)
inputequationindependentvariablerules=
Table[ToString[
independentvariables[[jjnd]]]->"x"<>ToString[jjnd],{jjnd,1,numberindependentvars}];

(*Construct a set of rules for converting the input dependent variables
to generic variables, y1, y2, y3,...,ym.*)
inputequationdependentvariablerules=
Table[ToString[
dependentvariables[[iind]]]->"y"<>ToString[iind],{iind,1,numberdependentvars}];

(*Join these two sets of rules.*)
inputequationvariablerules=
Join[inputequationindependentvariablerules,
inputequationdependentvariablerules];

(*Convert the input frozenstrings to strings in case the user fails to 
enter these as strings. The reason the table of frozenstrings is required 
is that when the user defined variables are converted to generic variables, 
pieces of function or constant names may also be converted 
if the characters happen to be common to the two names. Here a table of incorrect parameter names is actually created.*)
frozenstringsnewsymbols=
StringReplace[frozenstrings,inputequationvariablerules];

(*This command sets up a table of rules converting incorrect 
parameter names to correct ones.*)
frozenstringrules=
Table[frozenstringsnewsymbols[[kknd]]->ToString[frozenstrings[[kknd]]],{kknd,1,frozenstringslength}];

(*Now change input user variable names to generic names.*)
inputequationbeforefrozenstrings=
StringReplace[inputequationuservariables,inputequationvariablerules];

(*Correct any errors that may have been created in the parameter names. 
This is the input equation in string form and in terms of generic variable 
names that is now used throughout the program.*)
inputequation=
StringReplace[inputequationbeforefrozenstrings,frozenstringrules];

(*Be able to look at the equation in generic variables*)
equationgenericvariables=inputequation;

(*The same replacement of user defined names needs to be carried 
out in the table of rules input by the user.*)
rulesarraybeforefrozenstrings=
StringReplace[rulesarrayuservariables,inputequationvariablerules];

(*Similarly errors produced in the parameter names by the replacement
process need to be corrected.*)
rulesarraystrings=
StringReplace[rulesarraybeforefrozenstrings,frozenstringrules];

(*Convert the rules put in by the user as strings to expressions
that will be used by the program.*)
rulesarray=ToExpression[rulesarraystrings];

(*Clear the variable names for the 'xse' infinitesimals.First construct
a table of strings with the infinitesimal names.*)
xseclearstring1=Table["xse"<>ToString[j],{j,1,numberindependentvars}];

(*Remove the curly brackets.*)
xseclearstring=StringReplace[StringReplace[
ToString[xseclearstring1],"{"->""],"}"->""];

(*Now clear the 'xse' variable names.*)
ToExpression["ClearAll["<>xseclearstring<>"]"];

(*Clear the variable names for the 'eta' infinitesimals.First construct
a table of strings with the infinitesimal names.*)
etaclearstring1=Table["eta"<>ToString[i],{i,1,numberdependentvars}];

(*Remove the curly brackets.*)
etaclearstring=StringReplace[StringReplace[
ToString[etaclearstring1],"{"->""],"}"->""];

(*Now clear the 'eta' variable names.*)
ToExpression["ClearAll["<>etaclearstring<>"]"];

(*Call the function 'InvarianceConditionRulesApplied' to
produce a table of terms of the invariance condition. Note
that each term has a factor of 2 out in front. This is so
that the function 'Cases' will properly split each term
into a table of factors. Otherwise terms that happen not
to have a multiplying factor out in front get separated at
too low a level. For example y1[x1,x2] would be split into
{y1,[x1,x2]}. Whereas 2*y1[x1,x2] would be split up as
{2,y1[x1,x2]} which is what is wanted.*)
InvarianceConditionRulesApplied[
numberindependentvars,numberdependentvars,porderoftheequation,rorderofinfinitesimals,xseon,inputequation,rulesarray,internalrules];

(*The invariance condition with rules applied 
sorted into a table of terms, with the factor of two.*)
termsoftheInvarianceConditionRulesApplied1;

(*Number of terms in this table.*)
termsoftheInvarianceConditionRulesApplied1length=Length[
termsoftheInvarianceConditionRulesApplied1];

(*Split each term into a table of factors.*)
termsoftheInvarianceConditionRulesApplied1splittable=
Table[Cases[
termsoftheInvarianceConditionRulesApplied1[[indx]],x_],
{indx,1,
termsoftheInvarianceConditionRulesApplied1length}];

(*Now generate a table of the independent
variables upon which the infinitesimals are
assumed to depend.*)
GenerateVariableTable[numberindependentvars,numberdependentvars,rorderofinfinitesimals];

(*This is the result in the form of a table of strings.*)
variablestringtable1=variablestringtable;

(*Convert each item to an expression.*)
infinitesimalindependentvariables=ToExpression[
variablestringtable1];

(*Now generate a table of all the y-derivatives 
that can appear in the problem.*)
GenerateVariableTable[numberindependentvars,numberdependentvars,rorderofinfinitesimals+porderoftheequation];

(*This is the resulting table.*)
variablestringtable2=variablestringtable;

(*Convert to a table of expressions.*)
allvariables=ToExpression[
variablestringtable2];

(*Make a table of the y-derivatives which are free 
variables.These are the y-derivatives which the 
infinitesimals don't depend on.*)
freeyderivativevariables=Complement[
allvariables,infinitesimalindependentvariables];

(*The number of items in this table.*)
freeyderivativevariableslength=
Length[freeyderivativevariables];

(*Convert the list of free variables to a list of strings.*)
freeyderivativevariablesstrings=
Table[ToString[InputForm[freeyderivativevariables[[knd]]]],
{knd,1,freeyderivativevariableslength}];

(*Make a table of the number of factors in each term of
the invariance condition with rules applied.*)
lengthterms1table=Table[Length[
termsoftheInvarianceConditionRulesApplied1splittable[[ind]]],
{ind,1,
termsoftheInvarianceConditionRulesApplied1length}];

(*Convert the table of invariance condition term factors to
a table of strings.*)
termsoftheInvarianceConditionRulesApplied1splittablestrings=
Table[ToString[InputForm[
termsoftheInvarianceConditionRulesApplied1splittable[[
ind,jnd]]]],
{ind,1,
termsoftheInvarianceConditionRulesApplied1length},
{jnd,1,lengthterms1table[[ind]]}];

(*Now create a three dimensional true-false table matching each
term in the table of invariance condition factors with each term
in the table of free y-derivative variables.*)
truthtablefreevariables=Table[StringMatchQ[
termsoftheInvarianceConditionRulesApplied1splittablestrings[[
ind,jnd]],
"*"<>freeyderivativevariablesstrings[[knd]]<>"*"],
{ind,1,
termsoftheInvarianceConditionRulesApplied1length},
{jnd,1,lengthterms1table[[ind]]},
{knd,1,freeyderivativevariableslength}];

(*Convert this table to one's and zero's with one
corresponding to 'True'.*)
truthtablenumbers=Table[
If[truthtablefreevariables[[ind,jnd,knd]],1,0],
{ind,1,
termsoftheInvarianceConditionRulesApplied1length},
{jnd,1,lengthterms1table[[ind]]},
{knd,1,freeyderivativevariableslength}];

(*Now sum this three-dimensional table over the depth 
of the number of free derivative terms. This puts a
one at each position of the invariance condition 
corresponding to any of the free derivative terms.*)
freevariablelocations=
Table[Apply[Plus,truthtablenumbers[[ind,jnd]]],
{ind,1,
termsoftheInvarianceConditionRulesApplied1length},
{jnd,1,lengthterms1table[[ind]]}];

(*Now split the table of invariance condition terms into 
two tables. First the table of free y-derivative terms.*)
freeyderivativetable=Table[
If[freevariablelocations[[ind,jnd]]==1,
termsoftheInvarianceConditionRulesApplied1splittable[[
ind,jnd]],
1],{ind,1,
termsoftheInvarianceConditionRulesApplied1length},
{jnd,1,lengthterms1table[[ind]]}];

(*Second the table of infinitesimal terms.*)
infinitesimaltermtable1=Table[
If[freevariablelocations[[ind,jnd]]==0,
termsoftheInvarianceConditionRulesApplied1splittable[[
ind,jnd]],
1],{ind,1,
termsoftheInvarianceConditionRulesApplied1length},
{jnd,1,lengthterms1table[[ind]]}];

(*Create a one-dimensional table of y-derivative terms
multiplied together.*)
ytermstable=Table[Apply[Times,
freeyderivativetable[[ind]]],
{ind,1,
termsoftheInvarianceConditionRulesApplied1length}];

(*Create a one-dimensional table of infinitesimal
terms multiplied together.*)
infinitesimalproducttable=Table[Apply[Times,
infinitesimaltermtable1[[ind]]],
{ind,1,
termsoftheInvarianceConditionRulesApplied1length}];

(*Remove the factor of two inserted earlier. Let the 
tables 'ytermstable' and 'infinitesimaltable' be 
available as output, ie, don't include them within
the set of names protected by 'Module'.*)
infinitesimaltable=(1/2)*infinitesimalproducttable;

(*Check the number of terms in the table.*)
ytermstablngth=Length[ytermstable];

(*Identify the positions of yderivative terms in 
the table that match.*)
ytermsmatchtable=Table[
Position[
ytermstable,ytermstable[[ind]],{1}],
{ind,1,ytermstablngth}];

(*Make a table giving the number of times a given
y-derivative expression occurs.*)
ytermsmatchtableoflengths=Table[
Length[ytermsmatchtable[[ind]]],
{ind,1,ytermstablngth}];

(*Now sum the infinitesimal terms corresponding 
to matching y-derivative expressions.*)
infinitesimaltablesums=Flatten[
Table[
Sum[infinitesimaltable[[ytermsmatchtable[[ind,jnd]]]],
{jnd,1,ytermsmatchtableoflengths[[ind]]}],
{ind,1,ytermstablngth}]];

(*Produce a table containing indices of matching 
infinitesimal sums.*)
infinitesimaltablesumsmatchtbl=Table[
Position[
infinitesimaltablesums,infinitesimaltablesums[[ind]],{1}],
{ind,1,ytermstablngth}];

(*Pick off the first index.*)
infinitesimaltablesumsmatchtblfirstindices=Table[
infinitesimaltablesumsmatchtbl[[ind,1]],
{ind,1,ytermstablngth}];

(*Now use the function 'Union' to create a table of
indices where no two indices are the same.*)
finalinfinitesimalsumindices=Union[
infinitesimaltablesumsmatchtblfirstindices];
finalinfinitesimalsumindiceslength=Length[
finalinfinitesimalsumindices];

(*Now use the indices identified in the previous
step to pick off the unique terms in the table of
infinitesimal sums. These are the left-hand sides 
of the sought after determining equations of the 
group. The table 'determiningequationstable' is 
presented as output and therefore is not included 
with the names protected by 'Module'.*)
determiningequationstable1=
Flatten[
Table[
infinitesimaltablesums[[
finalinfinitesimalsumindices[[ind]]]],
{ind,1,finalinfinitesimalsumindiceslength}]];

(*The number of determining equations.*)
determiningequationstable1length=
Length[determiningequationstable1];

(*Convert each determining equation to a string and
append '==0' to each string.*)
determiningequationstable2=
Table[ToString[InputForm[
determiningequationstable1[[ind]]]]<>"==0",
{ind,1,determiningequationstable1length}];

(*Use the same set of indices to pick off the 
corresponding y-derivative factors. The product
of these two tables summed is the final form
of the invariance condition.The table
'yderivfactortable' is presented as
output and therefore is not included with the
names protected by 'Module'.*)
yderivfactortable=Flatten[
Table[
ytermstable[[
finalinfinitesimalsumindices[[ind]]]],
{ind,1,finalinfinitesimalsumindiceslength}]];

(*Convert each term in the determining equations table 
back to an expression.*)
determiningequationstable=
Table[ToExpression[determiningequationstable2[[ind]]],
{ind,1,determiningequationstable1length}];

(*The number of determining equations - again.*)
determiningequationstablelengthof=
Length[determiningequationstable];

(*Convert the determining equations table to a table of strings and
present the primary output of the function in that form.*)
determiningequationstableofstrings=
Table[ToString[InputForm[determiningequationstable[[ind]]]],
{ind,1,determiningequationstablelengthof}];

(*The variable name in the last step is just too long
and unwieldy to expect the user to get used to. Here is
a shorter name.This is the main output of the program.*)
determiningequations=determiningequationstableofstrings;

(*For convenience express the determining equations 
in terms of a set of dummy variables.*)
GenerateVariableTable[numberindependentvars,numberdependentvars,rorderofinfinitesimals];

(*Convert these strings to expressions.*)
infinitesimalindepvariabletable=
ToExpression[variablestringtable];

(*The number of dummy variables that need to be 
generated.*)
infinitesimalindepvariabletablelength=
Length[infinitesimalindepvariabletable];

(*Now generate the required number of dummy 
variables.*)
dummyvariabletable=
Table[ToExpression["z"<>ToString[indx]],
{indx,1,infinitesimalindepvariabletablelength}];

(*Generate the needed table of transformation rules,
first as strings.*)
dummyvariabletableofrulesstrings=
Table[ToString[
InputForm[
infinitesimalindepvariabletable[[indx]]]]<>
"->"<>ToString[dummyvariabletable[[indx]]],
{indx,1,infinitesimalindepvariabletablelength}];

(*Convert the strings to expressions.*)
ztableofrulesxytoz=
ToExpression[dummyvariabletableofrulesstrings];

(*Be able to view this table*)
ztableofrulesxy=ztableofrulesxytoz;

(*Now make the change to dummy variables.
This is the final output of the program.*)
zdeterminingequationstable=
determiningequationstable/.ztableofrulesxytoz;

(*The number of determining equations - yet again.*)
zdeterminingequationstablelengthof=
Length[zdeterminingequationstable];

(*Convert the above table to a table of strings and
present the primary output of the function in that form.*)
zdeterminingequationstableofstrings=
Table[ToString[InputForm[zdeterminingequationstable[[ind]]]],
{ind,1,zdeterminingequationstablelengthof}];

(*The variable name in the last step is just too long
and unwieldy to expect the user to get used to. Here is
a shorter z-name.This is a convenient form of the main output
of the program.*)
zdeterminingequations=zdeterminingequationstableofstrings;

(*Produce a table of rules that takes z variables to xy variables*)
dummyvariabletableofrulesstringsztoxy=
Table[ToString[dummyvariabletable[[indx]]]<>
"->"<>ToString[
InputForm[
infinitesimalindepvariabletable[[indx]]]],
{indx,1,infinitesimalindepvariabletablelength}];

(*Construct a set of rules for converting generic independent variables, x1, x2, x3,...,xn to input independent variables.*)
inputequationindependentvariablerulesxtoinput=
Table["x"<>ToString[jjnd]->ToString[
independentvariables[[jjnd]]],{jjnd,1,numberindependentvars}];

(*Construct a set of rules for converting generic dependent variables, y1, y2, y3,...,ym to input dependent variables.*)
inputequationdependentvariablerulesytoinput=
Table["y"<>ToString[iind]->ToString[
dependentvariables[[iind]]],{iind,1,numberdependentvars}];

(*Join these two sets of rules.*)
inputequationvariablerulesxytoinput=
Join[inputequationindependentvariablerulesxtoinput,
inputequationdependentvariablerulesytoinput];

(*If the number of independent or dependent variables is ten or more
the replacement will change the left two characters in y10, y11, etc 
with the user variable corresponding to y1. This happens because the
replacement begins with the first element in the replacement rules for 
each entry in the list dummyvariabletableofrulesstringsztoxy. To avoid 
this problem we reverse the order of the replacement rules in 
inputequationvariablerulesxytoinput so that y10, y11, etc are replaced 
before the y1 is replaced. The problem was encountered when analyzing 
the second order form of the restricted Euler equation which has 11 
dependent variables including Q and R. This has no effect on the operation 
of the package except that the list ztableofrules would have had incorrect 
entries for variables beyond the ninth. This change was made September 8, 2001.*)
inputequationvariablerulesxytoinputreverse=
Reverse[inputequationvariablerulesxytoinput];

(*Create a table of rules that take xy variables to input variables.*)
dummyvariabletableofrulesstringsxytoinput=
StringReplace[
dummyvariabletableofrulesstringsztoxy,
inputequationvariablerulesxytoinputreverse];

(*Convert the strings to expressions.*)
ztableofrules=
ToExpression[dummyvariabletableofrulesstringsxytoinput];

(*Inform the user that the output is ready.*)
Print["FindDeterminingEquations has finished executing. You can look at the output in the table zdeterminingequations. Each entry in this table is a determining equation in string format expressed in terms of z-variables. Rules for converting between z-variables and conventional variables are contained in the table ztableofrules. To view the determining equations in terms of conventional variables use the command ToExpression[zdeterminingequations]/.ztableofrules. There are two other items the user may wish to look at; the equation converted to generic (x1,x2,...,y1,y2,...) variables is designated equationgenericvariables and the various derivatives of the equation that appear in the invariance condition can be viewed in the table invarconditiontable. Rules for converting between z-variables and generic variables are contained in the table ztableofrulesxy."];])
(******END*OF*FUNCTION*FindDeterminingEquations*******)

(****** *BEGINNING*OF*FUNCTION*SolveDeterminingEquations****)
(*This function takes the determining equations contained in 
the table 'zdeterminingequationstable' and attempts to 
find a solution in the form of a power series. This will usually
capture many but not all of the solutions of the determining 
equations which can contain arbitrary functions. The presence
of such functions is usually indicated when the series solution 
doesn't truncate.The output of the function is presented in two
tables 'xseinfinitesimals' and 'etainfinitesimals'containing
algebraic expressions for the infinitesimals ordered according to
the conventional ordering of independent variables and dependent
variables respectively.*)
SolveDeterminingEquations[
independentvariablesstrings_,dependentvariablesstrings_,rorderofinfinitesimals_,xseon_,
zdeterminingequationstableofstrings_,polyorder_]:=

(*The function 'Module' protects variable names from the 
outside world.*)
(Module[{numberofvariables,zvariables,kindices,
kindicessumstring,powersofzs1,powersproductstring,
indicestablestring1,indicestablestring,powerofzs1string,
powertable1,powertablenull1,powertablenull,
lengthpowertable,zvariablesundefinedstring,
powertable,xsecoefficienttable,etacoefficienttable,
xseproductstep1,xsefunctionsumtable,
tableofxxsefunctionsstrings,
tableofxxsefunctions1,tableofxxsefunctions,
etaproductstep1,etafunctionsumtable,
tableofeetafunctionsstrings,
tableofeetafunctions,
zdeterminingequationstablelength,
zdeterminingequationstablestrings,
zdeterminingequationstablexxseeetastrings,
zdeterminingequationstablexxseeetastringsnoequals,
zdeterminingequationstablexxseeetaexpressions,
expandedzdeterminingequationstable,tableofcoefficients,
tableofcoefficientsall,
tableofcoefficientslength,tableofcoefficientsstrings,
unknowns1,unknowns,rulestable1,rulestable,i,j,ind,
xseinfinitesimals1,etainfinitesimals1,
zdeterminingequationstable,
xseclearstring1,xseclearstring,
etaclearstring1,etaclearstring,
xseinfinitesimalstrings,
etainfinitesimalstrings,
allunknowns1a,allunknowns1,groupparametersonesandzeros,
groupparameterswithzeros,groupparameters1,
parameterlength,infinitesimalgroups1,numberindependentvars,numberdependentvars,order,
independentvariables,dependentvariables},

(*Convert the list of independent variables to an expression*)
independentvariables=ToExpression[independentvariablesstrings];

(*Convert the list of dependent variables to an expression*)
dependentvariables=ToExpression[dependentvariablesstrings];

(*order is a protected name for the trial polynomial order.*)
order=polyorder;

(*The number of independent variables*)
numberindependentvars=Length[independentvariables];

(*The number of dependent variables*)
numberdependentvars=Length[dependentvariables];

(*Clear the variable names for the 'xse' infinitesimals. First construct
a table of strings with the infinitesimal names.*)
xseclearstring1=Table["xse"<>ToString[j],{j,1,numberindependentvars}];

(*Remove the curly brackets.*)
xseclearstring=StringReplace[StringReplace[
ToString[xseclearstring1],"{"->""],"}"->""];

(*Now clear the 'xse' variable names.*)
ToExpression["ClearAll["<>xseclearstring<>"]"];

(*Clear the variable names for the 'eta' infinitesimals. First construct
a table of strings with the infinitesimal names.*)
etaclearstring1=Table["eta"<>ToString[i],{i,1,numberdependentvars}];

(*Remove the curly brackets.*)
etaclearstring=StringReplace[StringReplace[
ToString[etaclearstring1],"{"->""],"}"->""];

(*Now clear the 'eta' variable names.*)
ToExpression["ClearAll["<>etaclearstring<>"]"];

(*First, convert the determining equation strings 
to expressions.*)
zdeterminingequationstable=
ToExpression[zdeterminingequationstableofstrings];

(*Call this function to generate a list of the
independent variables, dependent variables and
derivatives of dependent variables. The only
purpose is to determine the number of z variables.*)
GenerateVariableTable[numberindependentvars,numberdependentvars,rorderofinfinitesimals];

(*This is the output list of variables.*)
variablestringtable;

(*This is the number of z-variables.*)
numberofvariables=Length[variablestringtable];

(*Generate the list of z-variables.*)
zvariables=Table[ToExpression["z"<>ToString[ind]],
{ind,1,numberofvariables}];

(*Generate a list of k-indices to be associated with
a polynomial expansion of the infinitesimals in z-variables.*)
kindices=Table[ToExpression["k"<>ToString[ind]],
{ind,1,numberofvariables}];

(*Create a sum of the k-indices to be used a couple of lines
down as a test sum needed to create the polynomial expansions 
of the infinitesimals.*)
kindicessumstring=ToString[Apply[Plus,kindices]];

(*Now create a table of powers of the z-variables.*)
powersofzs1=Table["("<>
ToString[zvariables[[ind]]]<>"^"<>
ToString[kindices[[ind]]]<>")",{ind,1,numberofvariables}];

(*Multiply the various powers together.*)
powersproductstring=ToString[Apply[Times,powersofzs1]];

(*Create the table of index ranges needed in the polynomial expansion.*)
indicestablestring1=ToString[Table["{"<>
ToString[kindices[[numberofvariables+1-ind]]]<>","<>
"0"<>","<>ToString[order]<>"}",{ind,1,numberofvariables}]];

(*Convert the above table of index ranges to a string.*)
indicestablestring=StringReplace[
StringReplace[indicestablestring1,"{{"->"{"],"}}"->"}"];

(*Create a string which, when converted to an expression, will 
create powers of z and where the sum of powers exceeds 'order'
generate a 'Null'.*)
powerofzs1string=
"Table["<>
"If["<>
kindicessumstring<>"<="<>ToString[order]<>","<>
powersproductstring<>"],"<>
indicestablestring<>"]";

(*Flatten the above table.*)
powertable1=Flatten[ToExpression[powerofzs1string]];

(*Create a table of 'Nulls'.*)
powertablenull1="Table["<>
"If["<>
kindicessumstring<>"<=-1"<>","<>
powersproductstring<>"],"<>
indicestablestring<>"]";

(*Flatten this table of 'Null's.*)
powertablenull=Flatten[ToExpression[powertablenull1]];

(*Form a table of powers of z's by taking the Complement of the above
two tables.*)
powertable=Complement[powertable1,powertablenull];

(*Find the length of this table of powers. This is the number
of terms in the polynomial expansions to be generated below.*)
lengthpowertable=Length[powertable];

(*Create a table of coefficients for the 'xse' infinitesimals. If
xseon==0 set the xsecoefficients to zero.*)
xsecoefficienttable=
Table[If[xseon==1,ToExpression["a"<>ToString[j]<>ToString[ind-1]],0],
{j,1,numberindependentvars},{ind,1,lengthpowertable}];

(*Produce a table of expressions of products of the 'xse' coefficients
multiplied by the powers of z.*)
xseproductstep1=Table[xsecoefficienttable[[j,ind]]*
powertable[[ind]],
{j,1,numberindependentvars},{ind,1,lengthpowertable}];

(*Sum the terms in the above tables to produce the 
right hand sides for the 'xse' functions.*)
xsefunctionsumtable=
Table[Apply[Plus,xseproductstep1[[j]]],
{j,1,numberindependentvars}];

(*Create the string to be used as undefined z-variables
in the infinitesimal polynomial functions to be generated 
below.*)
zvariablesundefinedstring=
StringReplace[
StringReplace[
ToString[Table["z"<>ToString[ind]<>"_",
{ind,1,numberofvariables}]],
"{"->""],"}"->""];

(*Now create a table of string function definitions for the 'xse' 
infinitesimals expressed as polynomials with powers up 
to 'order'.*)
tableofxxsefunctionsstrings=Table["xxse"<>ToString[j]<>"["<>
zvariablesundefinedstring<>"]"<>"="<>
ToString[InputForm[xsefunctionsumtable[[j]]]],{j,1,numberindependentvars}];

(*Convert the 'xse' strings to expressions.*)
tableofxxsefunctions1=ToExpression[tableofxxsefunctionsstrings];

(*Set these functions to zero if xseon is not equal to one.*)
tableofxxsefunctions=
Table[If[xseon==1,tableofxxsefunctions1[[j]],0],{j,1,numberindependentvars}];

(*Now create a table of coefficients for the 'eta' infinitesimals.*)
etacoefficienttable=
Table[ToExpression["b"<>ToString[i]<>ToString[ind-1]],
{i,1,numberdependentvars},{ind,1,lengthpowertable}];

(*Produce a table of expressions of products of the 'eta' coefficients
multiplied by the powers of z.*)
etaproductstep1=Table[etacoefficienttable[[i,ind]]*
powertable[[ind]],
{i,1,numberdependentvars},{ind,1,lengthpowertable}];

(*Sum the terms in the above tables to produce the 
right hand sides for the 'eta' functions.*)
etafunctionsumtable=
Table[Apply[Plus,etaproductstep1[[i]]],
{i,1,numberdependentvars}];

(*Now create a table of string function definitions for the 'eta' 
infinitesimals expressed as polynomials with powers up 
to 'order'.*)
tableofeetafunctionsstrings=Table["eeta"<>ToString[i]<>"["<>
zvariablesundefinedstring<>"]"<>"="<>
ToString[InputForm[etafunctionsumtable[[i]]]],{i,1,numberdependentvars}];

(*Convert the 'eta' strings to expressions.*)
tableofeetafunctions=ToExpression[tableofeetafunctionsstrings];

(*The number of determining equations.*)
zdeterminingequationstablelength=
Length[zdeterminingequationstable];

(*Convert the determining equations to strings.*)
zdeterminingequationstablestrings=
Table[ToString[InputForm[
zdeterminingequationstable[[ind]]]],
{ind,1,zdeterminingequationstablelength}];

(*Change 'xse' to 'xxse' and 'eta' to 'eeta'. The 
reason for this is that the polynomial expansions
are defined so as to be immediately evaluated. 
Changing names slightly prevents the original table
of determining equations from being modified.*)
zdeterminingequationstablexxseeetastrings=
StringReplace[StringReplace[
zdeterminingequationstablestrings,"xse"->"xxse"],
"eta"->"eeta"];

(*Remove the double equals from each of the determining equations.*)
zdeterminingequationstablexxseeetastringsnoequals=
StringReplace[
zdeterminingequationstablexxseeetastrings,"== 0"->""];

(*Now convert the determining equations strings to 
expressions. It is important that this step come after
the polynomial expressions for the 'xse' and 'eta' 
functions. This is because the 'xse' and 'eta' functions
are defined so as to be evaluated immediately. Therefore
the substitution of these functions into the determining
equations will occur immediately.*)
zdeterminingequationstablexxseeetaexpressions=
Table[ToExpression[
zdeterminingequationstablexxseeetastringsnoequals[[ind]]],
{ind,1,zdeterminingequationstablelength}];

(*Now substitute the expansions into the table of 
determining equations with the ==0 removed.*)
expandedzdeterminingequationstable=
ExpandAll[
zdeterminingequationstablexxseeetaexpressions];

(*Now pick off the coefficients of the z-variables 
in the determining equations.*)
tableofcoefficientsall=
Flatten[CoefficientList[
expandedzdeterminingequationstable,zvariables]];

(*The list of coefficients typically involves a large
number of items that are simply zero. The function Union
removes redundant members from the list.*)
tableofcoefficients=Union[tableofcoefficientsall];

(*Find the number of coefficients involved.*)
tableofcoefficientslength=
Length[tableofcoefficients];

(*Convert each coefficient to a string.*)
tableofcoefficientsstrings=
Table[ToString[InputForm[tableofcoefficients[[ind]]]],
{ind,1,tableofcoefficientslength}];

(*Set each coefficient to zero and convert to an expression.*)
rulestable1=
Table[ToExpression[tableofcoefficientsstrings[[ind]]<>"==0"],
{ind,1,tableofcoefficientslength}];

(*Create a single table of all of the unknown polynomial
coefficients.*)
unknowns1=Join[xsecoefficienttable,etacoefficienttable];

(*Flatten this table so that internal brackets are removed.*)
unknowns2=Flatten[unknowns1];

unknowns=DeleteCases[unknowns2,0];

(***Be able to look at a table of the Solve equations.***)
solveequations=Table[tableofcoefficientsstrings[[ind]]<>"==0",
{ind,1,tableofcoefficientslength}];

(***Be able to look at the whole set of polynomial coefficients.***)
allunknowns=unknowns;

(***Print out the number of unknowns being solved for.***)
Print[StringForm[
"The number of unknown polynomial coefficients = ``",Length[allunknowns]]];

(***Print out the number of equations being solved.***)
Print[StringForm[
"The number of equations for the polynomial coefficients = ``",Length[rulestable1]]];

(*The expansion polynomials have been substituted into
the determining equations and the coefficients isolated
to form a set of algebraic equations for the polynomial
coefficients. Not every polynomial coefficient appears 
in this system of equations and so the Solve function will
usually return an error message that "Equations may not give
solutions for all "solve" variables." The Off function turns
this error message. There is no error here, the missing
coefficients simply remain free parameters in the infinitesimal.*)
Off[Solve::svars];

(*Solve for the constants in the 'xse' and 'eta polynomial 
expansions.*)
rulestable=Solve[rulestable1,unknowns];

(*Be able to look at the table of rules that has been produced.*)
coefficientrules=rulestable;

(*Create a table of the nonzero polynomial coefficients.*)
allunknowns1a=unknowns/.rulestable;

(*Get rid of internal brackets.*)
allunknowns1=Flatten[allunknowns1a];

(*To fully isolate the nonzero group parameters lets do a 
comparison with the full table of polynomial coefficients.
Generate an index limit equal to the total number of 
polynomial coefficients.*)
lengthunknownstable=Length[unknowns];

(*Create a table with ones at the positions of the 
nonzero polynomial coefficients and zeros elswhere.*)
groupparametersonesandzeros=Table[If[
      MemberQ[allunknowns1,
        unknowns[[ind]],Infinity],1,0],
    {ind,1,lengthunknownstable}];

(*Multiply the table of unknowns by the table of ones and zeros.*)
groupparameterswithzeros=unknowns*groupparametersonesandzeros;

(*Remove all zeros to produce the table of nonzero group parameters.*)
groupparameters1=DeleteCases[groupparameterswithzeros,0];

(*Be able to look at the table of group parameters.*)
groupparameters=groupparameters1;
    
(*Substitute these constants into the table of 'xse' infinitesimal 
functions. Here is the final result.*)
xseinfinitesimals1=Flatten[tableofxxsefunctions//.rulestable];

(*Insert "xsej==" in front of each of the xse expressions.The 
final output is presented as a table of strings.*)
xseinfinitesimalstrings=Table[
"xse"<>ToString[j]<>
"["<>zvariablesundefinedstring<>"]"<>"="<>
ToString[InputForm[xseinfinitesimals1[[j]]]],{j,1,numberindependentvars}];

(*Substitute these constants into the table of 'eta' infinitesimal 
functions. here is the final result.*)
etainfinitesimals1=Flatten[tableofeetafunctions//.rulestable];

(*Insert "etai==" in front of each of the xse expressions.The 
final output is presented as a table of strings.*)
etainfinitesimalstrings=Table[
"eta"<>ToString[i]<>
"["<>zvariablesundefinedstring<>"]"<>"="<>
ToString[InputForm[etainfinitesimals1[[i]]]],{i,1,numberdependentvars}];

(*Use a simpler name for the xse output table of strings.
This is the main output of the SolveDeterminingEquations function.*)
xsefunctions=xseinfinitesimalstrings;

(*Use a simpler name for the eta output table of strings.
This is the main output of the SolveDeterminingEquations function.*)
etafunctions=etainfinitesimalstrings;

(*The number of group parameters.*)
parameterlength=Length[groupparameters1];

(*Finally, here is a table of the infinitesimal groups.*)
infinitesimalgroups1=Table[
      Coefficient[{ToExpression[xsefunctions],
          ToExpression[etafunctions]},{groupparameters1[[ind]],
          groupparameters1[[ind]]}],{ind,1,parameterlength}];

(*Be able to look at the table of the infinitesimal groups.*)
infinitesimalgroups=infinitesimalgroups1;

(*Tell the user the output is ready.*)
Print["SolveDeterminingEquations has finished executing. You can look at the output in the tables xsefunctions and etafunctions. Each entry in these tables is an infinitesimal function in string format expressed in terms of z-variables and the group parameters. The output can also be viewed with the group parameters stripped away by looking at the table infinitesimalgroups. In either case you may wish to convert the z-variables to conventional variables using the table ztableofrules. 

Keep in mind that this function only finds solutions of the determining equations that are of algebraic form. The determining equations may admit solutions that involve transcendental functions and/or integrals. Note that arbitrary functions may appear in the infinitesimals and that these can be detected by running the package function SolveDeterminingEquations for several polynomial orders. If terms of ever increasing order appear, then an arbitrary function is indicated."];])
(*************END*OF*FUNCTION*SolveDeterminingEquations******)

(**********BEGINNING*OF*FUNCTION*MakeCommutatorTable***********)
(*This function is used to generate the commutator table of the Lie 
algebra of the groups contained in 'infinitesimalgroupsxy'.
The symbols used for 'independentvariables' and
'dependentvariables' must correspond to the symbols in
'infinitesimalgroupsxy'. The output is contained in the table
'commutatortable'. This is the last function in the package.*)
MakeCommutatorTable[independentvariablesstrings_, dependentvariablesstrings_, 
   infinitesimalgroupsxy_] := 
 (*The function 'Module' protects variable names from the outside world.*)  
   (Module[{nnumberindependentvars, mnumberdependentvars, numberofgroups, i, j, i1, j1, a, b,
   independentvariables,dependentvariables}, 
   
(*Convert the list of independentvariables to expressions.*)
independentvariables=ToExpression[independentvariablesstrings];

(*Convert the list of dependentvariables to expressions.*)
dependentvariables=ToExpression[dependentvariablesstrings];
    
(*The number of independent variables.*)
nnumberindependentvars = Length[independentvariables]; 
    
(*The number of dependent variables.*)
mnumberdependentvars = Length[dependentvariables]; 
    
(*The number of groups in the list infinitesimalgroupsxy.*)
numberofgroups = 
  Length[infinitesimalgroupsxy]; 
     
(*Construct the commutator table.*)
commutatortable = 
  Table[{Table[Sum[infinitesimalgroupsxy[[a,1,j1]]*
           D[infinitesimalgroupsxy[[b,1,j]], independentvariables[[j1]]], 
          {j1, 1, nnumberindependentvars}] + Sum[infinitesimalgroupsxy[[a,2,i1]]*
           D[infinitesimalgroupsxy[[b,1,j]], dependentvariables[[i1]]], 
          {i1, 1, mnumberdependentvars}] - Sum[infinitesimalgroupsxy[[b,1,j1]]*
           D[infinitesimalgroupsxy[[a,1,j]], independentvariables[[j1]]], 
          {j1, 1, nnumberindependentvars}] - Sum[infinitesimalgroupsxy[[b,2,i1]]*
           D[infinitesimalgroupsxy[[a,1,j]], dependentvariables[[i1]]], 
          {i1, 1, mnumberdependentvars}], {j, 1, nnumberindependentvars}], 
       Table[Sum[infinitesimalgroupsxy[[a,1,j1]]*D[infinitesimalgroupsxy[[b,2,
             i]], independentvariables[[j1]]], {j1, 1, nnumberindependentvars}] + 
         Sum[infinitesimalgroupsxy[[a,2,i1]]*D[infinitesimalgroupsxy[[b,2,
             i]], dependentvariables[[i1]]], {i1, 1, mnumberdependentvars}] - 
         Sum[infinitesimalgroupsxy[[b,1,j1]]*D[infinitesimalgroupsxy[[a,2,
             i]], independentvariables[[j1]]], {j1, 1, nnumberindependentvars}] - 
         Sum[infinitesimalgroupsxy[[b,2,i1]]*D[infinitesimalgroupsxy[[a,2,
             i]], dependentvariables[[i1]]], {i1, 1, mnumberdependentvars}], {i, 1, mnumberdependentvars}]}, 
      {a, 1, numberofgroups}, {b, 1, numberofgroups}];
      
      (*Tell the user the output is ready.*)
Print["MakeCommutatorTable has finished executing. You can look at the output in the table commutatortable. To present the output in the most readable form you may want view it as a matrix using MatrixForm[commutatortable]. Occasionally the entries in the commutatortable will have terms that cancel. To get rid of these terms use the function Simplify before viewing the table."];])
 (*****END*OF*FUNCTION*MakeCommutatorTable********************)





EndPackage[]