(* :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 was Version 2.0 of the package and was included on a CD with my book "Introduction to Symmetry Analysis" published in 2002 as part of the series "Cambridge Texts in Applied Mathematics". Updated versions of the package are available on my website at www.stanford.edu/~cantwell/.

Over the years the package has worked quite well except 
for cases involving Lie-Backlund transformations which, when forming the determining equations, often ran very slowly even for a seemingly modest calculation. The results were correct but the time needed to form the determining equations seemed overly long. I finally got around to addresssing this problem in December 2006. The culprit is in a potentially large number of calls to the Mathematica function "Position" near the end of the call to my function FindDeterminingEquations. At this point the program is sorting through an array containing the left-hand-sides of the determining equations plus a number 
of duplicates. The function "Position" is used to find the 
indices of the duplicates and in a couple of steps remove the duplicates. The problem is that the function "Position" is called a number of times equal to the number of items in the table containing the determining equations plus duplicates. This can be a very large number equal to a thousand or more. Now the duplicates are removed in a single step using the Mathematica set function "Union" one line ahead of the step that calls"Position" and that step is now comented out with beginning (* and ending*) along with a number of following steps that manipulated the list of indices that resulted from the multiple calls to "Position". FindDeterminingEquations now runs much faster. The changes (colored blue) were made final on December 22, 2006.

Another issue that I encountered in the process of making this correction has to do with the application of extended rules in the case of Lie-Backlund transformations. To begin with it is important to be aware that the package is designed to insert these rules into the array InvarianceConditionNoRules without inserting them into the infinitesimals themselves (ie. into the positions where the independent variables appear in the infinitesimals). 
So there is no possibility of an infinitely recursive loop. Therefore there is no need to use a differentiated form of the input equation when searching for Lie-Backlund groups. However it is important to set up the rules so that the highest derivative in the equation is on the left hand side. Otherwise when the rule is substituted into the invariance condition one of the independent derivatives might be replaced by a dependent derivative typically leading to too many determining equations and fewer groups. For example for the heat equation with r=1 the required rule is rulesarray={"D[u[x,t],x,x]->D[u[x,t],t]"} not rulesarray={"D[u[x,t],t]->D[u[x,t],x,x]"}. If the latter is substituted into the invariance condition then places where D[u[x,t],t](an independent variable) appears will be replaced by D[u[x,t],x,x] (a dependent variable). As a result the determining equations will be too restricive, there will be too many of them and fewer symmetries will be found.

In addition I have included a number of diagnostics that can be printed during the run to show the current status of the function FindDeterminingEquations. The printed line usually gives a byte count and time so the user can monitor the run if desired. At the moment most of these are commented out. The user may wish to un-comment these lines so that the status of FindDeterminingEquations can be tracked in some detail. This is particularly useful if the equations being analyzed are long and complicated or if one is searching for Lie-Backlund equations of order greater than one where the invariance condition can be quite long and the number of variables that the infinitesimals depend on is large. These diagnostic lines are in red so the user can easily find them in the text below and un-comment them if so desired.  

This is version 3.0 of the package dated December 24, 2006.

On December 9, 2009 this comment was added here and also in the call to InvarianceConditionRulesApplied below. On December 5, 2009 I got an e-mail pointing out that the package fails to find the symmetries of the equation
Ut+k*x^2*Uxx*Ux^2/(Ux-c*x*Uxx)^2=0. After examining this case it quickly became clear that package is not finding the proper set of determining equations when the rule(s) applied to the invariance condition contain derivatives in the denominator. When the rule Ut->-k*x^2*Uxx*Ux^2/(Ux-c*x*Uxx)^2 is inserted into the invariance condition there will of course be quite a few terms that will be in the form of fractions with various powers of (Ux-c*x*Uxx)^2 in the denominator. In order to generate the proper invariance condition is is necessary to put the entire invariance condition with rules applied over a common denominator and then remove the denominator. This is now done as follows.
Two lines after the Call to the function InvarianceConditionRulesApplied I do the following.
step 1 - note that the variable termsoftheInvarianceConditionRulesApplied is a list containing all the various terms of the invariance condition with rules applied and if the rules involve any denominators some or all of these terms will have denominators that have to be removed
step 2 - Use the function Apply[Plus,termsoftheInvarianceConditionRulesApplied1] to form the sum of these terms
step 3 - Use the function Together[Apply[Plus,termsoftheInvarianceConditionRulesApplied1]] to put the whole sum over a common denominator
step 4 - Use the function Numerator[Together[Apply[Plus,termsoftheInvarianceConditionRulesApplied1]]] to remove the denominator. The result has a 2 out in front of everything and this needs to be brought inside.
step 5 - Use the function Expand[Numerator[Together[Apply[Plus,termsoftheInvarianceConditionRulesApplied1]]]] to bring the 2 inside. This gets rid of the braces around the invariance condition.
step 6 - Use the function Apply[List,Expand[Numerator[Together[Apply[Plus,termsoftheInvarianceConditionRulesApplied1]]]]] to convert the sum to a list of individual terms. Call this list termsoftheInvarianceConditionRulesApplied1nodenominator. In reality all these steps are done using the single command given in step 6. 

Finally, it should be noted that if the original invariance condition with rules applied, termsoftheInvarianceConditionRulesApplied1, has no denominators it is already correct and this process will correctly leave the invariance condition unchanged.

On December 29, 2009 I made several changes to speed up the function SolveDeterminingEquations and reduce the amount of memory used. This function searches for a multivariate polynomial solution to the determining equations. The variables on which the infinitesimals are assumed to depend include all independent and dependent variables and, in the case of Lie-Backlund transformations, all derivatives of the dependent variables up to order r specified by the user. The order of the trial polynomial (the variable "order" in the package) is selected by the user. The current procedure for creating the list of powers of zvariables that appear in the polynomials does a search over all powers up to z1^order*z2^order*z3^order*...*zn^order where n is the number of variables in the list of zvariables, ie, a search that includes powers up to n*order. But only powers up to order are actually needed and items in the list that involve powers exceeding order are set to Null and then deleted. For calculations of point symmetries where the number of variables is not too large this does not present much of a problem but for Lie-Backlund symmetries where the number of variables can easily be considerably greater than ten, and where higher order polynomials are required, this intermediate step unnecessarily uses vastly more memory and time than is required. The new procedure sorts through powers of z variables only up to the required order and deletes any redundancies. This saves a considerable amount of time and memory for Lie-Backlund calculations. 

With the change described above in magenta this becomes version 4.1 of the package now Dated December 29, 2009.

In May 2010 I encountered a rather unusual problem when I was using the package to work out the Maxwell equations with currents and charges. When I used the dependent variable names j1, j2, j3 to denote the currents the package at first seemed to run OK but with some odd replacements. On the first run the array ztableofrules had the string j1 replaced by the string 1. When run a second time the equation in generic variables had 1 where it should have had z8 and ran incorrectly. The problem was that in the function GenerateInvarianceConditionRulesApplied indices j1, j2, j3,…  are used, and j1 ends up with the value 1 which later gets inserted into ztableofrules . When the package is run a second time j1 is still 1 and gets inserted into the input equation when it is expressed in generic variables. To prevent this the command Clear[“j*”] is used at the end of GenerateInvarianceConditionRulesApplied and at the end of FindDeterminingEquations.  

Finally in May 2010 I added the commands Unprotect["`*"]; and ClearAll["`*"]; at the beginning of the package so it can be loaded more than once during a given session.

With the May 2010 changes, the package is designated version 4.2 dated May 20, 2010.*)

(* :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 
             ψ=Ψ[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}={∂y1/∂x1,...,∂ym/∂xn},
second derivatives,      
       y{2}={∂(∂y1/∂x1)/∂x1,...,∂(∂ym/∂xn)/∂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ξj[x,y,y{1},y{2}...,y{r}]
 yi'=yi+sηi[x,y,y{1},y{2}...,y{r}]
 ∂yi'/∂xj'=∂yi/∂xj+sηij[x,y,y{1}...,y{r+1}]
 ∂(∂yi'/∂xj1')/∂xj2'=∂(∂yi/∂xj1)/∂xj2+sηij1j2[x,y,...]
 ......................................................
 ∂(...∂(∂(pyi'/∂xj1')/∂xj2'...)/∂xjp'=
 ∂(...∂(∂(pyi/∂xj1)/∂xj2...)/∂xjp+
             sη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 ξj and ηi.
The infinitesimals of the extended transformation are
given by the formulae,
        ηij=Dηi/Dxj-(∂yi/∂xk)Dξk/Dxj
     ηij1j2=Dηij1/Dxj-(∂(∂yi/∂xj1)/∂xk)Dξk/Dxj2
ηij1j2...jp=Dηij1...j(p-1)/Dxjp-
        (∂...(∂yi/∂xj1)/∂xj2...∂xj(p-1)∂xjk)Dξk/Dxjp      
where D( )/Dxj is the operator of total differentiation 
with respect to xj. The invariance condition is
                      XΨ=0
where the p-th order extended operator X is given by
  X=ξj∂/∂xj+ηi∂/∂yi+ηij1∂/∂yij1+...+ηij1j2..jp∂/∂yij1...jp.
The invariance condition, XΨ=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 ξ's can be set to zero without loss of generality. For
point groups the ξ'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.

  ξ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+...

  η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 ξ 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.";

(*Note that this is not an initialization cell.*)
Begin["`Private`"];

(*These commands were added May 20, 2010 to permit the package to be read in twice if needed.*)
Unprotect["`*"]; 
ClearAll["`*"];

(*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,j1,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];

(*(***Print out the status.***)
Print[StringForm[
"Start the function InvarianceConditionNoRules.
First 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. The memory in use = ``, The time used = ``",MemoryInUse[],TimeUsed[]]];*)

(*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];

(*(***Print out the status.***)
Print[StringForm[
"Call the function GenerateInfinitesimalTable to produce
a table of strings called allpinfinitesimaltableofstrings 
corresponding to the infinitesimals and extensions. The memory in use = ``, The time used = ``",MemoryInUse[],TimeUsed[]]];*)

(*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}];

(*(***Print out the status.***)
Print[StringForm[
"Convert the required infinitesimal strings to expressions. The memory in use = ``, The time used = ``",MemoryInUse[],TimeUsed[]]];*)

(*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;

(*(***Print out the status.***)
Print[StringForm[
"The table finalinvarianceconditiontable contains the invariance condition
in the form of a table where each term is a fully worked 
out term in the invariance condition. The memory in use = ``, The time used = ``",MemoryInUse[],TimeUsed[]]];*)

(*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;

rulesarrayexpandedtable1sep98check=rulesarrayexpandedtable1sep98;

(*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];

(*(***Print out the status.***)
Print[StringForm[
"The call to InvarianceConditionNoRules to produce a table of the terms that appear in the raw invariance condition is complete, the length of termsoftheInvarianceConditionNoRules is = ``, the byte count of termsoftheInvarianceConditionNoRules is = ``, the memory in use = ``, the time used is = ``",Length[termsoftheInvarianceConditionNoRules],ByteCount[termsoftheInvarianceConditionNoRules],MemoryInUse[],TimeUsed[]]];*)

(*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];

(*(***Print out the status.***)
Print[StringForm[
"Call MakeRulesArray to produce an extended table of rules to be applied to the invariance condition for cases with r>0, the memory in use = ``, the time used is = ``",MemoryInUse[],TimeUsed[]]];*)

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;

rulesarrayexpandedcheck=rulesarrayexpanded;

(*(***Print out the status.***)
Print[StringForm[
"The generation of an extended array of rules is complete. 
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. The memory in use = ``, The time used is = ``",MemoryInUse[],TimeUsed[]]];*)

(*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}];

(*(***Print out the status.***)
Print[StringForm[
"Here is the invariance condition with rules applied sorted into a table of terms 
called termsoftheInvarianceConditionRulesApplied1.The length of termsoftheInvarianceConditionRulesApplied1 = ``, The byte count of termsoftheInvarianceConditionRulesApplied1 = ``, the memory in use = ``, the time used = ``",Length[termsoftheInvarianceConditionRulesApplied1],ByteCount[termsoftheInvarianceConditionRulesApplied1],MemoryInUse[],TimeUsed[]]];*)

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

(*This command was added on May 20, 2010. I ran into an odd situation looking for the symmetries of Maxwell's equations with currents and charges where the choice of j1 for a user variable name caused problems. The problem is that by the end of this function j1 has been assigned the value 1. This caused the j1 to be replaced by 1 in several places in particular the list ztableofrules. Clearing all j's at this point prevents the problem.*)
Clear["j*"];])
(*****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,infinitesimaltablesumsunion,infinitesimaltablesumsmatchtbl,
infinitesimaltablesumsmatchtblfirstindices,
finalinfinitesimalsumindices,
determiningequationstable1,arethetwosetsofdetermningequationsthesame,
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},

(***Print out the status.***)
Print[StringForm[
"The function FindDetermining Equations has begun, the memory in use = ``, the time used = ``",MemoryInUse[],TimeUsed[]]];

(*Print[StringForm["The list of independent variables = ``",independentvariablesstrings]];*)

(*Print[StringForm["The list of dependent variables = ``",dependentvariablesstrings]];*)

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

(*Print[StringForm["The list of independent variables = ``",independentvariables]];*)

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

(*Print[StringForm["The list of dependent variables = ``",dependentvariables]];*)

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

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

(*Print[StringForm["The list of frozen variables = ``",frozenstrings]];*)

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

(*Print[StringForm["The number of frozen strings = ``",frozenstringslength]];*)

(*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}];

(*Print[StringForm["The rules for converting input independent variables = ``",inputequationindependentvariablerules]];*)

(*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}];

(*Print[StringForm["The rules for converting input dependent variables = ``",inputequationdependentvariablerules]];*)

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

(*Print[StringForm["The joined rules for converting input variables = ``",inputequationvariablerules]];*)

(*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];


(*Print[StringForm["The table of incorrect parameter names = ``",frozenstringsnewsymbols]];*)

(*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}];

(*Print[StringForm["The table of correct parameter names = ``",frozenstringrules]];*)

(*Print[StringForm["Input equation in user variable names, inputequationuservariables = ``",inputequationuservariables]];*)

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

(*Print[StringForm["Variable names converted to generic names, inputequationbeforefrozenstrings = ``",inputequationbeforefrozenstrings]];*)

(*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];

(*Print[StringForm["Variable names converted to generic names, final input equation in generic variables = ``",inputequation]];*)

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

(*Print[StringForm["The input equation in generic variables = ``",equationgenericvariables]];*)

(*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<>"]"];

(*Print[StringForm["firstj1 = ``",j1]];*)

(*(***Print out the status.***)
Print[StringForm[
"Now call the the function InvarianceConditionRulesApplied, the memory in use = ``, the time used = ``",MemoryInUse[],TimeUsed[]]];*)

(*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];

(*Print[StringForm["firstAj1 = ``",j1]];*)

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

termsoftheInvarianceConditionRulesApplied1check=termsoftheInvarianceConditionRulesApplied1;

(*This comment was added Dec 9, 2009.*)
(**On December 5, 2009 I got an e-mail pointing out that the package fails to find the symmetries of the equation
Ut+k*x^2*Uxx*Ux^2/(Ux-c*x*Uxx)^2=0. After examining this case it quickly became clear that package is not finding the proper set of determining equations when the rule(s) applied to the invariance condition contain derivatives in the denominator. When the rule Ut->-k*x^2*Uxx*Ux^2/(Ux-c*x*Uxx)^2 is inserted into the invariance condition there will of course be quite a few terms that will be in the form of fractions with various powers of (Ux-c*x*Uxx)^2 in the denominator. In order to generate the proper invariance condition is is necessary to put the entire invariance condition with rules applied over a common denominator and then remove the denominator. This is now done as follows.
Two lines after the Call to the function InvarianceConditionRulesApplied I do the following.
step 1 - note that the variable termsoftheInvarianceConditionRulesApplied is a list containing all the various terms of the invariance condition with rules applied and if the rules involve any denominators some or all of these terms will have denominators that have to be removed
step 2 - Use Apply[Plus,termsoftheInvarianceConditionRulesApplied1] to form the sum of these terms
step 3 - Use Together[Apply[Plus,termsoftheInvarianceConditionRulesApplied1]] to put the whole sum over a common denominator
step 4 - Use Numerator[Together[Apply[Plus,termsoftheInvarianceConditionRulesApplied1]]] to remove the denominator. The result has a 2 out in front of everything and this needs to be brought inside.
step 5 - Use Expand[Numerator[Together[Apply[Plus,termsoftheInvarianceConditionRulesApplied1]]]] to bring the 2 inside. This gets rid if the braces around the invariance condition.
step 6 - Use Apply[List,Expand[Numerator[Together[Apply[Plus,termsoftheInvarianceConditionRulesApplied1]]]]] to convert the sum to a list of individual terms. Call this list termsoftheInvarianceConditionRulesApplied1nodenominator. Allthis is done using the single command given in step 6. Finally if the original invariance condition with rules applied, termsoftheInvarianceConditionRulesApplied1, has no denominators it is already correct and this process will leave the invariance conditioned unchanged.**)

(*This line was added Dec 9, 2009.*)
termsoftheInvarianceConditionRulesApplied1nodenominator = Apply[List, Expand[Numerator[Together[Apply[Plus, termsoftheInvarianceConditionRulesApplied1]]]]];

(*Number of terms in this table. This line was modified Dec 9, 2009.*)
termsoftheInvarianceConditionRulesApplied1length=Length[
termsoftheInvarianceConditionRulesApplied1nodenominator];

(*(***Print out the status. This line was modified Dec 9, 2009.***)
Print[StringForm[
"The call to InvarianceConditionRulesApplied to produce a table of the terms that appear in the invariance condition is complete, the length of termsoftheInvarianceConditionRulesApplied1nodenominator = ``, the bytecount of termsoftheInvarianceConditionRulesApplied1nodenominator = ``, the memory in use = ``, the time used = ``",Length[termsoftheInvarianceConditionRulesApplied1nodenominator],ByteCount[termsoftheInvarianceConditionRulesApplied1nodenominator],MemoryInUse[],TimeUsed[]]];*)

(*Split each term into a table of factors. This line was modified Dec 9, 2009.*)
termsoftheInvarianceConditionRulesApplied1splittable=
Table[Cases[
termsoftheInvarianceConditionRulesApplied1nodenominator[[indx]],x_],
{indx,1,
termsoftheInvarianceConditionRulesApplied1length}];

termsoftheInvarianceConditionRulesApplied1splittablecheck=termsoftheInvarianceConditionRulesApplied1splittable;

(*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;


(*(***Print out the status.***)
Print[StringForm[
"The call to GenerateVariableTable to regenerate a table of the variables upon which the infinitesimals are assumed to depend is complete, the length of variablestringtable1 is = ``, the byte count of variablestringtable1 is = ``, the memory in use = ``, the time used = ``",Length[variablestringtable1],ByteCount[variablestringtable1],MemoryInUse[],TimeUsed[]]];*)

(*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;

variablestringtable2check=variablestringtable2;

(*(***Print out the status.***)
Print[StringForm[
"The call to GenerateVariableTable to regenerate a table of y-derivatives that can appear in the problem is complete, the length of variablestringtable2 is = ``, the byte count of variablestringtable2 is = ``, the memory in use = ``, the time used = ``",Length[variablestringtable2],ByteCount[variablestringtable2],MemoryInUse[],TimeUsed[]]];*)

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

(*(***Print out the status.***)
Print[StringForm[
"Check 1, the length of allvariables = ``, the byte count of allvariables = ``, the memory in use = ``, the time used = ``",Length[allvariables],ByteCount[allvariables],MemoryInUse[],TimeUsed[]]];*)

(*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];

freeyderivativevariablescheck=freeyderivativevariables;

(*(***Print out the status.***)
Print[StringForm[
"Check 2, the length of freeyderivativevariables is = ``, the bytecount of freeyderivativevariables is = ``, the memory in use = ``, the time used = ``",Length[freeyderivativevariables],ByteCount[freeyderivativevariables],MemoryInUse[],TimeUsed[]]];*)

(*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}];

(*(***Print out the status.***)
Print[StringForm[
"Check 3, the length of freeyderivativevariablesstrings is = ``, the byte count of freeyderivativevariablesstrings is = ``, the memory in use = ``, the time used = ``",Length[freeyderivativevariablesstrings],ByteCount[freeyderivativevariablesstrings],MemoryInUse[],TimeUsed[]]];*)

(*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}];

lengthterms1tablecheck=lengthterms1table;

(*(***Print out the status.***)
Print[StringForm[
"Check 4, the length of lengthterms1table is = ``, the byte count of lengthterms1table is = ``, the memory in use = ``, the time used = ``",Length[lengthterms1table],ByteCount[lengthterms1table],MemoryInUse[],TimeUsed[]]];*)

(*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]]}];

(*(***Print out the status.***)
Print[StringForm[
"Check 5, the length of termsoftheInvarianceConditionRulesApplied1splittablestrings is = ``, the byte count of termsoftheInvarianceConditionRulesApplied1splittablestrings is = ``, the memory in use = ``, the time used = ``",Length[termsoftheInvarianceConditionRulesApplied1splittablestrings],ByteCount[termsoftheInvarianceConditionRulesApplied1splittablestrings],MemoryInUse[],TimeUsed[]]];*)

(*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}];

truthtablefreevariablescheck=truthtablefreevariables;

(*(***Print out the status.***)
Print[StringForm[
"Check 6, the length of truthtablefreevariables is = ``, the length of truthtablefreevariables is = ``, the memory in use = ``, the time used = ``",Length[truthtablefreevariables],ByteCount[truthtablefreevariables],MemoryInUse[],TimeUsed[]]];*)

(*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}];

truthtablenumberscheck=truthtablenumbers;

(*(***Print out the status.***)
Print[StringForm[
"Check 7, the length of truthtablenumbers is = ``, the byte count of truthtablenumbers is = ``, the memory in use = ``, the time used = ``",Length[truthtablenumbers],ByteCount[truthtablenumbers],MemoryInUse[],TimeUsed[]]];*)

(*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]]}];

freevariablelocationscheck=freevariablelocations;

(*(***Print out the status.***)
Print[StringForm[
"Check 8, the length of freevariablelocations is = ``, the byte count of freevariablelocations is = ``, the memory in use = ``, the time used = ``",Length[freevariablelocations],ByteCount[freevariablelocations],MemoryInUse[],TimeUsed[]]];*)

(*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]]}];

freeyderivativetablecheck=freeyderivativetable;

(*(***Print out the status.***)
Print[StringForm[
"Check 9, the length of freeyderivativetable is = ``, the byte count of freeyderivativetable is = ``, the memory in use = ``, the time used = ``",Length[freeyderivativetable],ByteCount[freeyderivativetable],MemoryInUse[],TimeUsed[]]];*)

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

infinitesimaltermtable1check=infinitesimaltermtable1;

(*(***Print out the status.***)
Print[StringForm[
"Check 10, the length of infinitesimaltermtable1 is = ``, the byte count of infinitesimaltermtable1 is = ``, the memory in use = ``, the time used = ``",Length[infinitesimaltermtable1],ByteCount[infinitesimaltermtable1],MemoryInUse[],TimeUsed[]]];*)

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

ytermstablecheck=ytermstable;

(*(***Print out the status.***)
Print[StringForm[
"Check 11, the length of ytermstable is = ``, the byte count of ytermstable is = ``, the memory in use = ``, the time used = ``",Length[ytermstable],ByteCount[ytermstable],MemoryInUse[],TimeUsed[]]];*)

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

infinitesimalproducttablecheck=infinitesimalproducttable;

(*(***Print out the status.***)
Print[StringForm[
"Check 12, the length of infinitesimalproducttable is = ``, the byte count of infinitesimalproducttable is = ``, the memory in use = ``, the time used = ``",Length[infinitesimalproducttable],ByteCount[infinitesimalproducttable],MemoryInUse[],TimeUsed[]]];*)

(*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;

infinitesimaltablecheck=infinitesimaltable;

(*(***Print out the status.***)
Print[StringForm[
"Check 13, the length of infinitesimaltable is = ``, the byte count of infinitesimaltable is = ``, the memory in use = ``, the time used = ``",Length[infinitesimaltable],ByteCount[infinitesimaltable],MemoryInUse[],TimeUsed[]]];*)

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

ytermstablngthcheck=ytermstablngth;

(*(***Print out the status.***)
Print[StringForm[
"Check 14, the length of ytermstablngth is = ``, the byte count of ytermstablngth is = ``, the memory in use = ``, the time used = ``",Length[ytermstablngth],ByteCount[ytermstablngth],MemoryInUse[],TimeUsed[]]];*)

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

ytermsmatchtablecheck=ytermsmatchtable;

(*(***Print out the status.***)
Print[StringForm[
"Check 15, the length of ytermsmatchtable is = ``, the byte count of ytermsmatchtable is = ``, the memory in use = ``, the time used = ``",Length[ytermsmatchtable],ByteCount[ytermsmatchtable],MemoryInUse[],TimeUsed[]]];*)

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

ytermsmatchtableoflengthscheck=ytermsmatchtableoflengths;

(*(***Print out the status.***)
Print[StringForm[
"Check 16, the length of ytermsmatchtableoflengths is = ``, the byte count of ytermsmatchtableoflengths is = ``, the memory in use = ``, the time used = ``",Length[ytermsmatchtableoflengths],ByteCount[ytermsmatchtableoflengths],MemoryInUse[],TimeUsed[]]];*)

Print["The function FindDeterminingEquations is nearly complete. The invariance condition has been created with all rules applied. The final step in the generation of the determining equations is to sum together terms in the table of invariance condition terms (called infinitesimaltable) that are multiplied by the same combination of products of free y derivatives. The result is the table infinitesimaltablesums corresponding to matching y-derivative expressions. If the invariance condition is long as it often is this process could take a long time since it requires sorting through the table infinitesimaltable once for each possible combination of y derivative products. This is the rate limiting step in the function FindDeterminingEquations.Virtually all other steps are quite fast including the generation of the extended derivatives of the infinitesimals."];

(*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}]];

infinitesimaltablesumscheck=infinitesimaltablesums;

(*(***Print out the status.***)
Print[StringForm[
"Check 17, the length of infinitesimaltablesums is = ``, the byte count of infinitesimaltablesums is = ``, the memory in use = ``, the time used = ``",Length[infinitesimaltablesums],ByteCount[infinitesimaltablesums],MemoryInUse[],TimeUsed[]]];*)

(***Change made December 22, 2006. I am going to look at another way of generating the determining equations so I do not have to perform the ytermstablngth calls to the function Position. This potentially large number of calls can slow the generation of the determining equations hugely.***)
infinitesimaltablesumsunion=Union[infinitesimaltablesums];

infinitesimaltablesumsunioncheck=infinitesimaltablesumsunion;

(*(***Print out the status.***)
Print[StringForm[
"The length of infinitesimaltablesums is = ``, the byte count of infinitesimaltablesums is = ``, the memory in use = ``, the time used = ``",Length[infinitesimaltablesumsunioncheck],ByteCount[infinitesimaltablesumsunioncheck],MemoryInUse[],TimeUsed[]]];*)

(*Change (removal of this step) made December 22, 2006.This is the call that slows down the function FindDeterminingEquations so badly. Instead I just make one call to the Mathematica function "Union" one step above. (*Produce a table containing indices of matching 
infinitesimal sums.*)
infinitesimaltablesumsmatchtbl=Table[
Position[
infinitesimaltablesums,infinitesimaltablesums[[ind]],{1}],
{ind,1,ytermstablngth}];*)

(*Change made December 22, 2006. infinitesimaltablesumsmatchtblcheck=infinitesimaltablesumsmatchtbl;*)

(*Change made December 22, 2006. (***Print out the status.***)
Print[StringForm[
"Check 17 looks like a bottleneck. The memory in use = ``, the time used = ``",MemoryInUse[],TimeUsed[]]];*)

(*Change made December 22, 2006. (***Print out the status.***)
Print[StringForm[
"Check 18, the length of infinitesimaltablesumsmatchtbl is = ``, the byte count of infinitesimaltablesumsmatchtbl is = ``, the memory in use = ``, the time used = ``",Length[infinitesimaltablesumsmatchtbl],ByteCount[infinitesimaltablesumsmatchtbl],MemoryInUse[],TimeUsed[]]];*)

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

(*Change made December 22, 2006. infinitesimaltablesumsmatchtblfirstindicescheck=infinitesimaltablesumsmatchtblfirstindices;*)

(*Change made December 22, 2006. (***Print out the status.***)
Print[StringForm[
"Check 19, the length of infinitesimaltablesumsmatchtblfirstindices is = ``, the byte count of infinitesimaltablesumsmatchtblfirstindices is = ``, the memory in use = ``, the time used = ``",Length[infinitesimaltablesumsmatchtblfirstindices],ByteCount[infinitesimaltablesumsmatchtblfirstindices],MemoryInUse[],TimeUsed[]]];*)

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

(*Change made December 22, 2006. finalinfinitesimalsumindicescheck=finalinfinitesimalsumindices;*)

(*Change made December 22, 2006. (***Print out the status.***)
Print[StringForm[
"Check 20, the length of finalinfinitesimalsumindices is = ``, the byte count of finalinfinitesimalsumindices is = ``, the memory in use = ``, the time used = ``",Length[finalinfinitesimalsumindices],ByteCount[finalinfinitesimalsumindices],MemoryInUse[],TimeUsed[]]];*)

(*Change made December 22, 2006. (*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}]];*)

determiningequationstable1=infinitesimaltablesumsunion;

(*(***Print out the status.***)
Print[StringForm[
"Check 21, the length of determiningequationstable1 is = ``, the byte count of determiningequationstable1 is = ``, the memory in use = ``, the time used = ``",Length[determiningequationstable1],ByteCount[determiningequationstable1],MemoryInUse[],TimeUsed[]]];*)

determiningequationstable1check=determiningequationstable1;

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

(*(***Print out the status.***)
Print[StringForm[
"Check 22, the length of determiningequationstable1length is = ``, the byte count of determiningequationstable1length is = ``, the memory in use = ``, the time used = ``",Length[determiningequationstable1length],ByteCount[determiningequationstable1length],MemoryInUse[],TimeUsed[]]];*)

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

determiningequationstable2check=determiningequationstable2;

(*(***Print out the status.***)
Print[StringForm[
"Check 23, the length of determiningequationstable2 is = ``, the byte count of determiningequationstable2 is = ``, the memory in use = ``, the time used = ``",Length[determiningequationstable2],ByteCount[determiningequationstable2],MemoryInUse[],TimeUsed[]]];*)

(*Change made December 22, 2006. (*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}]];*)

yderivfactortablecheck=yderivfactortable;

(*(***Print out the status.***)
Print[StringForm[
"Check 24, the length of yderivfactortable is = ``, the byte count of yderivfactortable is = ``, the memory in use = ``, the time used = ``",Length[yderivfactortable],ByteCount[yderivfactortable],MemoryInUse[],TimeUsed[]]];*)

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

determiningequationstablecheck=determiningequationstable;

(*(***Print out the status.***)
Print[StringForm[
"Check 25, the length of determiningequationstable is = ``, the byte count of determiningequationstable is = ``, the memory in use = ``, the time used = ``",Length[determiningequationstable],ByteCount[determiningequationstable],MemoryInUse[],TimeUsed[]]];*)

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

(*(***Print out the status.***)
Print[StringForm[
"Check 26, the length of determiningequationstablelengthof is = ``, the byte count of determiningequationstablelengthof is = ``, the memory in use = ``, the time used = ``",Length[determiningequationstablelengthof],ByteCount[determiningequationstablelengthof],MemoryInUse[],TimeUsed[]]];*)

(*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}];

determiningequationstableofstringscheck=determiningequationstableofstrings;

(*(***Print out the status.***)
Print[StringForm[
"Check 27, the length of determiningequationstableofstrings is = ``, the byte count of determiningequationstableofstrings is = ``, the memory in use = ``, the time used = ``",Length[determiningequationstableofstrings],ByteCount[determiningequationstableofstrings],MemoryInUse[],TimeUsed[]]];*)

(*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;

(*(***Print out the status.***)
Print[StringForm[
"The determining equations have been prepared, the length of determiningequations = ``, the byte count of determiningequations = ``, the memory in use = ``, the time used = ``",Length[determiningequations],ByteCount[determiningequations],MemoryInUse[],TimeUsed[]]];*)

(*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;

(***Print out the status.***)
Print[StringForm[
"The determining equations have been expressed in terms of z-variables, the length of zdeterminingequations = ``, the byte count of zdeterminingequations = ``, the memory in use = ``, the time used = ``",Length[zdeterminingequations],ByteCount[zdeterminingequations],MemoryInUse[],TimeUsed[]]];

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

(*Print[StringForm["dummyvariabletableofrulesstringsztoxy = ``",dummyvariabletableofrulesstringsztoxy]];*)

(*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}];

(*Print[StringForm["inputequationindependentvariablerulesxtoinput = ``",inputequationindependentvariablerulesxtoinput]];*)

(*Print[StringForm["dependentvariables = ``",dependentvariables]];*)

(*Print[StringForm["dependentvariablesstrings = ``",dependentvariablesstrings]];*)

(*Convert the list of dependent variables to an expression - again! This refreshes the table of expressions dependentvariables. Change added May 20, 2010.*)
dependentvariables=ToExpression[dependentvariablesstrings];

(*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}];

(*Print[StringForm["inputequationdependentvariablerulesytoinput = ``",inputequationdependentvariablerulesytoinput]];*)

(*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];

(***Print out the status.***)
Print[StringForm[
"FindDeterminingEquations is done. The memory in use = ``, the time used = ``",MemoryInUse[],TimeUsed[]]];

(*This command was also (same as at the end of GenerateInvarianceConditionRulesApplied) added on May 20, 2010. I ran into an odd situation, looking for the symmetries of Maxwell's equations with currents and charges, where the choice of j1 for a user variable name caused problems. The problem at this point (at the end of FindDeterminingEquations) is that j1 has been assigned the value 1. If the notebook is run again the j1 will be replaced by 1 in several places in particular in the input equation in generic variables. Clearing all j's at this point prevents the problem.*)
Clear["j*"];

(*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},

(*Print[StringForm["Step1 = ``",TimeUsed[]]];*)

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

(*Print[StringForm["Step2 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step3 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step4 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step5 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step6 = ``",TimeUsed[]]];*)


(*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}];

(*Print[StringForm["Step7 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step8 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step9 = ``",TimeUsed[]]];*)


(*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}];

(*Print[StringForm["Step10 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step11 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step12 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step13 = ``",TimeUsed[]]];*)


(*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];

(*Print[StringForm["Step14 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step15 = ``",TimeUsed[]]];*)


(*Be able to view the output list of variables. Keep this, Dec 29, 2009.*)
variablestringtableforviewing=variablestringtable;

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

(*Print[StringForm["Step16 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step17 = ``",TimeUsed[]]];*)


(*Be able to view the list of zvariables. Keep this, Dec 29, 2009. *)
zvariablesforviewing=zvariables;

(***Comment out the*OLD Construction of powers of z*****)
(***************December 27, 2009***********************)
(*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}];

Print[StringForm["Step18 = ``",TimeUsed[]]];


(*Be able to view the list of kindices. Keep this*)
kindicesforviewing=kindices;

Print[StringForm["Step19 = ``",TimeUsed[]]];


(*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]];

Print[StringForm["Step20 = ``",TimeUsed[]]];


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

Print[StringForm["Step21 = ``",TimeUsed[]]];


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

Print[StringForm["Step22 = ``",TimeUsed[]]];


(*Be able to view the powersproductstring.  Keep this.*)
powersproductstringforviewing=powersproductstring;

(*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}]];

Print[StringForm["Step23 = ``",TimeUsed[]]];


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

Print[StringForm["Step24 = ``",TimeUsed[]]];


(*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'. This was a huge waste of time and memory, Dec 29, 2009.*)
powerofzs1string=
"Table["<>
"If["<>
kindicessumstring<>"<="<>ToString[order]<>","<>
powersproductstring<>"],"<>
indicestablestring<>"]";

Print[StringForm["Step25 = ``",TimeUsed[]]];


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

Print[StringForm["Step26 = ``",TimeUsed[]]];


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

Print[StringForm["Step27 = ``",TimeUsed[]]];


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

Print[StringForm["Step28 = ``",TimeUsed[]]];


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

(*End Comment out of the OLD construction of powers of z*)
(****************December 27, 2009***********************)

(*Print["We are ready to construct the powers of zvariables that appear in the multivariate polynomial expansions."];*)

(*Print[StringForm["Time = ``",TimeUsed[]]];*)

orderforviewing=order;

numberofvariablesforviewing=numberofvariables;

zvariablesforviewing=zvariables;

(*Print[StringForm["Polynomial order=``",orderforviewing]];*)

(*Print[StringForm["List of zvariables=``",zvariablesforviewing]];*)

(*Print[StringForm["Number of zvariables =``",numberofvariablesforviewing]];*)

(*********Insert NEW construction of powers of z*******)
(***************December 27, 2009**********************)

(*This step constructs a string that is in the right form to create and join tables of z variables of order one up to "order". This did not work until I used the names zvariablesforviewing instead of zvariables and the index name alpha instead of j. I am not sure what the problem was but it would not use the string zvariables which is a protected name. Dec 30, 2009.*)
powertablenewdec2709a = Table[StringJoin["Table[" , StringJoin[Table[StringReplace["zvariablesforviewing[[kalpha1]]*", "alpha1" -> ToString[alpha1]], {alpha1, 1, orderindex - 1}], StringReplace["zvariablesforviewing[[kpolyorderindex]]", "polyorderindex" -> ToString[orderindex]]] , Table[StringReplace[",{kalpha2,1,numberofvariablesforviewing}", "alpha2" -> ToString[alpha2]], {alpha2, 1, orderindex}] , "]"], {orderindex, 1, order}];

(*Print[StringForm["powertablenewdec2709a=``",powertablenewdec2709a]];*)

(*Convert the previous string to an expression. The expression gets evaluated and creates a table of all the powers of zvariables that are needed for the multivariate trial polynomials.*)
powertablenewdec2709b = ToExpression[powertablenewdec2709a];

(*Print[StringForm["powertablenewdec2709b=``",powertablenewdec2709b]];*)

(*This inserts a one in the first position of the table of zvariable powers.*)
powertablenewdec2709withone={1,powertablenewdec2709b};

(*Print[StringForm["powertablenewdec2709withone=``",powertablenewdec2709withone]];*)

(*This removes all brackets within the table of zvariable powers and removes redundancies.*)
powertable = Union[Flatten[powertablenewdec2709withone]];

(*Print[StringForm["powertable=``",powertable]];*)

Print[StringForm["The variable powertablelength is the number of terms required for each multivariate polynomial used for the infinitesimals. This number is determined by the choice of polynomial order and the number of zvariables. The time needed to solve the determining equations increases as powertable increases. powertablelength = ``",Length[powertable]]];

(*Print["The terms required for the infinitesimal polynomials have been constructed."];*)

(*Print[StringForm["Time = ``",TimeUsed[]]];*)


(********End of NEW construction of powers of z*********)
(**************December 27, 2009************************)

(*Be able to view the powertable. Keep this.*)
powertableforviewing=powertable;

(*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];

(*Print[StringForm["Step30 = ``",TimeUsed[]]];*)


(*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}];

(*Print[StringForm["Step31 = ``",TimeUsed[]]];*)


(*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}];

(*Print[StringForm["Step32 = ``",TimeUsed[]]];*)


(*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}];

(*Print[StringForm["Step33 = ``",TimeUsed[]]];*)


(*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}]],
"{"->""],"}"->""];

(*Print[StringForm["Step34 = ``",TimeUsed[]]];*)


(*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}];

(*Print[StringForm["Step35 = ``",TimeUsed[]]];*)


(*Be able to view the tableofxxsefunctionsstrings. Keep this. Dec 29, 2009*)
(*tableofxxsefunctionsstringsforviewing=tableofxxsefunctionsstrings;*)

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

(*Print[StringForm["Step36 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step37 = ``",TimeUsed[]]];*)


(*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}];

(*Print[StringForm["Step38 = ``",TimeUsed[]]];*)


(*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}];

(*Print[StringForm["Step39 = ``",TimeUsed[]]];*)


(*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}];

(*Print[StringForm["Step40 = ``",TimeUsed[]]];*)


(*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}];

(*Print[StringForm["Step41 = ``",TimeUsed[]]];*)


(*Be able to view the tableofeetafunctionsstrings. Keep this.  Dec 29, 2009.*)
(*tableofeetafunctionsstringsforviewing=tableofeetafunctionsstrings;*)

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

(*Print[StringForm["Step42 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step43 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step44 = ``",TimeUsed[]]];*)


(*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"];

(*Print[StringForm["Step45 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step46 = ``",TimeUsed[]]];*)


(*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}];

(*Be able to look at zdeterminingequationstablexxseeetaexpressions. Dec 29, 2009*)
(*zdeterminingequationstablexxseeetaexpressionsforviewing=zdeterminingequationstablexxseeetaexpressions;*)

(*Print["Look at expandedzdeterminingequationstableforviewing before Union is applied."];*)

(*Print[StringForm["zdeterminingequationstablexxseeetaexpressionslength = ``",Length[zdeterminingequationstablexxseeetaexpressions]]];*)


(*Print[StringForm["Step47 = ``",TimeUsed[]]];*)


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


(*Be able to look at all the expansions of determining equations before Union is applied. Dec 29, 2009*)
(*expandedzdeterminingequationstableforviewing=expandedzdeterminingequationstable;*)

(*Print["Look at expandedzdeterminingequationstableforviewing before Union is applied."];*)

(*Print[StringForm["expandedzdeterminingequationstablelength = ``",Length[expandedzdeterminingequationstable]]];*)

Print["The polynomial expansions have been substituted into the determining equations. It is now time to collect the coefficients of various powers of zvariables into a table called table of coefficientsall. This step uses the function CoefficientList and is a fairly time consuming procedure."];

(*Print[StringForm["Time = ``",TimeUsed[]]];*)

Print[StringForm["The memory in use = ``, The time = ``",MemoryInUse[],TimeUsed[]]];


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

(*Be able to look at all the coefficients before Union is applied. Dec 29, 2009*)
(*tableofcoefficientsallforviewing=tableofcoefficientsall;*)

(*Print["Look at tableofcoefficientsallforviewing before Union is applied."];*)

(*Print[StringForm["The table of coefficients collected from the determining equations has been constructed. Because many of the powers of zvariables end up with zero as the coefficient this table will contain a large number of zeros. In the case of Lie-Backlund transformations of high order, the number can be in the many millions. So it is useful to check the length of this table. The length of the list tableofcoefficientsall is = ``",Length[tableofcoefficientsall]]];*)

(*Print[StringForm["Time = ``",TimeUsed[]]];*)

(*Print["The next step is to remove all the zeros in this table using the function DeleteCases. This can take a long time."];*)


(*The tableofcoefficientsall includes a huge number of zeros (millions if this is a high 
order Lie-Backlund calculation with r=3 or more). Before manipulating the table it is 
extremely important to remove all the zeros to reduce the size of the calculation. The 
previous version of the software did this with the function Union which was slow and 
left a single zero in the first slot of the table. In this step added December 23, 2009 we 
remove ALL zeros using the function DeleteCases. Union is used a little further down to 
remove any remaining redundancies.*)

tableofcoefficientsnozeros=DeleteCases[tableofcoefficientsall,_0];

(*Be able to look at tableofcoefficientsnozeros after DeleteCases removing all zeros is applied.*)
(*tableofcoefficientsnozeros=tableofcoefficientsnozeros;*)

(*Print["Look at tableofcoefficientsnozeros before Union is applied."];*)

(*Print[StringForm["The number of coefficients with all the zeros removed is = ``",Length[tableofcoefficientsnozeros]]];*)

(*Print[StringForm["Time = ``",TimeUsed[]]];*)

(*Print["The next step is to use the function Union to remove any remaining redundancies in the table of coefficients."];*)

(*The next step is the one replaced by the process just completed as of Dec 29, 2009..*)
(*The list of coefficients typically involves a large
number of items that are simply zero. The function Union
removes redundant members from the list this step takes an incredibly long time and requires a lot of virtual memeory which slows things down immensely.*)
(*tableofcoefficients=Union[tableofcoefficientsall];*)

(*Now with all the zeros removed, the function Union will work on the remaining coefficients to remove any further redundancies. This takes quite a bit of time, Dec 29, 2009.*)
tableofcoefficients=Union[tableofcoefficientsnozeros];

(*Be able to look at tableofcoefficients after Union is applied.*)
(*tableofcoefficientsforviewing=tableofcoefficients;*)

(*Print["Look at tableofcoefficients after Union is applied."];*)

(*Print[StringForm["The length of the table of coefficients with remaining redundancies removed = ``",Length[tableofcoefficients]]];*)

(*Print[StringForm["Time = ``",TimeUsed[]]];*)


(*Dec 19, 2009. I still have to define tableofcoefficient. 
For Burgers eqn with r=3 and order =3 I end up with 94583296 
equations to solve. Turn back on the Union function 
but test order 2 first.  The Union function reduces the
number of equations from 94353296 down to only 1973 - an
incredibly effective reduction.*)
(*tableofcoefficients=tableofcoefficientsall;*)

(*Print[StringForm["Step50 = ``",TimeUsed[]]];*)


(*The number of remaining coefficients after all the reductions.*)
tableofcoefficientslength=Length[tableofcoefficients];

(*Print[StringForm["Step51 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step52 = ``",TimeUsed[]]];*)

(*Print["At this point the coefficients are converted to expressions. This is the set of equations that will be solved for the nonzero polynomial coefficients that will define the symmetries of the original equation(s)."];*)

(*Print[StringForm["Time = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step53 = ``",TimeUsed[]]];*)


(*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];

(*Print[StringForm["Step54 = ``",TimeUsed[]]];*)


unknowns=DeleteCases[unknowns2,0];

(*Print[StringForm["Step55 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step56 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step57 = ``",TimeUsed[]]];*)


allunknownslength=Length[allunknowns];

(*Print[StringForm["Step58 = ``",TimeUsed[]]];*)


allunknownsexpressions=ToExpression[allunknowns];

(*Print[StringForm["Step59 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step60 = ``",TimeUsed[]]];*)


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

(*Print[StringForm["Step61 = ``",TimeUsed[]]];*)


(*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];

(*Print[StringForm["Step62 = ``",TimeUsed[]]];*)

(***December 20, 2009. I looked at using the commands below in brown to reduce the determining equations by picking out those that involved setting one of the polynomial coefficients to zero and using these as rules to apply to the rest of the equations for the unknown coefficients. Seven levels of reduction were tried. In the end this made no improvement over the single step of just using the mathematica command Solve on the whole system. But note that the reduction with the mathematica commands DeleteCases and Union that are used above to get rid of redundant equations are extremely effective at reducing the number of equations to be solved. Even though this idea seemed to produce no improvement I am going to simply comment out the commands rather than remove them. One never knows when the approach might come in handy.***)

(****Comment out the seven reductions below************)
(*solveequations1=solveequations;


solveequations1length=Length[solveequations1];


(************FirstReduction*******************************)
solveequationsnoequalszero1=StringReplace[solveequations1,{"==0"->""}];


solveequationsnoequalszeroexpressions1=ToExpression[solveequationsnoequalszero1];


solveequationsnoequalszeroexpressionsoneterm1=Cases[solveequationsnoequalszeroexpressions1,x_*y_];


solveequationsnoequalszeroexpressionsonetermtruefalse1=Table[MemberQ[solveequationsnoequalszeroexpressionsoneterm1,allunknownsexpressions[[i]],Infinity],{i,1,allunknownslength}];


solveequationsnoequalszeroexpressionsonetermtruefalse1length=Length[solveequationsnoequalszeroexpressionsonetermtruefalse1];


solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero1=solveequationsnoequalszeroexpressionsonetermtruefalse1/.{False->0,True->1};


coefficientsthatarezero1=Table[solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero1[[i]]*allunknownsexpressions[[i]],{i,1,allunknownslength}];


coefficientsthatarezeroonly1=DeleteCases[coefficientsthatarezero1,0];


coefficientsthatarezeroonly1length=Length[coefficientsthatarezeroonly1];


coefficientsthatarezeroonlystrings1=Table[ToString[coefficientsthatarezeroonly1[[i]]]<>"==0",{i,1,coefficientsthatarezeroonly1length}];


coefficientsthatarezeroonlystringsrules1=Table[ToString[coefficientsthatarezeroonly1[[i]]]<>"->0",{i,1,coefficientsthatarezeroonly1length}];


coefficientsthatarezeroonlystringsrulesexpressions1=Table[ToExpression[coefficientsthatarezeroonlystringsrules1[[i]]],{i,1,coefficientsthatarezeroonly1length}];


(******************************************************)
solveequations2expressions=DeleteCases[Table[ToExpression[solveequations1[[i]]],{i,1,solveequations1length}]/.coefficientsthatarezeroonlystringsrulesexpressions1,True];


solveequations2length=Length[solveequations2expressions];


solveequations2=Table[ToString[solveequations2expressions[[i]]],{i,1,solveequations2length}];


(********************************************************)
(*************************SecondReduction****************)
solveequationsnoequalszero2=StringReplace[solveequations2,{"==0"->""}];


solveequationsnoequalszeroexpressions2=ToExpression[solveequationsnoequalszero2];


solveequationsnoequalszeroexpressionsoneterm2=Cases[solveequationsnoequalszeroexpressions2,x_*y_];


solveequationsnoequalszeroexpressionsonetermtruefalse2=Table[MemberQ[solveequationsnoequalszeroexpressionsoneterm2,allunknownsexpressions[[i]],Infinity],{i,1,allunknownslength}];


solveequationsnoequalszeroexpressionsonetermtruefalse2length=Length[solveequationsnoequalszeroexpressionsonetermtruefalse2];


solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero2=solveequationsnoequalszeroexpressionsonetermtruefalse2/.{False->0,True->1};


coefficientsthatarezero2=Table[solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero2[[i]]*allunknownsexpressions[[i]],{i,1,allunknownslength}];


coefficientsthatarezeroonly2=DeleteCases[coefficientsthatarezero2,0];


coefficientsthatarezeroonly2length=Length[coefficientsthatarezeroonly2];


coefficientsthatarezeroonlystrings2=Table[ToString[coefficientsthatarezeroonly2[[i]]]<>"==0",{i,1,coefficientsthatarezeroonly2length}];


coefficientsthatarezeroonlystringsrules2=Table[ToString[coefficientsthatarezeroonly2[[i]]]<>"->0",{i,1,coefficientsthatarezeroonly2length}];


coefficientsthatarezeroonlystringsrulesexpressions2=Table[ToExpression[coefficientsthatarezeroonlystringsrules2[[i]]],{i,1,coefficientsthatarezeroonly2length}];


(********************************************************)
(********************************************************)
solveequations3expressions=DeleteCases[Table[ToExpression[solveequations2[[i]]],{i,1,solveequations2length}]/.coefficientsthatarezeroonlystringsrulesexpressions2,True];


solveequations3length=Length[solveequations3expressions];


solveequations3=Table[ToString[solveequations3expressions[[i]]],{i,1,solveequations3length}];


(*********************************************************)
(*******************ThirdReduction************************)
solveequationsnoequalszero3=StringReplace[solveequations3,{"==0"->""}];


solveequationsnoequalszeroexpressions3=ToExpression[solveequationsnoequalszero3];


solveequationsnoequalszeroexpressionsoneterm3=Cases[solveequationsnoequalszeroexpressions3,x_*y_];


solveequationsnoequalszeroexpressionsonetermtruefalse3=Table[MemberQ[solveequationsnoequalszeroexpressionsoneterm3,allunknownsexpressions[[i]],Infinity],{i,1,allunknownslength}];


solveequationsnoequalszeroexpressionsonetermtruefalse3length=Length[solveequationsnoequalszeroexpressionsonetermtruefalse3];


solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero3=solveequationsnoequalszeroexpressionsonetermtruefalse3/.{False->0,True->1};


coefficientsthatarezero3=Table[solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero3[[i]]*allunknownsexpressions[[i]],{i,1,allunknownslength}];


coefficientsthatarezeroonly3=DeleteCases[coefficientsthatarezero3,0];


coefficientsthatarezeroonly3length=Length[coefficientsthatarezeroonly3];


coefficientsthatarezeroonlystrings3=Table[ToString[coefficientsthatarezeroonly3[[i]]]<>"==0",{i,1,coefficientsthatarezeroonly3length}];


coefficientsthatarezeroonlystringsrules3=Table[ToString[coefficientsthatarezeroonly3[[i]]]<>"->0",{i,1,coefficientsthatarezeroonly3length}];


coefficientsthatarezeroonlystringsrulesexpressions3=Table[ToExpression[coefficientsthatarezeroonlystringsrules3[[i]]],{i,1,coefficientsthatarezeroonly3length}];


(********************************************************)
(********************************************************)
solveequations4expressions=DeleteCases[Table[ToExpression[solveequations3[[i]]],{i,1,solveequations3length}]/.coefficientsthatarezeroonlystringsrulesexpressions3,True];


solveequations4length=Length[solveequations4expressions];


solveequations4=Table[ToString[solveequations4expressions[[i]]],{i,1,solveequations4length}];


(********************************************************)
(********************FourthReduction*********************)
solveequationsnoequalszero4=StringReplace[solveequations4,{"==0"->""}];


solveequationsnoequalszeroexpressions4=ToExpression[solveequationsnoequalszero4];


solveequationsnoequalszeroexpressionsoneterm4=Cases[solveequationsnoequalszeroexpressions4,x_*y_];


solveequationsnoequalszeroexpressionsonetermtruefalse4=Table[MemberQ[solveequationsnoequalszeroexpressionsoneterm4,allunknownsexpressions[[i]],Infinity],{i,1,allunknownslength}];


solveequationsnoequalszeroexpressionsonetermtruefalse4length=Length[solveequationsnoequalszeroexpressionsonetermtruefalse4];


solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero4=solveequationsnoequalszeroexpressionsonetermtruefalse4/.{False->0,True->1};


coefficientsthatarezero4=Table[solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero4[[i]]*allunknownsexpressions[[i]],{i,1,allunknownslength}];


coefficientsthatarezeroonly4=DeleteCases[coefficientsthatarezero4,0];


coefficientsthatarezeroonly4length=Length[coefficientsthatarezeroonly4];


coefficientsthatarezeroonlystrings4=Table[ToString[coefficientsthatarezeroonly4[[i]]]<>"==0",{i,1,coefficientsthatarezeroonly4length}];


coefficientsthatarezeroonlystringsrules4=Table[ToString[coefficientsthatarezeroonly4[[i]]]<>"->0",{i,1,coefficientsthatarezeroonly4length}];


coefficientsthatarezeroonlystringsrulesexpressions4=Table[ToExpression[coefficientsthatarezeroonlystringsrules4[[i]]],{i,1,coefficientsthatarezeroonly4length}];


(********************************************************)
(********************************************************)
solveequations5expressions=DeleteCases[Table[ToExpression[solveequations4[[i]]],{i,1,solveequations4length}]/.coefficientsthatarezeroonlystringsrulesexpressions4,True];


solveequations5length=Length[solveequations5expressions];


solveequations5=Table[ToString[solveequations5expressions[[i]]],{i,1,solveequations5length}];


(*******************************************************)
(********************FifthReduction*********************)
solveequationsnoequalszero5=StringReplace[solveequations5,{"==0"->""}];


solveequationsnoequalszeroexpressions5=ToExpression[solveequationsnoequalszero5];


solveequationsnoequalszeroexpressionsoneterm5=Cases[solveequationsnoequalszeroexpressions5,x_*y_];


solveequationsnoequalszeroexpressionsonetermtruefalse5=Table[MemberQ[solveequationsnoequalszeroexpressionsoneterm5,allunknownsexpressions[[i]],Infinity],{i,1,allunknownslength}];


solveequationsnoequalszeroexpressionsonetermtruefalse5length=Length[solveequationsnoequalszeroexpressionsonetermtruefalse5];


solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero5=solveequationsnoequalszeroexpressionsonetermtruefalse5/.{False->0,True->1};


coefficientsthatarezero5=Table[solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero5[[i]]*allunknownsexpressions[[i]],{i,1,allunknownslength}];


coefficientsthatarezeroonly5=DeleteCases[coefficientsthatarezero5,0];


coefficientsthatarezeroonly5length=Length[coefficientsthatarezeroonly5];


coefficientsthatarezeroonlystrings5=Table[ToString[coefficientsthatarezeroonly5[[i]]]<>"==0",{i,1,coefficientsthatarezeroonly5length}];


coefficientsthatarezeroonlystringsrules5=Table[ToString[coefficientsthatarezeroonly5[[i]]]<>"->0",{i,1,coefficientsthatarezeroonly5length}];


coefficientsthatarezeroonlystringsrulesexpressions5=Table[ToExpression[coefficientsthatarezeroonlystringsrules5[[i]]],{i,1,coefficientsthatarezeroonly5length}];


(********************************************************)
(********************************************************)
solveequations6expressions=DeleteCases[Table[ToExpression[solveequations5[[i]]],{i,1,solveequations5length}]/.coefficientsthatarezeroonlystringsrulesexpressions5,True];


solveequations6length=Length[solveequations6expressions];


solveequations6=Table[ToString[solveequations6expressions[[i]]],{i,1,solveequations6length}];


(********************************************************)
(*************SixthReduction*****************************)
solveequationsnoequalszero6=StringReplace[solveequations6,{"==0"->""}];


solveequationsnoequalszeroexpressions6=ToExpression[solveequationsnoequalszero6];


solveequationsnoequalszeroexpressionsoneterm6=Cases[solveequationsnoequalszeroexpressions6,x_*y_];


solveequationsnoequalszeroexpressionsonetermtruefalse6=Table[MemberQ[solveequationsnoequalszeroexpressionsoneterm6,allunknownsexpressions[[i]],Infinity],{i,1,allunknownslength}];


solveequationsnoequalszeroexpressionsonetermtruefalse6length=Length[solveequationsnoequalszeroexpressionsonetermtruefalse6];


solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero6=solveequationsnoequalszeroexpressionsonetermtruefalse6/.{False->0,True->1};


coefficientsthatarezero6=Table[solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero6[[i]]*allunknownsexpressions[[i]],{i,1,allunknownslength}];


coefficientsthatarezeroonly6=DeleteCases[coefficientsthatarezero6,0];


coefficientsthatarezeroonly6length=Length[coefficientsthatarezeroonly6];


coefficientsthatarezeroonlystrings6=Table[ToString[coefficientsthatarezeroonly6[[i]]]<>"==0",{i,1,coefficientsthatarezeroonly6length}];


coefficientsthatarezeroonlystringsrules6=Table[ToString[coefficientsthatarezeroonly6[[i]]]<>"->0",{i,1,coefficientsthatarezeroonly6length}];


coefficientsthatarezeroonlystringsrulesexpressions6=Table[ToExpression[coefficientsthatarezeroonlystringsrules6[[i]]],{i,1,coefficientsthatarezeroonly6length}];


(********************************************************)
(********************************************************)
solveequations7expressions=DeleteCases[Table[ToExpression[solveequations6[[i]]],{i,1,solveequations6length}]/.coefficientsthatarezeroonlystringsrulesexpressions6,True];


solveequations7length=Length[solveequations7expressions];


solveequations7=Table[ToString[solveequations7expressions[[i]]],{i,1,solveequations7length}];


(*******************************************************)
(****************SeventhReduction***********************)
solveequationsnoequalszero7=StringReplace[solveequations7,{"==0"->""}];


solveequationsnoequalszeroexpressions7=ToExpression[solveequationsnoequalszero7];


solveequationsnoequalszeroexpressionsoneterm7=Cases[solveequationsnoequalszeroexpressions7,x_*y_];


solveequationsnoequalszeroexpressionsonetermtruefalse7=Table[MemberQ[solveequationsnoequalszeroexpressionsoneterm7,allunknownsexpressions[[i]],Infinity],{i,1,allunknownslength}];


solveequationsnoequalszeroexpressionsonetermtruefalse7length=Length[solveequationsnoequalszeroexpressionsonetermtruefalse7];


solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero7=solveequationsnoequalszeroexpressionsonetermtruefalse7/.{False->0,True->1};


coefficientsthatarezero7=Table[solveequationsnoequalszeroexpressionsonetermtruefalseoneorzero7[[i]]*allunknownsexpressions[[i]],{i,1,allunknownslength}];


coefficientsthatarezeroonly7=DeleteCases[coefficientsthatarezero7,0];


coefficientsthatarezeroonly7length=Length[coefficientsthatarezeroonly7];


coefficientsthatarezeroonlystrings7=Table[ToString[coefficientsthatarezeroonly7[[i]]]<>"==0",{i,1,coefficientsthatarezeroonly7length}];


coefficientsthatarezeroonlystringsrules7=Table[ToString[coefficientsthatarezeroonly7[[i]]]<>"->0",{i,1,coefficientsthatarezeroonly7length}];


coefficientsthatarezeroonlystringsrulesexpressions7=Table[ToExpression[coefficientsthatarezeroonlystringsrules7[[i]]],{i,1,coefficientsthatarezeroonly7length}];


(*******************************************************)
(*******************************************************)
solveequations8expressions=DeleteCases[Table[ToExpression[solveequations7[[i]]],{i,1,solveequations7length}]/.coefficientsthatarezeroonlystringsrulesexpressions7,True];


solveequations8length=Length[solveequations8expressions];


solveequations8=Table[ToString[solveequations8expressions[[i]]],{i,1,solveequations8length}];


(*******************************************************)
solutionrules8=Solve[solveequations8expressions,allunknownsexpressions];

allcoefficients=Join[coefficientsthatarezeroonlystringsrulesexpressions1,coefficientsthatarezeroonlystringsrulesexpressions2,coefficientsthatarezeroonlystringsrulesexpressions3,coefficientsthatarezeroonlystringsrulesexpressions4,coefficientsthatarezeroonlystringsrulesexpressions5,coefficientsthatarezeroonlystringsrulesexpressions6,coefficientsthatarezeroonlystringsrulesexpressions7,Flatten[solutionrules8]];

rulestablereducedsevenlevels=allcoefficients;*)

(*Print[StringForm["Step63 = ``",TimeUsed[]]];*)

Print["Now it we are ready to use the function Solve to find the nonzero polynomial coefficients corresponding to the symmetries of the input equation(s). This can take a while."];

(*Print[StringForm["Time = ``",TimeUsed[]]];*)

Print[StringForm["The memory in use = ``, The time = ``",MemoryInUse[],TimeUsed[]]];

(*Solve for the constants in the 'xse' and 'eta polynomial 
expansions. This was commented out and is not needed with the addition above in brown. Dec 29, 2009, the seven level reduction above turned out to make no improvement so I just sticking with the original approach of using the function Solve on the full (reduced) system.*)
rulestable=Solve[rulestable1,unknowns];

Print["Solve has finished."];

(*Print[StringForm["Time = ``",TimeUsed[]]];*)

(*Be able to look at the table of rules that has been produced. Commented out Dec 29, 2009.*)
(*coefficientrules=rulestablereducedsevenlevels;*)

coefficientrules=rulestable;

(*Print[StringForm["Step65 = ``",TimeUsed[]]];*)

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

allunknowns1a=unknowns/.rulestable;

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

(*Print[StringForm["Step66 = ``",TimeUsed[]]];*)

(*To fully isolate the nonzero group parameters let's 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];

(*Print[StringForm["Step67 = ``",TimeUsed[]]];*)

(*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}];
    
(*Print[StringForm["Step68 = ``",TimeUsed[]]];*)    

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

(*Print[StringForm["Step69 = ``",TimeUsed[]]];*)

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

(*Print[StringForm["Step70 = ``",TimeUsed[]]];*)

(*Be able to look at the table of group parameters.*)
groupparameters=groupparameters1;

(*Print[StringForm["Step71 = ``",TimeUsed[]]];*)
    
(*Substitute these constants into the table of 'xse' infinitesimal functions. Here is the final result.*)
(*xseinfinitesimals1=Flatten[tableofxxsefunctions//.rulestablereducedsevenlevels];*)

xseinfinitesimals1=Flatten[tableofxxsefunctions//.rulestable];

(*Print[StringForm["Step72 = ``",TimeUsed[]]];*)

(*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}];

(*Print[StringForm["Step73 = ``",TimeUsed[]]];*)

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

etainfinitesimals1=Flatten[tableofeetafunctions//.rulestable];

(*Print[StringForm["Step74 = ``",TimeUsed[]]];*)

(*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}];

(*Print[StringForm["Step75 = ``",TimeUsed[]]];*)

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

(*Print[StringForm["Step76 = ``",TimeUsed[]]];*)

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

(*Print[StringForm["Step77 = ``",TimeUsed[]]];*)

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

(*Print[StringForm["Step78 = ``",TimeUsed[]]];*)

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

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

Print["The function SolveDeterminingEquations is finished executing."];

(*Print[StringForm["Time = ``",TimeUsed[]]];*)

Print[StringForm["The memory in use = ``, The time = ``",MemoryInUse[],TimeUsed[]]];

(*Tell the user the output is ready.*)
Print["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 polynomial 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********************)

(*Note that this is not an initialization cell.*)
End[]

Protect[GenerateVariableTable,SecondTerm,PthInfinitesimal,GenerateInfinitesimalTable,InvarianceConditionNoRules,InvarianceConditionRulesApplied,MakeRulesArray,FindDeterminingEquations,SolveDeterminingEquations,MakeCommutatorTable];

EndPackage[]