
(* Time-stamp: Monday, July 20, 2009, 15:50 Mountain Daylight Time *)

(*****************************************************************************)
(*                                                                           *)
(*               *** M A T H E M A T I C A   P R O G R A M ***               *)
(*                                                                           *)
(*            HOMOTOPY METHOD FOR INVERSION OF DERIVATIVES IN 1-D            *)
(*                      AND DIVERGENCES IN 2-D AND 3-D                       *)
(*                                                                           *)
(* Program name: homotopyOperator                                            *)
(*                                                                           *)
(* Purpose: To apply the homotopy method so that a 1-dimensional (1-D)       *)
(*          differential expression can be integrated, or so that a          *)
(*          known divergence in 2-D or 3-D can be inverted.                  *)
(*                                                                           *)
(* Input to homotopyOperator:                                                *)
(*          Any differential expression with one, two or three independent   *)
(*          variables. Independent variables can appear explicitly.          *)
(*          A list of independent variables.                                 *)
(*          A list of parameters if any exist in the expression.             *)
(*                                                                           *)
(* Output : Integrated form of a one-dimensional expression, or an inverted  *)
(*          divergence with most (if not all) curl terms removed.            *)
(*                                                                           *)
(* Language : Mathematica 5.0, 5.1, 5.2, 6.0, and 7.0                        *)
(*                                                                           *)
(* Authors: Douglas Poole and Willy Hereman                                  *)
(*          Department of Mathematical and Computer Sciences                 *)
(*          Colorado School of Mines                                         *)
(*          Golden, CO 80401-1887, U.S.A.                                    *)
(*                                                                           *)
(* Version 1.1: 15 May, 2009                                                 *)
(*                                                                           *)
(*****************************************************************************)

(* Clear["Global`*"];                                                        *)

(* ------------------------------------------------------------------------- *)

(*****************************************************************************)
(* commentinter[]: prints a welcome message                                  *)
(*****************************************************************************)

commentinter[] := Block[{},
Print["*********************************************************"];
Print["        WELCOME TO THE MATHEMATICA PROGRAM FOR           "];
Print["   THE HOMOTOPY METHOD FOR INTEGRATION BY PARTS IN 1-D   "];
Print["     AND THE INVERSION OF DIVERGENCES IN 2-D AND 3-D     "];
Print["             by DOUGLAS POOLE and WILLY HEREMAN          "];
Print["            Version 1.1 released on May 15, 2009         "];
Print["                                                         "];
Print["*********************************************************"];
Print[" "];
Print["Loading code homotopy.m of July 20, 2009."]
]; (* end commentinter *)

(**********************************************************************)
(* The list of Options controls print  and debug statements.          *)
(* The option is given as a rule in the opening function of the       *)
(* program.                                                           *)
(**********************************************************************)

(* ##########     Function: variationalDerivativeMultiD      ########## *)

(************************************************************************)
(* variationalDerivativeMultiD[expression, dependentVar, indVarList]    *)
(* Purpose:  To find the variational derivative of the given expression *)
(*           in order to determine if the expression in integrable.     *)
(* Input:    An expression                                              *)
(*           A number for the dependent variable                        *)
(* Output:   A  partial derivative of the expression with respect       *)
(*           to terms over the jet space                                *)
(*           The order of the expression with respect to the            *)
(*           independent variable                                       *)
(* Created:  23 May, 2006 by DP                                         *)
(* Code is in File:  newhomotopyoper/hovarder.m                         *)
(* Last Modified: 10 March, 2008 15:00 by DP at CSM                     *)
(************************************************************************)

variationalDerivativeMultiD[expression_, dependentVar_, indVarList_] :=
Module[{lengthIndVars, derivativeList, orderList, higherOrderList, dummyderiv,
    depVarList, partialWRTpartialList, combinedList, variationalDerivative,
    termsWithDerList, i, j, indVar1 = Part[indVarList, 1], indVar2, indVar3,
    variationalDerivativeList},

  If[debugVariationalDerivative, printVDMD=Print, Clear[printVDMD],
      Clear[printVDMD]];
  printVDMD["debug code: VDMD, Function: variationalDerivativeMultiD,"<>
      " File: hovarder.m"];

  printVDMD["At VDMD IN, the given expression is:"];
  printVDMD[expression];
  printVDMD["At VDMD IN, the current dependent variable is ", dependentVar];
  printVDMD["At VDMD IN, the independent variables are ", indVarList];

  (* Assign a standard form to the independent variables.               *)
  lengthIndVars = Length[indVarList];
  If[lengthIndVars > 1, indVar2 = Part[indVarList, 2]];
  If[lengthIndVars > 2, indVar3 = Part[indVarList, 3]];

  (* Create a list of patial derivatives existing in the expression.    *)
  termsWithDerList = Cases[expression, Derivative[__][u[dependentVar]]
      [Apply[Sequence, indVarList]], {0, Infinity}];
  printVDMD["At VDMD 1, all derivatives in the expression are given in"<>
      " the following list:"];
  printVDMD[termsWithDerList];

  (* Check for terms that have no derivatives.                          *)
  depVarList = Cases[expression, u[dependentVar][Apply[Sequence, indVarList]],
      {0, Infinity}];

  (* Put depVarList and derivative lists together.                      *)
  combinedList = Union[depVarList, termsWithDerList];
  printVDMD["At VDMD 2, derivative terms and non derivative terms"<>
      " contained in the given expression are shown in the following list:"];
  printVDMD[combinedList];

  (* Create a list of the orders of all derivatives that exist in the   *)
  (* expression.                                                        *)
  higherOrderList = Map[Head, Map[Head, termsWithDerList]] /.
      Derivative -> List;

  (* Include order if terms with no derivatives exist.                  *)
  orderList = If[depVarList =!= {}, Union[{Table[0, {i, 1, lengthIndVars}]},
      higherOrderList], Union[higherOrderList], Union[higherOrderList]];
  printVDMD["At VDMD 3, a list of orders for all terms is:"];
  printVDMD[orderList];

  (* Create a list of the partial derivatives of the expression with    *)
  (* respect to the partial derivative list created above.              *)
  partialWRTpartialList = Map[D[expression, #] &, combinedList];
  printVDMD["At VDMD 4, the partial derivatives of the expression with"<>
      " respect to the terms given in the term list above are given in"<>
      " the following list:"];
  printVDMD[partialWRTpartialList];
  derivativeList = Table[{Part[{indVar1, indVar2, indVar3}, i],
      Part[orderList, j, i]},
      {j, 1, Length[orderList]}, {i, 1, lengthIndVars}];

  (* Compute a list of terms that will add together to form the         *)
  (* variational derivative by taking a series of derivatives the       *)
  (* previous list.                                                     *)
  variationalDerivativeList = Map[(-1)^Apply[Plus, #] &, orderList]*
      ReleaseHold[Thread[dummyderiv[partialWRTpartialList,
      Map[Hold[Apply[Sequence, #]] &, derivativeList]]]] /. dummyderiv -> D;
  printVDMD["At VDMD 5, after all derivatives have been computed,"<>
      " the list of terms becomes:"];
  printVDMD[variationalDerivativeList];

  (* Add the terms in the variational derivative list.                  *)
  variationalDerivative = If[Cases[variationalDerivativeList,
      _?(Head[#] === Sin || Head[#] === Cos || Head[#] === Tan ||
      Head[#] === Sinh || Head[#] === Cosh || Head[#] === Tanh &),
      {0, Infinity}] =!= {},
      TrigExpand[Apply[Plus, variationalDerivativeList]],
      Expand[Apply[Plus, variationalDerivativeList]],
      variationalDerivativeList];
  variationalDerivative = If[variationalDerivative =!= 0,
      Factor[variationalDerivative], variationalDerivative,
      variationalDerivative];
  printVDMD["At VDMD 6, the variational derivative (result from"<>
      " applying the zero Euler operator) for u[", dependentVar, "] is:"];
  printVDMD[variationalDerivative];

  Clear[lengthIndVars, derivativeList, orderList, higherOrderList, dummyderiv,
    partialWRTpartialList, combinedList, variationalDerivativeList,
    depVarList, termsWithDerList, i, j, indVar1, indVar2, indVar3];

  Return[variationalDerivative]
]; (* end Module variationalDerivativeMultiD *)

(* ##########         Function: homotopyStripper           ########## *)

(**********************************************************************)
(* homotopyStripper[givenexpr, indVarsList, paramswithweight]         *)
(* Purpose: Cancels the numerical factors of each term of the         *)
(*          argument and puts these into a list after cancellation    *)
(* Input:   An expression                                             *)
(*          A list of independent variables                           *)
(*          A list of weighted parameters is optional                 *)
(* Output:  A list of the terms of the expression with all constant   *)
(*          coefficients removed.                                     *)
(* Adapted From:  data/newzeala/reu2004/reu2004-0811/software/        *)
(*                whlstrip.m and conservationlaws/dpnstrip.m          *)
(* Code is in File:  newhomotopyoper/honstrip.m                       *)
(* Major Changes Made: 22 May, 2006, 14:30 by DP at CSM               *)
(* Last Modified:  10 March, 2007, 12:31 by DP at CSM                 *)
(**********************************************************************)

homotopyStripper[givenexpr_, indVarsList_, paramswithweight_:{}] :=
Module[{expnList, expnFnList, coefficientList, singleNoCoefficientList,
    strippedTermsList, expandExpr, n1, rules1, rules2, rules3, printHST},

  If[debugHomotopyStripper, printHST = Print, Clear[printHST],
      Clear[printHST]];
  printHST["debug code: HST, Function: homotopyStripper, File: honstrip.m"];

  (* If input expression is zero, deal with it separately.              *)
  If[givenexpr === 0, strippedTermsList = {0};
  Return[strippedTermsList]];

  (* Input expression must be in expanded form.                         *)
  expandExpr = Expand[givenexpr];

  (* Create a list of terms in the given expression.                    *)
  expnList = If[Head[expandExpr] === Plus, Apply[List, expandExpr],
      {expandExpr}, {expandExpr}];
  printHST["At HST 1, a list of terms in the given expression is:"];
  printHST[expnList];

  (* Separate the coefficients and parameters from the terms            *)
  expnFnList = Apply[Power, Map[FactorList, expnList], {2, 2}];
  printHST["At HST 2, each term in expnList has been changed into a"<>
      " list of factors:"];
  printHST[expnFnList];

  (* Create a list of coefficients - not including any weighted         *)
  (* parameters. Two rules are developed to remove all unwanted parts.  *)
  rules1 = _[n1_, ___] /; MemberQ[indVarsList, n1] -> {};
  rules2 = Map[Rule[#, {}] &, paramswithweight];
  rules3 = Map[Rule[#, {}] &, indVarsList];
  printHST["At HST 4, the rules for removing coefficients and parameters:"];
  printHST["rules1 = ", rules1, "\nrules2 = ", rules2 ,"\nrules3 = ", rules3];

  coefficientList = expnFnList /. rules1;
  coefficientList = coefficientList /. rules2;
  coefficientList = coefficientList /. rules3;
  printHST["At HST 4, a list of coefficients where all unwanted parts"<>
      " have been changed to {}:"];
  printHST[coefficientList];

  (* The {} is used first to avoid getting zeros in denominators when   *)
  (* rational terms are given in the problem.                           *)
  coefficientList = coefficientList /. {} -> 1;
  printHST["At HST 5, a list of coefficients where {} is changed to 1:"];
  printHST[coefficientList];

  (* Combine the coefficients for each term.                            *)
  singleNoCoefficientList = Apply[Times, coefficientList, {1, 1}];
  printHST["At HST 6, a combined coefficient list is:"];
  printHST[singleNoCoefficientList];

  (* Form a list of terms without coefficients, but including weighted  *)
  (* parameters.                                                        *)
  strippedTermsList = If[Length[expnFnList] > 1,
      Thread[Times[expnList, 1/singleNoCoefficientList]],
      expnList/singleNoCoefficientList];
  strippedTermsList = Factor[strippedTermsList];
  printHST["At HST 7, a list of stripped terms is:"];
  printHST[strippedTermsList];

  Clear[expnList, expnFnList, coefficientList, singleNoCoefficientList,
    expandExpr, n1, rules1, rules2, rules3, printHST];

  Return[strippedTermsList];
]; (* end Module homotopyStripper *)

(* ##########        Function: applyConstantsToVector        ########## *)

(************************************************************************)
(* applyConstantsToVector[vector, lengthIndVar]                         *)
(* Purpose:  To apply constants k[i] through the components of a vector *)
(*           so that the indicies run from 1 to n through all           *)
(*           components.                                                *)
(* Input:    An expression in rational form                             *)
(*           The number of independent variables                        *)
(* Output:   The numerator of the expression                            *)
(* Created:  3 September, 2007 by DP                                    *)
(* Code is in File:  newhomotopyoper/hoaconve.m                         *)
(* Last Modified: 10 March, 2008, 14:56 by DP at CSM                    *)
(************************************************************************)

applyConstantsToVector[vector_, lengthIndVar_] :=
Module[{i1 = 1, lengthComp, homotopyResultList = {}, unknownConstants,
    homotopyResultPart, termsWithConstants, expnWithConstants, printACTV},

  If[debugApplyConstantsToVector, printACTV = Print, Clear[printACTV],
      Clear[printACTV]];
  printACTV["debug code: ACTV, Function: applyConstantsToVector,"<>
      " File: hoaconve.m"];

  printACTV["At ACTV IN, the vector given is:"];
  printACTV[vector];
  printACTV["At ACTV IN, the number of independent variables is ",
      lengthIndVar];

  (* Find the length of each component.                                 *)
  lengthComp[1] = If[Head[Part[vector, 1]] === Plus,
      Length[Part[vector, 1]], 1];
  lengthComp[2] = If[lengthIndVar > 1,
      If[Head[Part[vector, 2]] === Plus, Length[Part[vector, 2]], 1], 0];
  lengthComp[3] = If[lengthIndVar === 3,
      If[Head[Part[vector, 3]] === Plus, Length[Part[vector, 3]], 1], 0];
  printACTV["At ACTV 1, component 1 has length ", lengthComp[1], " "<>
      "component 2 has length ", lengthComp[2], ", and component 3 has"<>
      " length ", lengthComp[3]];

  (* Write all terms into a single list.                                *)
  While[i1 <= lengthIndVar,
    homotopyResultPart = If[Head[Part[vector, i1]] === Plus,
        Apply[List, Part[vector, i1]],
        Part[vector, i1]];
    printACTV["At ACTV  2, part ", i1, " of the homotopy result with "<>
        "", Length[homotopyResultPart], " terms has been converted to"<>
        " list form:"];
    printACTV[homotopyResultPart];
    homotopyResultList =
      Flatten[Append[homotopyResultList, homotopyResultPart]];
    printACTV["At ACTV 3, a full list of terms up to part ", i1,
      " of the homotopy result:"];
    printACTV[homotopyResultList];
    i1++
  ]; (* end While *)

  (* Set the list of unknown constants and match the constants with the  *)
  (* terms in the homotopy term list and with the ratings for each term. *)
  unknownConstants = Table[k[i], {i, 1, Length[homotopyResultList]}];

  (* Attach constants to the terms produced by the homotopy operator.    *)
  termsWithConstants = Thread[Times[homotopyResultList, unknownConstants]];
  printACTV["At ACTV 4, homotopy terms with their constants attached:"];
  printACTV[termsWithConstants];

  (* Remake the homotopy components from the list of homotopy terms with *)
  (* constants so that a divergence can be found.                        *)
  expnWithConstants[1] =
      Apply[Plus , Take[termsWithConstants, lengthComp[1]]];
  expnWithConstants[2] = If[lengthIndVar > 1, Apply[Plus,
      Take[termsWithConstants, {lengthComp[1] + 1,
      lengthComp[1] + lengthComp[2]}]], {}, {}];
  expnWithConstants[3] = If[lengthIndVar === 3, Apply[Plus,
      Take[termsWithConstants, {lengthComp[1] + lengthComp[2] + 1,
      Length[termsWithConstants]}]], {}, {}];
  printACTV["At ACTV 5, the homotopy result with constants is:"];
  printACTV[Flatten[{expnWithConstants[1], expnWithConstants[2],
      expnWithConstants[3]}]];

  Clear[i1, lengthComp, homotopyResultList, homotopyResultPart,
      termsWithConstants, printACTV];

  Return[{expnWithConstants[1], expnWithConstants[2],
          expnWithConstants[3], unknownConstants}]
]; (* end Module applyConstantsToVector *)

(* ##########        Function: expressionVarReassign         ########## *)

(************************************************************************)
(* expressionVarReassign[expression, indepVariableList]                 *)
(* Purpose:  To reassign independent and dependent variables to fit     *)
(*           the format used the all homotopy operator functions.       *)
(* Input:    An expression                                              *)
(*           A list of independent variables                            *)
(* Output:   An expression in the correct format to run in all          *)
(*           functions called by the homotopy operator.                 *)
(* Created:  17 March, 2006 by DP                                       *)
(* Code is in File:  newhomotopyoper/hoexvar2.m                         *)
(* Last Modified: 20 March, 2009, 8:35 by DP at CSM                     *)
(************************************************************************)

expressionVarReassign[expression_, indepVariableList_] :=
Module[{i, depVarList, expressionTerms, n0, n1, n2, n3, n4, integrVarList,
    newexpression, revisedExpression, depVarAssignRule, positionIntegVars,
    depVarAssignRuleInvert, lengthIndVarList = Length[indepVariableList],
    indepCheckList, indList, freeIndVars, newNumDepVars, reverseRule,
    printEVR},

  If[debugExpressionVarReassign, printEVR = Print, Clear[printEVR],
       Clear[printEVR]];
  printEVR["debug code: EVR, Function: expressionVarReassign,"<>
      " File: hoexvar2.m"];

  printEVR["At EVR IN, the expression given is:"];
  printEVR[expression];
  printEVR["At EVR IN, the independent variable list is ", indepVariableList];

  (* Modify input expression and variable list to match format used in  *)
  (* the homotopy operator functions.                                   *)

  (* Input Check 1:  Are the independent variables entered into a list  *)
  (* in the correct form?                                               *)
  If[Flatten[indepVariableList] =!= indepVariableList,
    Print["Independent variables must be entered as a single List, {a,b,c}."];
    Abort[]
  ]; (* end If Flatten[indepVariableList] *)

  (* Input Check 2:  Are there more than three independent variables    *)
  (* given for integration?                                             *)
  If[lengthIndVarList > 3,
    Print["This program can integrate expressions with respect to one"<>
        " independent variable, or invert a divergence on two or three"<>
        " independent variables. You have given ", lengthIndVarList, ""<>
        " independent variables, which is beyond the capabilities of"<>
        " this program."];
    Abort[]
  ]; (* end If lengthIndVarList *)

  (* Make lists of dependent variables and independent variables by     *)
  (* taking them from the input expression. At the same time check the  *)
  (* input for consistency to avoid obscure error messages in later     *)
  (* functions.                                                         *)

  (* Separate the terms in the given expression, then isolate functions.*)
  expressionTerms = Flatten[Expand[expression] /. Plus -> List /.
      Times -> List /. Power -> List //. _[n1_[n2__]] -> n1[n2]];
  expressionTerms = If[Head[expressionTerms] =!= List,
      {expressionTerms}, expressionTerms, expressionTerms];
  printEVR["At EVR 1, terms in the given expression put into a list:"];
  printEVR[expressionTerms];

  expressionTerms = Flatten[Map[FactorList, expressionTerms]];
  printEVR["At EVR 2, after mapping FactorList to each of the terms in the"<>
      " previous list:"];
  printEVR[expressionTerms];

  expressionTerms = Union[DeleteCases[expressionTerms,
     (n1_ /; NumericQ[n1] || Head[n1] === Symbol) | (_[n2_] /; NumericQ[n2])]];
  printEVR["At EVR 3, the list after taking out all non differential"<>
      " function terms:"];
  printEVR[expressionTerms];

  (* From expressionTerms, make a list of independent variables.        *)
  indepCheckList = Union[expressionTerms /.
      {n2_[n1__] /; n2 =!= List -> List[n1], List[List[n3__]] -> List[n3]}];
  printEVR["At EVR 4, a list of independent variables found in the"<>
      " expression:"];
  printEVR[indepCheckList];

  (* Input Check 3:  Are the independent variables put into the         *)
  (* functions in a consistent manner?                                  *)
  If[Length[indepCheckList] > 1,
    Print["The independent variables have not been placed in functions in"<>
        " a consistent manner.  That is, one function may show u[x,y] while"<>
        " another shows u[y,x] or u[y,z]."];
    Abort[]
  ]; (* end If Length[indepCheckList] *)

  indepCheckList = Flatten[indepCheckList];

  (* Input Check 4:  Do the independent variables given for integration *)
  (* match the independent variables in the expression?                 *)
  If[Not[Part[Union[Map[MemberQ[indepCheckList, #] &,
                        indepVariableList]], 1]],
    Print["Independent variables given for integration do not match the"<>
        " independent variables given in the differential function."];
    Abort[]
  ]; (* end If Not Union Map *)

  freeIndVars = Length[indepCheckList] - lengthIndVarList;
  printEVR["At EVR 5, there are ", freeIndVars, " free independent variables"<>
      " in the expression.  These variables exist in the expression,"<>
      " but are not being used in the integration process."];

  (* From expressionTerms, make a list of dependent variables.          *)
  indList = Apply[Sequence, indepCheckList];
  depVarList =
      Union[expressionTerms /. n0_[n1_[indList]] /; n0 =!= List -> n1 /.
      {Derivative[__][n2_][indList] -> n2, Derivative[__][n3_] -> n3,
      n4_[indList]  -> n4}];
  printEVR["At EVR 6, a list dependent variables found in the expression:"];
  printEVR[depVarList];

  (* Input Check 5:  Are there any variables appearing in both the      *)
  (* independent variable list and the dependent variable list at       *)
  (* the same time?                                                     *)
  If[Intersection[indepCheckList, depVarList] =!= {},
    Print["The variable(s), ", Intersection[indepCheckList, depVarList],
        " occur as both independent and dependent. The program does not"<>
        " know how to interpret this."];
    Abort[]
  ]; (* end If Complement *)

  (* Input Check 6:  Check that the order of the variables given for    *)
  (* integration matches the order they are found in the independent    *)
  (* variable list.                                                     *)
  positionIntegVars = Flatten[Map[Position[indepCheckList, #] &,
      indepVariableList], 1];
  printEVR["At EVR 7, the positions that the integrating independent"<>
      " variables are found in the full independent variable list:"];
  printEVR[positionIntegVars];

  If[Sort[positionIntegVars] =!= positionIntegVars,
    integrVarList = Extract[indepCheckList, Sort[positionIntegVars]],
  (* else *)
    integrVarList = indepVariableList
  ]; (* end If Sort *)
  printEVR["At EVR 8, the list of variables given for integration is now:"];
  printEVR[integrVarList];

  (* Build a list of rules reassigning the dependent variables.          *)
  (* The expression with reassigned variables will be passed to          *)
  (* other functions in the process of running the homotopy operator.    *)
  depVarAssignRule = Table[Part[depVarList, i] -> u[i],
      {i, 1, Length[depVarList]}];
  depVarAssignRule = Apply[Rule, DeleteCases[Apply[List, depVarAssignRule, 1],
      {_, _}?(Part[#, 1] === Part[#, 2] &)], 1];
  depVarAssignRuleInvert =
      Table[Reverse[Part[depVarAssignRule, i]], {i, 1,
          Length[depVarAssignRule]}];
  printEVR["At EVR 9, a list of rules for reassigning dependent variables:"];
  printEVR[depVarAssignRule];
  printEVR["At EVR 10, list of rules to undo the reassignment of dependent"<>
      " variables:"];
  printEVR[depVarAssignRuleInvert];

  newexpression = If[depVarAssignRule =!= {}, expression /. depVarAssignRule,
      expression, expression];
  printEVR["At EVR 11, the expression with the dependent variables replaced:"];
  printEVR[newexpression];
  printEVR["At EVR 12, the expression to be used from this point on:"];
  printEVR[newexpression];

  (* Rewrite the expression, eliminating any free independent variables *)
  (* reassigning terms with derivatives on the free independent         *)
  (* variables, so they do not interfere with the homotopy operator.    *)
  If[freeIndVars > 0,
    {revisedExpression, newNumDepVars, reverseRule} =
        changeIndependentParameters[newexpression, Length[depVarList],
        indepCheckList, integrVarList];
    printEVR["At EVR 13, the expression after passing through"<>
        " removeFreeIndpendentVariables is:"];
    printEVR[revisedExpression];
    printEVR["At EVR 14, there are now ", newNumDepVars, " dependent"
        " variables in the expression."];
    printEVR["At EVR 15, the rule(s) for reversing the renaming of terms"<>
        " containing derivatives on auxillary variables:"];
    printEVR[reverseRule],

  (* else *)
    revisedExpression = newexpression;
    newNumDepVars = Length[depVarList];
    reverseRule = {};
  ]; (* end If freeIndVars > 0 *)

  printEVR["At EVR OUT, the expression in appropriate form is:"];
  printEVR[revisedExpression];

  Clear[i, depVarList, expressionTerms, n0, n1, n2, n3, n4, newexpression,
      depVarAssignRule, positionIntegVars, lengthIndVarList, indList,
      freeIndVars,  printEVR];

  Return[{revisedExpression, depVarAssignRuleInvert, newNumDepVars,
      indepCheckList, integrVarList, reverseRule}]
]; (* end Module expressionVarReassign *)

(* ##########      Function: changeIndependentParameters     ########## *)

(************************************************************************)
(* changeIndependentParameters[expression, depVarList, indepVarList,    *)
(*                  integVarList]                                       *)
(* Purpose: To rewrite the given expression in terms of the variable(s) *)
(*          required in the integration process.  New dependent         *)
(*          variables will replace derivatives on variables not         *)
(*          involved in the integration process.                        *)
(*          e. g. If D[u[1][x,y],x,y] is to be integrated with respect  *)
(*          to x only, the expression will become D[u[2][x],x] with     *)
(*          u[2][x] = D[u[1][x,y],y]                                    *)
(* Input:   An expression                                               *)
(*          A list of dependent variables                               *)
(*          A list of all independent variables found in the given      *)
(*          expression                                                  *)
(*          A list of variables involved in integration                 *)
(* Output:  A revised expression.                                       *)
(*          A revised number of dependent variables                     *)
(*          Rules for rewriting the result in the same form as given    *)
(* Created: 10 July, 2007 by DP                                         *)
(* Code is in File:  newhomotopyoper/hochinpa.m                         *)
(* Last Modified: 10 March, 2008, 14:49 by DP at CSM                    *)
(************************************************************************)

changeIndependentParameters[expression_, noDepVars_, indepVarList_,
    integrVarList_] :=
Module[{listOfDiffTerms, lengthAllIndVars = Length[indepVarList],
    positionFreeVars, patternForFreeVars, patternForAllIndVars, newexpression,
    patternForIntegrVars, patternFreeVarsResult, patternIntegrVarsResult,
    patternForFreeDerivs, pattern, set, dd1, dd2, m1, m2, j,
    noOfDerivsListFreeVars, noOfDerivsListIntegrVars, conditions,
    termsWDerivsOnFreeIndVars, termsWDerivsSeparated, termsWFreeIndVarDers,
    ruleForFreeIndVarTerms, termsWDerivsOnIntegrVars, ruleToChangeExpn,
    freeIndVars, numNewDepVars, ruleToReverseChange, printCIP},

  If[debugChangeIndependentParameters, printCIP = Print, Clear[printCIP],
      Clear[printCIP]];
  printCIP["debug code: CIP, Function: changeIndependentParameters,"<>
      " File: hochinpa.m"];

  printCIP["At CIP IN, the given expression is:"];
  printCIP[expression];
  printCIP["At CIP IN, the number of dependent variables is ", noDepVars];
  printCIP["At CIP IN, the list of independent variables found in the"<>
      " expression is ", indepVarList];
  printCIP["At CIP IN, a list of variables to be used in the integration"<>
      " process is ", integrVarList];

  (* Make a list of differential terms in the expression.               *)
  listOfDiffTerms = Union[Cases[expression,
      _[___, Part[integrVarList, 1], ___], {0, Infinity}]];
  printCIP["At CIP 1, a list of differential terms found in the"<>
      " expression:"];
  printCIP[listOfDiffTerms];

  (* Make a list of free independent variables in the expression.       *)
  freeIndVars = Complement[indepVarList, integrVarList];

  (* Find the positions of the free independent variables in the full   *)
  (* independent variable list, then rewrite the list of free variables *)
  (* in the order in which they were found.                             *)
  positionFreeVars =
          Flatten[Map[Position[indepVarList, #] &, freeIndVars], 1];
  freeIndVars = Extract[indepVarList, Sort[positionFreeVars]];
  printCIP["At CIP 2, the free variables are found at these positions"<>
      " in the full independent variable list:"];
  printCIP[positionFreeVars];
  printCIP["At CIP 3, a list of free independent variables"<>
      " (not involved in the integration process):"];
  printCIP[freeIndVars];

  (* Create pattern maps for rules so that free variables can be        *)
  (* removed and terms containing derivative on free variables can be   *)
  (* renamed and kept track of.                                         *)
  patternForFreeVars = Table[If[MemberQ[positionFreeVars, {j}],
      pattern[ToExpression["n" <> ToString[j]], Blank[]], Blank[]],
      {j, 1, lengthAllIndVars}];
  printCIP["At CIP 4, the pattern for the free variables found"<>
      " the full independent variable list:"];
  printCIP[patternForFreeVars];

  patternForAllIndVars = Table[pattern[ToExpression[
      "n" <> ToString[j]], Blank[]], {j, 1, lengthAllIndVars}];
  patternForAllIndVars = patternForAllIndVars /. pattern -> Pattern;
  printCIP["At CIP 5, the pattern for all independent variables:"];
  printCIP[patternForAllIndVars];

  patternForIntegrVars = Table[If[! MemberQ[positionFreeVars, {j}],
      pattern[ToExpression["n" <> ToString[j]], Blank[]], Blank[]],
      {j, 1, lengthAllIndVars}];
  printCIP["At CIP 6, the pattern for all independent variables found"<>
      " in the integration list:"];
  printCIP[patternForIntegrVars];

  patternFreeVarsResult = DeleteCases[patternForFreeVars /.
      pattern -> Sequence, Verbatim[Blank[]]];
  printCIP["At CIP 7, the result after the pattern for the free variables"<>
      " has been applied:"];
  printCIP[patternFreeVarsResult];

  patternIntegrVarsResult = DeleteCases[patternForIntegrVars /.
      pattern -> Sequence, Verbatim[Blank[]]];
  printCIP["At CIP 8, the result after the pattern for the integrating"<>
      " variables has been applied:"];
  printCIP[patternIntegrVarsResult];

  noOfDerivsListFreeVars = Thread[List[freeIndVars, patternFreeVarsResult]];
  printCIP["At CIP 9, each free derivative can be placed with the number"<>
      " of derivatives operating on it:"];
  printCIP[noOfDerivsListFreeVars];

  noOfDerivsListIntegrVars = Thread[List[integrVarList,
      patternIntegrVarsResult]];
  printCIP["At CIP 10, each derivative on an integrating variable can be"<>
      " placed with the number of derivatives operating on it:"];
  printCIP[noOfDerivsListIntegrVars];

  patternForFreeDerivs = patternForFreeVars /. pattern -> Pattern;
  printCIP["At CIP 11, the pattern to be placed in Derivative to separate"<>
      " derivatives on free independent variables:"];
  printCIP[patternForFreeDerivs];

  patternForIntegrVars = patternForIntegrVars/. pattern -> Pattern;
  printCIP["At CIP 12, the pattern to be placed in Derivative to separate"<>
      " derivatives on the integrating variables:"];
  printCIP[patternForIntegrVars];

  conditions = Apply[Or, Map[Greater[#, 0] &, patternFreeVarsResult]];
  printCIP["At CIP 13, check for cases in Derivative where derivatives"<>
      " on free variables exist:"];
  printCIP[conditions];

  termsWDerivsOnFreeIndVars = Union[Cases[listOfDiffTerms,
      Derivative[Apply[Sequence, patternForFreeDerivs]][u[_]][__] /;
      Evaluate[conditions]]];
  printCIP["At CIP 14, a list of differential terms that contain"<>
      " derivatives on free independent variables:"];
  printCIP[termsWDerivsOnFreeIndVars];

 (* Separate the derivatives on free independent variables from the     *)
 (* part with derivatives on integrating variables.                     *)
  termsWDerivsSeparated = termsWDerivsOnFreeIndVars /.
      Derivative[Apply[Sequence, patternForAllIndVars]][u[m1_]][m2__] ->
      dd1[dd2[u[m1][m2], Apply[Sequence, noOfDerivsListFreeVars]],
      Apply[Sequence, noOfDerivsListIntegrVars]] /. dd2 -> D;
  printCIP["At CIP 15, derivatives on integrating variables have been"<>
      " separated from derivatives on free independent variables:"];
  printCIP[termsWDerivsSeparated];

  termsWFreeIndVarDers = Union[Map[Part[#, 1] &, termsWDerivsSeparated]];
  printCIP["At CIP 16, the terms with derivatives on free independent"<>
      " variables have been isolated and contain only the free variable"<>
      " derivatives:"];
  printCIP[termsWFreeIndVarDers];

  (* Reassign the parts with derivatives on free variables to new       *)
  (* dependent variables.                                               *)

  numNewDepVars = Length[termsWFreeIndVarDers];
  ruleForFreeIndVarTerms = Table[set[u[j][Apply[Sequence, indepVarList]],
      Part[termsWFreeIndVarDers, j - noDepVars]],
      {j, noDepVars + 1, numNewDepVars + noDepVars}];
  ruleForFreeIndVarTerms = ruleForFreeIndVarTerms /.
      set[m1_, m2_] -> Rule[m2, m1];
  printCIP["At CIP 17, the rule for assigning terms with derivatives"<>
      " on free independent variables to new dependent variables:"];
  printCIP[ruleForFreeIndVarTerms];

  termsWDerivsOnIntegrVars = termsWDerivsSeparated /.
      ruleForFreeIndVarTerms /. dd1 -> D;
  printCIP["At CIP 18, all terms have derivatives with respect to"<>
      " integrating variables only:"];
  printCIP[termsWDerivsOnIntegrVars];

  ruleToChangeExpn = Thread[Rule[termsWDerivsOnFreeIndVars,
      termsWDerivsOnIntegrVars]];
  printCIP["At CIP 19, rules for reassgning the terms in the expression:"];
  printCIP[ruleToChangeExpn];

  newexpression = expression /. ruleToChangeExpn;
  printCIP["At CIP 20, the expression after new dependent variables"<>
       " have replaced all terms containing derivatives on free variables:"];
  printCIP[newexpression];

  newexpression = newexpression /. {u[m1_][Apply[Sequence, indepVarList]] ->
          u[m1][Apply[Sequence, integrVarList]],
      Derivative[Apply[Sequence, patternForIntegrVars]][u[m1_]]
          [Apply[Sequence, indepVarList]] ->
      Derivative[Apply[Sequence, patternIntegrVarsResult]][u[m1]]
          [Apply[Sequence, integrVarList]]};
  printCIP["At CIP 21, the expression after all free variables have been"<>
      " removed:"];
  printCIP[newexpression];

  ruleToReverseChange = If[ruleForFreeIndVarTerms =!= {},
      ruleForFreeIndVarTerms /. Rule[m1_, m2_] -> Rule[m2, m1],
      {{}}];
  printCIP["At CIP 22, the rule(s) required to undo the changes made"<>
    " to the expression in removeFreeIndpendentVariables:"];
  printCIP[ruleToReverseChange];

  Clear[listOfDiffTerms, lengthAllIndVars, positionFreeVars,
      patternForFreeVars, patternForAllIndVars, patternForIntegrVars,
      patternFreeVarsResult, patternIntegrVarsResult, patternForFreeDerivs,
      pattern, set, dd1, dd2, m1, m2, j, noOfDerivsListFreeVars,
      noOfDerivsListIntegrVars, conditions, termsWDerivsOnFreeIndVars,
      termsWDerivsSeparated, termsWFreeIndVarDers, ruleForFreeIndVarTerms,
      termsWDerivsOnIntegrVars, ruleToChangeExpn, freeIndVars, printCIP];

  Return[{newexpression, numNewDepVars + noDepVars, ruleToReverseChange}]
]; (* end Module changeIndependentParameters *)

(* ##########      Function: returnAuxIndependentVars        ########## *)
(************************************************************************)
(* returnAuxIndependentVars[homotopyResult, fullIndVarList,             *)
(*                          integrateVarList, reverseRules]             *)
(* Purpose:  To replace renamed dependent variables with the derivative *)
(*           on auxillary variables removed in                          *)
(*           removeFreeIndpendentVariables, and reinsert the auxillary  *)
(*           variables back into the expression, returning it to its    *)
(*           original form.                                             *)
(* Input:    The result returned by the homotopy operator               *)
(*           A list of all independent variables original expression    *)
(*           A list of rules for reversing the changes made in the      *)
(*               function removeFreeIndpendentVariables                 *)
(* Output:   A revised expression.                                      *)
(* Created:  12 July, 2007 by DP                                        *)
(* Code is in File:  newhomotopyoper/hochinpa.m                         *)
(* Last Modified: 10 March, 2008, 13:04 by DP at CSM                    *)
(************************************************************************)

returnAuxIndependentVars[homotopyResult_, fullIndVarList_, integrateVarList_,
    reverseRules_] :=
Module[{positionIntegVars, patternForIntegVars, initialPattern, finalPattern,
    newExpression, addDepVarNos, auxDepTerms, auxDepTermsRev, ruleAuxToOldForm,
    pattern, dd1, dd2, j, m1, m2, m3, thread, list, printRAIV},

  If[debugReturnAuxIndependentVars, printRAIV = Print, Clear[printRAIV],
      Clear[printRAIV]];
  printRAIV["debug code: RAIV, Function: returnAuxIndependentVars,"<>
      " File: hochinpa.m"];

  printRAIV["At RAIV IN, the result from the homotopy operator:"];
  printRAIV[homotopyResult];
  printRAIV["At RAIV IN, a full list of independent variables is ",
      fullIndVarList];
  printRAIV["At RAIV IN, a list of independent variables used in the"<>
      " integration process is ", fullIndVarList];
  printRAIV["At RAIV IN, a list of rules for reversing the changes made"<>
      " when removing the auxillary independent variables is:"];
  printRAIV[reverseRules];

  (* Put all auxillary independent variables back into the homotopy result. *)
  positionIntegVars = Flatten[Map[Position[fullIndVarList, #] &,
      integrateVarList], 1];
  printRAIV["At RAIV 1, the positions of the variables used in the"<>
      " integration process as they occur in the full variable list:"];
  printRAIV[positionIntegVars];

  patternForIntegVars = Table[If[MemberQ[positionIntegVars, {j}],
      pattern[ToExpression["n" <> ToString[j]], Blank[]], 0],
      {j, 1, Length[fullIndVarList]}];
  printRAIV["At RAIV 2, setting up a pattern so that the derivatives are"<>
      " placed in the correct position with respect to the full independent"<>
      " variable list:"];
  printRAIV[patternForIntegVars];

  initialPattern = DeleteCases[patternForIntegVars, 0] /. pattern -> Pattern;
  printRAIV["At RAIV 3, derivative numbers will be picked up according to"<>
      " this list:"];
  printRAIV[initialPattern];

  finalPattern = DeleteCases[patternForIntegVars /. pattern -> Sequence,
      Verbatim[Blank[]]];
  printRAIV["At RAIV 4, derivative numbers will be placed according to"<>
      " this list:"];
  printRAIV[finalPattern];

  newExpression = homotopyResult /.
      {Derivative[Apply[Sequence, initialPattern]][u[m1_]][__] ->
      Derivative[Apply[Sequence, finalPattern]][u[m1]]
          [Apply[Sequence, fullIndVarList]],
      u[m2_][__] -> u[m2][Apply[Sequence, fullIndVarList]]};
  printRAIV["At RAIV 5, the result from the homotopy operator with all"<>
      " independent variables:"];
  printRAIV[newExpression];

  (* Apply reverse rules to rename terms with derivatives on auxillary  *)
  (* variables.                                                         *)
  If[reverseRules =!= {{}},
    addDepVarNos = Map[Part[#, 1] &, reverseRules] /. u[m1_][__] -> m1;
    printRAIV["At RAIV 6, a list of numbers on dependent variables to"<>
        " replace terms with derivatives on auxillary variables:"];
    printRAIV[addDepVarNos];

    auxDepTerms = Cases[newExpression, Derivative[__][u[m1_]][__] /;
        MemberQ[addDepVarNos, m1], {0, Infinity}];
    printRAIV["At RAIV 7, all terms that have auxillary dependent variables:"];
    printRAIV[auxDepTerms];

    (* Restructure the auxillary dependent variable terms so that their *)
    (* derivatives are included in the original form.                   *)
    auxDepTermsRev = auxDepTerms /. Derivative[m1__][m2_][m3__] ->
        dd2[m2[m3], thread[list[{m3}, {m1}]]] /. thread -> Thread /.
        list -> List;
    printRAIV["At RAIV 8, after restructuring the derivative format"<>
        " on the auxillary dependent variable terms:"];
    printRAIV[auxDepTermsRev];

    auxDepTermsRev = auxDepTermsRev /. reverseRules;
    printRAIV["At RAIV 9, the auxillary dependent terms after they"<>
        " have been translated back to the original dependent variables"<>
        " with derivatives on auxillary variables:"];
    printRAIV[auxDepTerms];

    auxDepTermsRev = ReleaseHold[auxDepTermsRev /.
        dd2[m1__, m2_] -> dd1[m1, Hold[Apply[Sequence, m2]]]] /. dd1 -> D;
    auxDepTerms = auxDepTerms /. reverseRules;
    printRAIV["At RAIV 10, after all derivatives on the auxillary"<>
        " dependent terms have been applied:"];
    printRAIV[auxDepTerms];

    ruleAuxToOldForm = Thread[Rule[auxDepTerms, auxDepTermsRev]];
    printRAIV["At RAIV 11, final replacement rules to return terms to the"<>
        " form given in the original expression:"];
    printRAIV[ruleAuxToOldForm];

    newExpression = newExpression /. ruleAuxToOldForm /. reverseRules;
    printRAIV["At RAIV 12, the result given by the homotopy operator"<>
        " in the form given in the original expression:"];
    printRAIV[newExpression]
  ]; (* end If reverseRules *)

  Clear[positionIntegVars, patternForIntegVars, initialPattern, finalPattern,
    addDepVarNos, auxDepTerms, auxDepTermsRev, ruleAuxToOldForm, pattern,
    dd1, dd2, j, m1, m2, m3, thread, list, printRAIV];

  Return[newExpression]
]; (* end Module returnAuxIndependentVars *)

(* ##########          Function: integrabilityTest           ########## *)

(************************************************************************)
(* integrabilityTest1D[expression, dependentVar, listOfIndepVars,       *)
(*                     (optional) parameters]                           *)
(* Purpose:  To determine if the expression is integrable by applying   *)
(*           the zero-Euler operator.  If the expression is not         *)
(*           integrable and parameters have been given, then the        *)
(*           function will try to identify values for the parameters    *)
(*           to make the expression integrable.                         *)
(* Input:    An expression                                              *)
(*           A number for the dependent variable                        *)
(*           A list of independent variables                            *)
(*           A list of parameters (optional)                            *)
(* Output:   A  partial derivative of the expression with respect to    *)
(*           terms over the jet space                                   *)
(*           The order of the expression with respect to the            *)
(*           independent variable                                       *)
(* Created:  22 March, 2006 by DP                                       *)
(* Code is in File:  newhomotopyoper/hointest.m                         *)
(* Last Modified: 13 May, 2009 at 13:26 by DP                           *)
(************************************************************************)

integrabilityTest[expression_, noOfDepVar_, listOfIndepVars_,
    parameters_:{}] :=
Module[{j, variationalDerivativeResult, testBasis, parameterSolutions,
    parameterList, intelligibleAnswer1 = 0, intelligibleAnswer2 = 0,
    paramList, answer1, answer2, answer3, newparamvalues = {}, secondTest,
    printIT, lthiv, freeVars, lps, lpm, lpa, lpf, lpmm, lpaa, lpff},

  If[debugIntegrabilityTest, printIT = Print, Clear[printIT], Clear[printIT]];
  printIT["debug code: IT, Function: integrabilityTest,"<>
      " File: integrabilityTest.m"];

  printIT["At IT IN, the given expression is:"];
  printIT[expression];
  printIT["At IT IN, the number of dependent variables is ", noOfDepVar];
  printIT["At IT IN, the list of independent variables is ", listOfIndepVars];
  printIT["At IT IN, parameters given to the function include:"];
  printIT[parameters];

  (* Check parameter list for proper format.                            *)
  parameterList = If[Length[parameters] === 0 && parameters =!= {},
      {parameters}, parameters, parameters];
  printIT["At IT 1, the list of parameters given for this problem is"];
  printIT[parameterList];

  (* Find the zero-Euler operator for the expression for each        *)
  (* component of the dependent variable.                               *)
  variationalDerivativeResult = Table[variationalDerivativeMultiD[
      expression, j, listOfIndepVars], {j, 1, noOfDepVar}];
  printIT["At IT 2, results from applying the zero-Euler operator to"<>
      " each component:"];
  printIT[variationalDerivativeResult];

  testBasis = Table[0, {j, 1, noOfDepVar}];
  If[variationalDerivativeResult === testBasis,
    printIT["At IT 3, the zero-Euler operator applied to the expression"<>
        " is zero for all components.  This expression is integrable,"<>
        " and the integration process will begin immediately."];

    If[globalHomotopyVerbose,
      lthiv = Length[listOfIndepVars];
      printOUT["The zero-Euler operator applied to the given expresion has"<>
        " produced zero for all dependent variables."<>
        "  The expression is a "<> If[lthiv === 1, "total derivative.",
          "divergence."] <>"  The "<> If[lthiv === 1, "total derivative",
          "divergence"] <>" will be "<> If[lthiv === 1, "integrated",
          "inverted"] <>" with respect to "<> If[lthiv === 1, Apply[ToString,
          listOfIndepVars], StringReplace[ToString[Insert[listOfIndepVars,
          "and", lthiv]], {"{" -> "", "}" -> "", ", and," -> " and"}]] <>"."];
      printOUT["------------------------------------------"];
    ]; (* end globalHomotopyVerbose *)

    Return[{"Integrable", {}}]
  ]; (* end If variationalDerivativeResult *)

  If[variationalDerivativeResult =!= testBasis && parameterList === {},
    printOUT["The zero-Euler operator applied to the expression has not",
        " produced zero for all dependent variables.  The expression is",
        " not integrable. The integration process will not be activated."];
    Clear[j, variationalDerivativeResult, testBasis, parameterSolutions,
        parameterList, printIT, lthiv];
    Return[{"Not Integrable", {}}]
  ]; (* end If variationalDerivativeResult *)

  If[variationalDerivativeResult =!= testBasis && parameterList =!= {},
    printIT["At IT 4, looking for values for the parameters that were given."];
    parameterSolutions = parameterTest[variationalDerivativeResult,
        listOfIndepVars, parameters];
    printIT["At IT 5, the list of solutions given by the variational"<>
        " derivative is:"];
    printIT[parameterSolutions];
    If[parameterSolutions === {}, Print["There are no values for the"<>
        " parameter(s) for which the expression is integrable.  The"<>
        " integration process will not be activated."];
      Return[{"Not Integrable", {}}]
    ]; (* end If parameterSolutions *)

    lps = If[Length[parameterSolutions] > 1, "values or expressions",
        "a value or an expression"];
    lpm = If[Length[parameterList] > 1, "parameters", "parameter"];
    lpa = If[Length[parameterList] > 1, "are", "is"];
    lpf = If[Length[parameterList] > 1, "free variables", "a free variable"];
    Print["The zero-Euler operator applied to the expression has not"<>
        " produced zero for all dependent variables.  The expression may"<>
        " become integrable if the "<> lpm <>" "<> lpa <>" set to the"<>
        " value(s) given in the following list:"];
    Print[parameterSolutions];
    freeVars = Complement[parameterList,
        Union[Map[First, parameterSolutions]]];
    If[freeVars =!= {},
      lpmm = If[Length[freeVars] > 1, "parameters", "parameter"];
      lpaa = If[Length[freeVars] > 1, "are", "is"];
      lpff = If[Length[freeVars] > 1, "free variables", "a free variable"];
      Print["The "<>lpmm<>" ", freeVars, " "<>lpaa<>" "<>lpff<>
          " and may be set to any value."]
    ]; (* end If freeVars *)

    While[intelligibleAnswer1 === 0,
      answer1 = ToString[Input["See the information about parameters on"<>
          " the main screen.  The program suggests "<> lps <>" for the "<>
          lpm <>" given in the list, "<> ToString[parameterList] <>
          ", which may allow the integration process to continue.  Do you"<>
          " want to set "<> lps <>" for the "<> lpm <>" and continue"<>
          " with the problem?  (Yes or No)"]];

      If[ToLowerCase[ToString[answer1]] === "qq", Abort[]];

      If[MemberQ[{"N", "NO", "NOT", "NONE", "NN"}, ToUpperCase[answer1]],
        intelligibleAnswer1 = 1;
        Print["The expression is not integrable. The integration process"<>
            " will not be activated."];
        Clear[j, variationalDerivativeResult, testBasis, parameterSolutions,
            parameterList, intelligibleAnswer1, intelligibleAnswer2,
            paramList, answer1, answer2, answer3, newparamvalues, secondTest,
            printIT, lthiv, freeVars, lps, lpm, lpa, lpg];
        Return[{"Not Integrable", {}}]
      ]; (* end If answer1 === No *)

      If[MemberQ[{"YES", "Y", "YE", "YY", "YS", "YEA" "YEP"},
          ToUpperCase[answer1]],
        intelligibleAnswer1 = 1;
        paramList = parameterList;
        printIT["At IT 6, the list of parameters to be revised:"];
        printIT[parameterList];
        While[paramList =!= {},
          currentParam = First[paramList];
          answer2 = Input["If you wish to change the parameter {"<>
              ToString[currentParam] <>"} enter a value or expression"<>
              " for the parameter. The value or expression should come"<>
              " from the solution list on the main screen.  \n\nIf you do"<>
              " not wish to change the parameter "<> ToString[currentParam] <>
              ", then enter  "<> ToString[currentParam] <>"."];
          Print["The parameter ", currentParam, " has been reset to ",
              currentParam, " = ", answer2];
          answer3 = Input["You have set "<> ToString[currentParam] <> " to"<>
              " the value or expression shown in the main notebook.  Enter"<>
              " yes if you wish to keep "<> ToString[currentParam] <>" as"<>
              " it is, or enter no if you wish to make further changes to "<> 
              ToString[currentParam] <> ".  Enter qq if you wish to quit."];
          If[ToLowerCase[ToString[answer3]] === "qq",
            Print["The expression is not integrable. The integration process"<>
                " will not be activated."];
            Clear[j, variationalDerivativeResult, parameterSolutions,
                intelligibleAnswer1, intelligibleAnswer2, secondTest, 
                paramList, answer1, answer2, answer3, newparamvalues,
                printIT, lthiv, lps, lpm, lpf, testBasis, parameterList, lpa];
          Return[{"Not Integrable", {}}]
          ]; (* end If ToLowerCase[answer3] *)
         If[MemberQ[{"y", "yy", "yyy", "yes"}, ToLowerCase[ToString[answer3]]],
            paramList = Drop[paramList, 1];
            printIT["At IT 7, the list of parameters to be revised:"];
            printIT[paramList];
            newparamvalues = Append[newparamvalues,
                Rule[currentParam, answer2]];
            printIT["At IT 8, the list of rules for revising the parameters:"];
            printIT[newparamvalues]
          ] (* end If MemberQ *)
        ] (* end While parlist *)
      ] (* end If MemberQ === Yes *)
    ]; (* end While intelligilbeAnswer1 *)

    printOUT["The values given are  ", TableForm[newparamvalues,
        TableDirections -> {Row, Column}]];
    printOUT["After a recheck for integrability, the integration process will",
        " begin."];
    secondTest = integrabilityTest[expression /. newparamvalues,
        noOfDepVar, listOfIndepVars]
  ]; (* end If variationalDerivativeResult *)

  Part[secondTest, 2] = newparamvalues;

  Clear[j, variationalDerivativeResult, testBasis, parameterSolutions,
      parameterList, intelligibleAnswer1, intelligibleAnswer2,
      paramList, answer1, answer2, answer3, newparamvalues,
      printIT, lthiv, freeVars, lps, lpm, lpa, lpf];

  Return[secondTest]
]; (* end Module integrabilityTest *)

(* ##########            Function: parameterTest             ########## *)

(************************************************************************)
(* parameterTest[variationalDerivativeResult, dependentVar, parameters] *)
(* Purpose:  To try to identify values for the parameters that will     *)
(*           make the expression integrable.                            *)
(* Input:    An expression                                              *)
(*           A number for the dependent variable                        *)
(*           A list of parameters                                       *)
(* Output:   Suggested values for parameters, or a statement that no    *)
(*           such parameter exists.                                     *)
(* Created:  22 March, 2006 by DP                                       *)
(* Code is in File:  newhomotopyoper/hopatest.m                         *)
(* Last Modified: 19 March, 2009, 17:31 by DP at CSM                    *)
(************************************************************************)

parameterTest[variationalDerivative_, indepvars_, parameterList_] :=
Module[{varDerList, whereIsU, termsNoU, termsWithU, positionFalse, coefEqns,
      extraDepVars,removeRule, parameterSolutions, coef, ext, cmpl, plus,
      printPT},

  If[debugParameterTest, printPT = Print, Clear[printPT], Clear[printPT]];
  printPT["debug code: PT, Function: parameterTest,"<>
      " File: parameterTest.m"];

  printPT["At PT IN, the variational derivative given is:"];
  printPT[variationalDerivative];
  printPT["At PT IN, the list of independent variables is ", indepvars];
  printPT["At PT IN, a list of parameters given by the user:"];
  printPT[parameterList];

  (* Make a list of terms in for each component of the variational      *)
  (* derivative.                                                        *)
  varDerList = Map[List, Expand[variationalDerivative], {1}] /.
      Plus -> Sequence;
  printPT["At PT 1, a list of terms in the variational derivative:"];
  printPT[varDerList];

  (* Separate terms in the list of variational derivative terms to that *)
  (* constants are separated from all other terms.                      *)
  whereIsU = Map[MemberQ[#, u[_][__] | Derivative[__][u[_]][__],
      {0, Infinity}] &, varDerList, {2}];
  printPT["At PT 2, True corresponds to a term in the variational"<>
      " derivative that contains u[_], False corresponds to a term with"<>
      " no u[_]: "];
  printPT[whereIsU];

  positionFalse = Map[Position[#, False] &, whereIsU];
  termsNoU = Thread[ext[varDerList, positionFalse]] /. ext -> Extract;
  printPT["At PT 3, all terms that do not have u[_]:"];
  printPT[termsNoU];

  termsWithU = Apply[Plus, Thread[cmpl[varDerList, termsNoU]] /. 
      cmpl -> Complement, {1}];
  zeroCheck = Map[SameQ[#, 0] &, termsWithU];
  Table[If[Part[zeroCheck, i1], 
    Part[termsWithU, i1] = u[10000][Apply[Sequence, indepvars]], 
    Part[termsWithU, i1] = Part[termsWithU, i1]], {i1, 1, Length[varDerList]}];
  printPT["At PT 4, all terms in the variational derivative that have u[_]:"];
  printPT[termsWithU];

  (* Set up a system of coefficient equations where the coefficients    *)
  (* contain all parameters.                                            *)
  strippedterms = Map[Union, Map[homotopyStripper[#, indepvars] &,
      termsWithU]];
  printPT["At PT 5, a list stripped terms in the variational derivative: "];
  printPT[strippedterms];

  (* All coefficients are grouped into expressions.                     *)
  coefEqns = Flatten[Thread[coef[Expand[variationalDerivative],
      strippedterms]] /. coef -> Coefficient];
  printPT["At PT 6, a list of coefficients for every term in the"<>
     " stripped list:"];
  printPT[coefEqns];

  (* Any extra u[_]s are removed from coefficient expressions.          *)
  coefEqns = DeleteCases[coefEqns, n1_ /; MemberQ[Denominator[n1], 
      u[_][__] | Derivative[__][u[_]][__], {0, Infinity}], {1, Infinity}];
  printPT["At PT 7, all denominators with u[_]s have been removed from the"<>
      "coefficient expressions:"];
  printPT[coefEqns];

  (* Any extra u[_]s are removed from coefficient expressions.          *)
  extraDepVars = Union[Flatten[Table[Cases[coefEqns,
      Nest[_, indepvars, i1] /. List -> Sequence,
      {0, Infinity}], {i1, 1, 4}]]];
  printPT["At PT 8, all u[_]s and functions of u[_] in the coefficient"<>
      " expressions:"];
  printPT[extraDepVars];
  
  removeRule = Map[rule[#, 0] &, extraDepVars] /. rule -> Rule;
    printPT["At PT 9, the set of rules for removing any extra u[_]s: "];
    printPT[removeRule];

  Off[Power::infy];
  coefEqns = Flatten[coefEqns /. Plus -> plus /. removeRule];
  printPT["At PT 10, all u[_] are removed from the coefficient"<>
      " expressions.  Note: any u[_] in denominators are replaced by"<>
      " ComplexInfinity:"];
  printPT[coefEqns];
  On[Power::infy];
  
  coefEqns = DeleteCases[coefEqns, ComplexInfinity, {0, Infinity}] /.
      plus -> Plus;
  printPT["At PT 11, all terms containing Complex Infinity have been"<>
      " removed from the coefficient expressions"];
  printPT[coefEqns];

  coefEqns = Union[coefEqns, Apply[Plus, termsNoU, {1}]];
  printPT["At PT 12, a list of all coefficient expressions from the"<>
      " variational derivative:"];
  printPT[coefEqns];

  parameterSolutions = Solve[Map[Equal[#, 0] &, coefEqns], parameterList];
  printPT["At PT 13, the solution to the coefficient equations:"];
  printPT[parameterSolutions];

  Clear[varDerList, whereIsU, termsNoU, termsWithU, coefEqns, extraDepVars,
     positionFalse, removeRule, coef, ext, cmpl, plus, printPT];

  Return[Union[Flatten[parameterSolutions]]]
]; (* end Module parameterTest *)

(* ##########           Function: degreeOperator             ########## *)

(************************************************************************)
(* degreeOperator[expression, dependentVar, IndepVars]                  *)
(* Purpose:  To determine the overall degree of a differential function *)
(* Input:    An expression                                              *)
(*           The number of dependent variables                          *)
(*           A list of independent variables                            *)
(* Output:   A revised expression with information about its degree     *)
(*           contained in its coefficient                               *)
(* Created:  11 June, 2007 by DP                                        *)
(* Code is in File:  newhomotopyoper/hodegopr.m                         *)
(* Last Modified: 24 March, 2008 17:28 by DP at CSM                     *)
(************************************************************************)

degreeOperator[term_, numDependentVar_Integer, IndepVars_List] :=
Module[{numIndVar, termsWithDerList, depVarList, combinedList,
    partialWRTpartialList, counter = 1, finalDegree = {}},
    
  printDO["At DO IN 2, the term given is:"];
  printDO[term];
  printDO["At DO IN 2, the number of dependent variables is ",
      numDependentVar];
  printDO["At DO IN 2, the list of independent variables contains ",
      IndepVars];

  numIndVar = Length[IndepVars];

  (* Add results for all dependent variables.                           *)
  While[counter <= numDependentVar,

    (* Create a list of patial derivatives existing in the expression.  *)
    termsWithDerList = Cases[term,
        Derivative[__][u[counter]][__], {0, Infinity}];
    printDO["At DO 2, all derivatives in the expression are given in"<>
        " the following list:"];
    printDO[termsWithDerList];

    (* Check for terms that have no derivatives.                        *)
    depVarList = Cases[term, u[counter][__], {0, Infinity}];

    (* Put depVarList and derivative lists together. *)
    combinedList = Union[depVarList, termsWithDerList];
    printDO["At DO 3, derivative terms and non derivative terms"<>
        " contained in the given expression are shown in the following list:"];
    printDO[combinedList];

    (* Create a list of the partial derivatives of the expression with  *)
    (* respect to the partial derivative list created above. *)
    partialWRTpartialList = Map[D[term, #] &, combinedList];
    printDO["At DO 4, the partial derivatives of the expression with"<>
        " respect to the terms given in the term list above are given in"<>
        " the following list:"];
    printDO[partialWRTpartialList];

    (* Multiply the partial term to the partial derivative just found.  *)
    partialWRTpartialList = Thread[Times[partialWRTpartialList, combinedList]];
    partialWRTpartialList = Expand[partialWRTpartialList];
    printDO["At DO 5, the partial derivatives multiplied to the partial"<>
        " expression:"];
    printDO[partialWRTpartialList];

    finalDegree = Flatten[Append[finalDegree, partialWRTpartialList]];
    printDO["At DO 6, a list of all expressions computed so far:"];
    printDO[finalDegree];
    counter ++
  ]; (* end While *)

  finalDegree = Factor[Apply[Plus, finalDegree]];

  Clear[numIndVar, termsWithDerList, depVarList, combinedList,
    partialWRTpartialList, counter];

  Return[finalDegree]
]; (* end Module degreeOperator *)

degreeOperExpression[expression_, numDepVars_Integer, indVarsList_List] :=
Module[{expandedExpression, degreeOperResult},

  If[debugDegreeOperator, printDO = Print, Clear[printDO], Clear[printDO]];

  printDO["debug code: DO, Function: degreeOperator, File: hodegopr.m"];

  printDO["At DO IN 1, the given expression is:"];
  printDO[expression];
  printDO["At DO IN 1, the number of dependent variables is ", numDepVars];
  printDO["At DO IN 1, the list of independent variables contains ",
      indVarsList];

  expandedExpression = Expand[expression];

  expandedExpression = If[Head[expandedExpression] === Plus,
      Apply[List, expandedExpression], {expression}, {expression}];
  printDO["At DO 1, a list of terms in the expression:"];
  printDO[expandedExpression];

  degreeOperResult = Map[degreeOperator[#, numDepVars, indVarsList] &,
      expandedExpression];
  printDO["At DO 7, the expression returned by the degree operator:"];
  printDO[degreeOperResult];

  Return[{expandedExpression, degreeOperResult}]
]; (* end Module degreeOperExpression *)

(* ##########       Function: rationalExpressionCheck        ########## *)

(************************************************************************)
(* rationalExpressionCheck[expressionGiven, indepVars, numDepVars]      *)
(* Purpose:  To check an expression for terms that fall in the kernel   *)
(*           of the degree operator, then to separate those terms from  *)
(*           the rest of the expression.  The separate pieces will be   *)
(*           run through the homotopy operator independently            *)
(* Input:    An expression                                              *)
(*           A list of terms in the denominator                         *)
(* Output:   The combined result from the homotopy operator             *)
(* Created:  27 March, 2008 by DP                                       *)
(* Code is in File:  newhomotopyoper/horatchk.m                         *)
(* Last Modified: 1 April, 2008, 18:03 by DP at CSM                     *)
(************************************************************************)

rationalExpressionCheck[expressionGiven_, indepVars_, numDepVars_] :=
Module[{expnDenom, revisedExpn, degree, zeroDegreeTerms, shiftedTerms,
    afterHomotopy, degreeChk, printREC},

  If[debugRationalExpressionCheck, printREC = Print, Clear[printREC],
      Clear[printREC]];
  printREC["debug code: REC, Function: rationalExpressionCheck,"<>
      " File: horatchk.m"];

  printREC["At REC IN, the given expression is:"];
  printREC[expressionGiven];
  printREC["At REC IN, the list of independent variables is ", indepVars];
  printREC["At REC IN, the number of dependent variables is ", numDepVars];

  (* Check the expression for an overall degree of zero.  If the degree *)
  (* is zero, the homotopy operator will not work until changes are     *)
  (* made.                                                              *)
  expnDenom = MemberQ[expressionGiven, Power[n1__, n2_] /;
      MemberQ[n1, u[_][__] | Derivative[__][u[_]][__], {0, Infinity}] &&
      n2 < 0, {0, Infinity}];
  printREC["At REC 1, True indicates that the expression is rational: ",
      expnDenom];

  If[expnDenom,
    {revisedExpn, degree} = degreeOperExpression[expressionGiven,
        numDepVars, indepVars];
    printREC["At REC 2, the results after applying the degree operator:"];
    printREC[degree];
    printREC["At REC 3, if terms in the expression have denominators, the"<>
        " result here is the expanded expression in a list.  If there are no"<>
        " denominators, the original expression is given."];
    printREC[revisedExpn];

    (* If the expression has a term with a degree of zero, shift one of *)
    (* the terms in the denominator of the degree zero term.            *)
    If[MemberQ[degree, 0],
      zeroDegreeTerms = Extract[revisedExpn, Position[degree, 0, {0,1}]];
      printREC["At REC 4, the terms in the expression that have degree 0:"];
      printREC[zeroDegreeTerms];

      shiftedTerms = shiftDenominatorTerm[Apply[Plus, zeroDegreeTerms],
          indepVars, numDepVars];
      printREC["At REC 5, an attempt to apply a shift to one term in the"<>
          " denominator:"];
      printREC[shiftedTerms];

      revisedExpn = {Apply[Plus, Complement[revisedExpn, zeroDegreeTerms]],
          shiftedTerms};
      printREC["At REC 6,the original expression has been split into two"<>
          " pieces, a piece not in the kernel of the degree operator"<>
          " unchanged and a piece in the kernel of the degree operator"<>
          " with denominators shifted:"];
      printREC[revisedExpn];

      (* Recheck the degree before moving on.                           *)
      degreeChk = degreeOperExpression[Part[revisedExpn, 2], numDepVars,
          indepVars];
      printREC["At REC 7, the results after applying the degree operator"<>
          " to the part in the kernel of the degree operator:"];
      printREC[Part[degreeChk, 2]];

      (* If degree sends 0 to homotopyFinalIntegration, Mathematica's   *)
      (* Integrate will be called for 1D functions.                     *)
      degreeChk = If[MemberQ[Part[degreeChk, 2], 0], 0, 1];
      printREC["At REC 8, the degree check now shows ", degreeChk, ".  If"<>
          " the degree check is 1, the homotopy operator will work."];

      revisedExpn = If[degreeChk === 0, {expressionGiven}, revisedExpn];
      printREC["At REC 9, the given expression is shown if the degree"<>
          " check shows 0; otherwise the expression in REC 6 is shown:"];
      printREC[revisedExpn];

      (* Run the homotopy operator independently on the terms not in    *)
      (* the kernel of the degree operator, the on the terms in the     *)
      (* kernel of the degree operator.                                 *)
      afterHomotopy = Map[homotopyFinalIntegration[
          #, numDepVars, indepVars, degreeChk] &, revisedExpn];
      printREC["At REC 10, the result returned from the homotopy operator:"];
      printREC[afterHomotopy];

      (* Undo any shifts applied to terms in the kernel of the degree   *)
      (* operator.                                                      *)
      afterHomotopy = afterHomotopy /. Log[shift] -> 0 /. shift -> 0;
      printREC["At REC 11, the results returned from the homotopy operator"<>
          " after all shifts have been removed:"];
      printREC[afterHomotopy];

      afterHomotopy = MapThread[Plus, afterHomotopy, 1];
      printREC["At REC 12, the combined results returned from the homotopy"<>
          " operator:"];
      printREC[afterHomotopy],

    (* else *)
      degreeChk = 1;
      afterHomotopy = homotopyFinalIntegration[expressionGiven, numDepVars,
          indepVars, degreeChk];
      printREC["At REC 13, the result returned from the homotopy operator:"];
      printREC[afterHomotopy]
    ], (* end If MemberQ[degree] *)

  (* else *)
    degreeChk = 1;
    afterHomotopy = homotopyFinalIntegration[expressionGiven, numDepVars,
        indepVars, degreeChk];
    printREC["At REC 14, the result returned from the homotopy operator:"];
    printREC[afterHomotopy]
  ]; (* end If expnDenom *)

  (* Clear all variables not being passed back.                         *)
  Clear[expnDenom, revisedExpn, degree, zeroDegreeTerms, shiftedTerms,
      degreeChk, printREC];

  Return[afterHomotopy]
]; (* end Module rationalExpressionCheck *)

(* ##########      Function: shiftDenominatorTerm        ########## *)

(************************************************************************)
(* shiftDenominatorTerm[expressionGiven, indepVars, numDepVars]         *)
(* Purpose:  To shift a term in the denominator of a rational expression*)
(*           away from the origin.  This alleviates the homogenous      *)
(*           degree problem.                                            *)
(* Input:    An expression                                              *)
(*           A list of terms in the denominator                         *)
(* Output:   A list of degrees corresponding to each term in the        *)
(*               expression                                             *)
(* Created:  19 March, 2008 by DP                                       *)
(* Code is in File:  newhomotopyoper/hoshift.m                          *)
(* Last Modified: 22 May, 2008, 10:29 by DP at CSM                      *)
(************************************************************************)

(* The function tuples replaces Mathematics's Tuples which is not       *)
(* available in versions before 5.1.                                    *)
tuples[listOflists_List] :=
Module[{end = Length[listOflists], tableOfParts, e, i1, j1, tupledList},
  tableOfParts = {Table[part[listOflists, j1, e[j1]], {j1, 1, end}],
      Apply[Sequence, Table[{e[i1], 1, Length[Part[listOflists, i1]]},
      {i1, 1, end}]]};
  tupledList = Flatten[Apply[Table, tableOfParts] /. part -> Part, end - 1];
  Return[tupledList]
]; (* end Module tuples *)

shiftDenominatorTerm[expressionGiven_, indepVars_, numDepVars_] :=
Module[{expression, zero, denominatorList, denominatorTermList, orderList,
    commonTerms, varDerResult = 1, termsToShift, shiftedExpression, pls,
    shiftedTerms, printSHD},

  If[debugShiftDenominatorTerm, printSHD = Print, Clear[printSHD],
      Clear[printSHD]];
  printSHD["debug code: SHD, Function: shiftDenominatorTerm,"<>
      " File: hoshift.m"];

  printSHD["At SHD IN, the given expression is:"];
  printSHD[expressionGiven];
  printSHD["At SHD IN, the list of independent variables is ", indepVars];
  printSHD["At SHD IN, the number of dependent variables is ", numDepVars];

  (* Initializations *)
  Clear[shift];
  zero = Table[0, {i1, 1, numDepVars}];
  expression = Expand[expressionGiven];
  printSHD["At SHD 1, the given expression after Expand has been applied:"];
  printSHD[expression];

  (* Get the denominators of the terms with a degree of zero.           *)
  denominatorList = Map[Denominator, Apply[List, expression]];
  printSHD["At SHD 2, a list of all denominators in the expression:"];
  printSHD[denominatorList];

  denominatorTermList = Union[Map[Cases[#, u[_][__] |
      Derivative[__][u[_]][__], {0, Infinity}] &, denominatorList]];
  denominatorTermList = DeleteCases[denominatorTermList, {}];
  printSHD["At SHD 3, jet space terms found in the denominators of all"<>
      " degree zero terms:"];
  printSHD[denominatorTermList];

  (* Sort the terms in the denominator by order.                        *)
  orderList = denominatorTermList /. {u[_][__] -> 0,
      Derivative[n1__][u[_]][__] -> pls[n1]} /. pls -> Plus;
  printSHD["At SHD 4, the orders corresponding to the jet space terms in"<>
       " the denominator list:"];
  printSHD[orderList];

  denominatorTermList = Sort[Thread[List[orderList, denominatorTermList]]];
  printSHD["At SHD 5, the denominators paired with their order sorted in"<>
      " decending order:"];
  printSHD[denominatorTermList];

  denominatorTermList = Map[Part[#, 2] &, denominatorTermList];
  printSHD["At SHD 6, the denominator list sorted in by order:"];
  printSHD[denominatorTermList];

  commonTerms = Apply[Intersection, denominatorTermList];
  printSHD["At SHD 7, all denominators contain these terms:"];
  printSHD[commonTerms];

  denominatorTermList = DeleteCases[denominatorTermList,
      Apply[Alternatives, commonTerms], {0, 3}];
  printSHD["At SHD 8, list of terms in the denominator after all common"<>
      " terms have been removed:"];
  printSHD[denominatorTermList];

  While[commonTerms =!= {} && varDerResult =!= zero ,
    termsToShift = First[commonTerms];
    printSHD["At SHD 9, this term will be shifted:"];
    printSHD[termsToShift];

    shiftedTerms = Rule[termsToShift, termsToShift + shift];
    printSHD["At SHD 10, the rule for making the shift on the term:"];
    printSHD[shiftedTerms];

    shiftedExpression = expression /. shiftedTerms;
    printSHD["At SHD 11, the expanded given expression with the shift:"];
    printSHD[shiftedExpression];

    varDerResult = Factor[Table[variationalDerivativeMultiD[
        shiftedExpression, i1, indepVars], {i1, 1, numDepVars}]];
    printSHD["At SHD 12, the variational derivative applied to the"<>
        " shifted expression:"];
    printSHD[varDerResult];

    commonTerms = If[commonTerms =!= {},
        Drop[commonTerms, 1], commonTerms];
  ]; (* end While *)

  denominatorTermList = If[varDerResult === zero,
      denominatorTermList, tuples[denominatorTermList]];
  printSHD["At SHD 13, all groups of denominator terms that can be"<>
      " shifted:"];
  printSHD[denominatorTermList];

  While[Flatten[denominatorTermList] =!={} && varDerResult =!= zero,
    termsToShift = First[denominatorTermList];
    printSHD["At SHD 14, the terms in the following list will be shifted:"];
    printSHD[termsToShift];

    shiftedTerms = Map[Rule[#, # + shift] &, termsToShift];
    printSHD["At SHD 15, the rules for making the shift on the terms:"];
    printSHD[shiftedTerms];

    shiftedExpression = expression /. shiftedTerms;
    printSHD["At SHD 16, the expanded given expression:"];
    printSHD[shiftedExpression];

    varDerResult = Factor[Expand[Table[variationalDerivativeMultiD[
        shiftedExpression, i1, indepVars], {i1, 1, numDepVars}]]];
    printSHD["At SHD 17, the variational derivative applied to the"<>
        " shifted expression:"];
    printSHD[varDerResult];

    denominatorTermList = If[denominatorTermList =!= {},
        Drop[denominatorTermList, 1], denominatorTermList];
  ]; (* end While *)

  shiftedExpression = If[varDerResult =!= zero,
      expression, shiftedExpression];

  (* Clear all variables not being passed back.                         *)
  Clear[expression, zero, denominatorList, denominatorTermList, orderList,
    commonTerms, varDerResult, termsToShift, shiftedTerms, printSHD];

  Return[shiftedExpression]
]; (* end Module shiftDenominatorTerm *)

(* ##########      Function: homotopyFinalIntegration        ########## *)

(************************************************************************)
(* homotopyFinalIntegration[expression, noOfDependentVar, indepVarList, *)
(*     degCheck]                                                        *)
(* Purpose:  To integrate an expression involving partial derivatives   *)
(*           in 1D and to invert the divergence operator in 2D and 3D.  *)
(* Input:    An expression                                              *)
(*           The number of dependent variables                          *)
(*           A list of independent variables                            *)
(* Output:   The inverted expression.                                   *)
(* Created:  6 March, 2006 by DP                                        *)
(* Code is in File:  newhomotopyoper/hofinint.m                         *)
(* Last Modified: 7 June, 2009, 16:05 by DP at CSM                      *)
(************************************************************************)

homotopyFinalIntegration[expression_, noOfDependentVar_, indepVarList_,
    degCheck_Integer] :=
Module[{numIndVars = Length[indepVarList], zero, lambdaPAR, expression2,
    indVar, uReplaceuLambdaRule, integrandList, sumOfDepComponents, lambdaSum,
    lambdaResult, homotopyAtOne, homotopyAtZero, interimResult, constantsk,
    divInterimResult, termsWithConstants, termsWOConstants, coefSystem,
    coefSystemSoln, homotopyResult, checkResult, dd, m1, m2, m3, m4, m5,
    printHFI},

  If[debugHomotopyFinalIntegration, printHFI = Print, Clear[printHFI],
      Clear[printHFI]];
  printHFI["debug code: HFI, Function: homotopyFinalIntegration,"<>
      " File: hofinint.m"];

  printHFI["At HFI IN, the given expression is:"];
  printHFI[expression];
  printHFI["At HFI IN, the number of dependent variables is ",
      noOfDependentVar];
  printHFI["At HFI IN, the list of independent variables given is ",
      indepVarList];
  printHFI["At HFI IN, the degree check is ", degCheck];

  (* If the expression is zero, return 0 vector.                        *)
  zero = Table[0, {i1, 1, numIndVars}];
  printHFI["At HFI 1, zero = ", zero];
  If[expression === 0, Return[zero]];

  (* Check to see if lambda is a parameter in the expression.           *)
  (* If so, lambda to another parameter, lambdaPAR.                     *)
  expression2 = If[MemberQ[expression, lambda, {0, Infinity}],
      expression /. lambda -> lambdaPAR, expression, expression];
  printHFI["At HFI 2, the expression after a check for lambdas:"];
  printHFI[expression2];

  indVar[1] = Part[indepVarList, 1];
  If[numIndVars > 1, indVar[2] = Part[indepVarList, 2]];
  If[numIndVars > 2, indVar[3] = Part[indepVarList, 3]];

  (* The uReplaceuLambda rule takes u from an expression and replaces   *)
  (* it with lambda*u.                                                  *)
  uReplaceuLambdaRule = {u[m1_][m2__] -> lambda*u[m1][m2],
      Derivative[m3__][u[m4_]][m5__] -> lambda*Derivative[m3][u[m4]][m5]};

  If[degCheck =!= 0,
    (* Use the integrand functions to determine the integrands for the  *)
    (* homotopy operator.                                               *)
    integrandList = Switch[numIndVars,
      1, Table[{integrandLeftDerList1D[expression2, j1, indepVarList]},
          {j1, 1, noOfDependentVar}],
      2, Table[integrandLeftDerList2D[expression2, j1, indepVarList],
          {j1, 1, noOfDependentVar}],
      3, Table[integrandLeftDerList3D[expression2, j1, indepVarList],
          {j1, 1, noOfDependentVar}]];
    printHFI["At HFI 3, a list of pieces produced by the integrand operators"<>
        " is given as follows.  The innermost nested list gives the"<>
        " integrand pieces for each dependent variable.  The outermost"<>
        " list{s} are groupings over each independent variable."];
    printHFI[integrandList];

    If[globalHomotopyVerbose && TimeUsed[] - timeKeeper > 5.0*numIndVars,
      printOUT["Beginning to sum the expressions determined for each"<>
          " dependent variable."]
    ]; (* end If globalHomotopyVerbose *)

    (* Add the terms produced by the integrand operator produced by     *)
    (* each dependent variable.                                         *)
    sumOfDepComponents = Apply[Plus, integrandList];
    printHFI["At HFI 4, the previous list has been refined by adding"<>
        " all of the dependent parts together.  The result consists of one"<>
        " expression for each independent variable."];
    printHFI[sumOfDepComponents];

    If[globalHomotopyVerbose,
      printOUT["Proceding to the integration step."];
      If[TimeUsed[] - timeKeeper > 6.0*numIndVars,
        printOUT["All jet space coordinates are being replaced by "<>
            ToString[\[Lambda]] <>"*coordinate."]
      ] (* end If TimeUsed *)
    ]; (* end If globalHomotopyVerbose *)

    (* Apply the lambda replace rule.                                   *)
    lambdaSum = sumOfDepComponents /. uReplaceuLambdaRule;
    printHFI["At HFI 5, the result when u is replaced with u*lambda:"];
    printHFI[lambdaSum];

    (* Divide all terms by lambda and integrate with respect to lambda. *)
    lambdaSum = Expand[lambdaSum/lambda];
    printHFI["At HFI 6, the result after all terms have been divided by"<>
        " lambda:"];
    printHFI[lambdaSum];

    If[globalHomotopyVerbose && TimeUsed[] - timeKeeper > 6.0*numIndVars,
      printOUT["Beginning integration over "<> ToString[\[Lambda]] <>"."]
    ]; (* end If globalHomotopyVerbose *)

    lambdaResult = Map[Integrate[#, lambda] &, lambdaSum, {1}];
    lambdaResult = Factor[Expand[lambdaResult]];
    printHFI["At HFI 7, the result after integrating with respect to lambda:"];
    printHFI[lambdaResult],

  (* else *)
    If[numIndVars === 1,
      lambdaResult = {Integrate[expression, indVar[1]]};
      printOUT["Mathematica's Integrate function has been called"<>
          " because a zero degree term existed in the expression to be"<>
          " integrated."];
      Return[lambdaResult],
    (* else *)
      printOUT["The homotopy operator is unable to invert the divergence"<>
          " on the given expression because at least one term in the"<>
          " expression has an overall degree of zero."];
      Return[{}]
    ]; (* end If numIndVars *)
  ]; (* end If degCheck *)

  If[MemberQ[lambdaResult, _Integrate, {0, Infinity}],
    printOUT["The homotopy operator was unable to perform the auxillary"<>
        " integration required in the final step because the function is"<>
        " to complicated to integrate.  It is unable to provide a result."];
    lambdaResult = {}
  ]; (* end If MemberQ *)

  If[globalHomotopyVerbose,
    If[TimeUsed[] - timeKeeper > 6.0*numIndVars,
      printOUT["Evaluating the result from "<> ToString[\[Lambda]] <>" = 0"<>
          " to "<> ToString[\[Lambda]] <>" = 1."];
    ]; (* end If TimeUsed *)
    printOUT["------------------------------------------"]
  ]; (* end If globalHomotopyVerbose *)

  homotopyAtOne = If[MemberQ[lambdaResult, Sin[_] | Cos[_] | Tan[_] | Cot[_] |
      Sec[_] | Csc[_] | Sinh | Cosh[_] | Tanh[_] | Coth[_] | Sech[_] | Csc[_],
      {0, Infinity}] === {},
    Expand[TrigExpand[lambdaResult /. {lambda -> 1}]],
    If[MemberQ[lambdaResult, Log[_], {0, Infinity}],
      Expand[PowerExpand[lambdaResult /. {lambda -> 1}]],
      Expand[lambdaResult /. {lambda -> 1}]]];
  printHFI["At HFI 9, the result from the homotopy operator after"<>
      " lambda -> 1:"];
  printHFI[homotopyAtOne];

  (* Check the result after lambda -> 1 by comparing its divergence     *)
  (* with the given expression.  If they match, return and finish.      *)
  checkResult = divergenceCheck[expression2, homotopyAtOne, indepVarList];
  printHFI["At HFI 10, a check on the integration result over lambda when"
      " lambda -> 1:"];
  printHFI[checkResult];
  If[checkResult === 0,
    (* Clear all variables not being passed back.                       *)
    Clear[numIndVars, zero, expression2, indVar, uReplaceuLambdaRule, 
          lambdaSum, integrandList, sumOfDepComponents, lambdaResult, 
          checkResult, printHFI];
    Return[homotopyAtOne /. lambdaPAR -> lambda]
  ]; (* end If checkResult] *)

  Off[Power::infy];
  Off[\[Infinity]::indet];

  homotopyAtZero = If[MemberQ[lambdaResult, lambda, {0, Infinity}],
       Expand[lambdaResult /. lambda -> 0], zero];
  printHFI["At HFI 11, the integration result over lambda when lambda -> 0:"];
  printHFI[homotopyAtZero];

  (* Mathematica may produce "constants" in terms of the dependent      *)
  (* variables when it integrates the integrand expression over lambda. *)
  (* These extra terms must be removed before reverting from lambda     *)
  (* space back to the independent variable space.                      *)

  (* First adjust the expression for any lambdas in the denominator.    *)
  choices = ComplexInfinity | Indeterminate | Infinity;
  If[MemberQ[homotopyAtZero, choices, {0, Infinity}],
    numHomotopyAtZero = Expand[Map[Numerator, lambdaResult]];
    printHFI["At HFI 12, the numerators of the result after integration"<>
        " over lambda:"];
    printHFI[numHomotopyAtZero];

    headNumParts = Map[Head, numHomotopyAtZero];
    printHFI["At HFI 13, the heads of the components of the result after"<>
        " integration:"];
    printHFI[headNumParts];

    numHomotopyAtZero = Map[Level[#, {1}] &, numHomotopyAtZero];
    printHFI["At HFI 14, the numerators as a list of terms:"];
    printHFI[numHomotopyAtZero];

    denHomotopyAtZero = Map[Denominator, lambdaResult];
    printHFI["At HFI 15, a list of denominators of the result after"<> 
              " integration over lambda:"];
    printHFI[denHomotopyAtZero];

    splitHomotopyAtZero = Thread[tt[numHomotopyAtZero, 1/denHomotopyAtZero]] /.
        tt -> Times;
    printHFI["At HFI 16, a list of terms of the result after integration"<>
        " over lambda:"];
    printHFI[splitHomotopyAtZero];

    splitHomotopyAtZero = splitHomotopyAtZero /. lambda -> 0;
    printHFI["At HFI 17, the list of terms of the result after integration"<>
        " over lambda after lambda has been replaced with zero:"];
    printHFI[splitHomotopyAtZero];

    splitHomotopyAtZero = splitHomotopyAtZero /. {ComplexInfinity -> 0,
        Indeterminate -> 0};
    printHFI["At HFI 18, the list of terms of the result after integration"<>
        " over lambda after ComplexInfinity has been removed:"];
    printHFI[splitHomotopyAtZero];

    splitHomotopyAtZero = splitHomotopyAtZero /.
        Times[n1_, __] /; n1 === Infinity || n1 === -Infinity -> 0;
    printHFI["At HFI 19, the list of terms of the result after integration"<>
        " over lambda after +/-Infinity has been removed:"];
    printHFI[splitHomotopyAtZero];

    homotopyAtZero = Thread[app[headNumParts, splitHomotopyAtZero]] /.
        app -> Apply;
    printHFI["At HFI 20, the reconstructed result for lambda = 0:"];
    printHFI[homotopyAtZero]
  ]; (* end If MemberQ *)

  (* Remove numerical constants.                                         *)
  homotopyAtZero = If[MemberQ[homotopyAtZero, u[_][__] |
     Derivative[__][u[_]][__], {0, Infinity}], homotopyAtZero, zero];
  printHFI["At HFI 22, the integration result over lambda when lambda -> 0"<>
      " is changed to zero if the result contains only numerical constants:"];
  printHFI[homotopyAtZero];

  (* Second, check the expression with lambda -> 1 minus the expression *)
  (* with lambda -> 0.  If the divergence on the result matches the     *)
  interimResult = Expand[PowerExpand[Factor[homotopyAtOne - homotopyAtZero]]];
  printHFI["At HFI 23, the integration result with lambda -> 1 minus the"<>
      " integration result with lambda -> 0 before the shift is undone:"];
  printHFI[interimResult];

  checkResult = divergenceCheck[expression2, interimResult, indepVarList];
  printHFI["At HFI 24, a check on the integration result on the difference"<>
      " where lambda -> 1, then lambda -> 0 have been applied:"];
  printHFI[checkResult];
  If[checkResult === 0,
    On[Power::infy];
    On[\[Infinity]::indet];

    (* Clear all variables not being passed back.                       *)
    Clear[numIndVars, zero, expression2, indVar, uReplaceuLambdaRule, 
        lambdaSum, integrandList, sumOfDepComponents, lambdaResult, 
        homotopyAtOne, homotopyAtZero, checkResult, printHFI];

    Return[interimResult /. lambdaPAR -> lambda]
  ]; (* end If checkResult *)

  (* Third, it may be necessary to evaluate terms produced when         *)
  (* lambda -> 0.  Not all terms in this case can removed, so           *)
  (* comparisons must be made with the given expression.                *)
  If[homotopyAtZero =!= zero,
    homotopyAtZero = applyConstantsToVector[homotopyAtZero, numIndVars];
    constantsk = Part[homotopyAtZero, 4];
    homotopyAtZero = If[numIndVars === 1, Part[homotopyAtZero, 1],
        Table[Part[homotopyAtZero, j], {j, 1, numIndVars}]];
    printHFI["At HFI 25, zero terms with constants:"];
    printHFI[homotopyAtZero];

    interimResult = Expand[homotopyAtOne - homotopyAtZero];
    printHFI["At HFI 26, zero terms with constants subtracted from the"<>
        " result where lambda -> 1:"];
    printHFI[interimResult];

    divInterimResult = If[numIndVars === 1,
        D[Part[interimResult, 1], indVar[1]],
        If[numIndVars === 2,
            divergence[Flatten[{interimResult, 0}], indVar[1], indVar[2], dd],
        divergence[interimResult, indVar[1], indVar[2], indVar[3]]]];
    printHFI["At HFI 27, the revised homotopy result with constants on"<>
        " zero terms after the divergence has been applied:"];
    printHFI[divInterimResult];
    (* Compare the derivative/divergence of the result with the given   *)
    (* expression.                                                      *)
    divInterimResult = Factor[divInterimResult - expression];
    printHFI["At HFI 28, the difference between the revised homotopy"<>
        " result with constants and the original expression:"];
    printHFI[divInterimResult];

    divInterimResult = Expand[Numerator[divInterimResult]];
    printHFI["At HFI 29, the numerator of the difference:"];
    printHFI[divInterimResult];

    termsWithConstants = Apply[Plus,
        Cases[divInterimResult, Times[k[_], ___], {0, Infinity}]];
    printHFI["At HFI 30, the terms in the difference containing k[i]:"];
    printHFI[termsWithConstants];

    termsWOConstants = -Complement[divInterimResult, termsWithConstants];
    printHFI["At HFI 31, the terms in the difference without k[i]:"];
    printHFI[termsWOConstants];

    coefSystem = makeCoefficientEquations[termsWithConstants,
        constantsk, termsWOConstants, indepVarList, "caseTWO"];
    printHFI["At HFI 32, a system of equations formed from the coefficients:"];
    printHFI[coefSystem];

    coefSystemSoln = Flatten[Solve[coefSystem, constantsk]];
    printHFI["At HFI 33, the solution to the system of equations formed"<>
        " from the coefficients:"];
    printHFI[coefSystemSoln];

    homotopyResult = interimResult /. coefSystemSoln /. k[_] -> 0;
    homotopyResult = Factor[homotopyResult];
    printHFI["At HFI 34, the revised homotopy result (after removing"<>
        " constants left from Mathematica's integration over lambda):"];
    printHFI[homotopyResult];
  ]; (* end If homotopyAtZero *)

  On[Power::infy];
  On[\[Infinity]::indet];

  (* Clear all variables not being passed back.                         *)
  Clear[numIndVars, zero, expression2, indVar, uReplaceuLambdaRule, lambdaSum,
      integrandList, sumOfDepComponents, lambdaResult, homotopyAtOne,
      homotopyAtZero, interimResult, constantsk, divInterimResult, printHFI,
      termsWithConstants, termsWOConstants, coefSystem, coefSystemSoln];

  Return[homotopyResult /. lambdaPAR -> lambda]
]; (* end Module homotopyFinalIntegration *)

(* ##########        Function: integrandLeftDerList1D        ########## *)

(************************************************************************)
(* integrandLeftDerList1D[expression, dependentVar, indepVarList]       *)
(* Purpose:  To find the integrand for the 1-D homotopy operator        *)
(*           using the form where the derivatives of the dependent      *)
(*           varialble have been propagated outside the inner most sum  *)
(*           (after Kruskal et. al.) In this version, operators are     *)
(*           determined by manipulating  lists                          *)
(* Input:    An expression                                              *)
(*           A number for the dependent variable                        *)
(*           A list of independent variables                            *)
(* Output:   The result of applying the higher Euler operators to the   *)
(*           expression.                                                *)
(* Created:  22 February, 2006 by DP                                    *)
(* Code is in File:  newhomotopyoper/hoingr1d.m                         *)
(* Last Modified: 9 April, 2008, 12::40 by DP at CSM                    *)
(************************************************************************)

integrandLeftDerList1D[expression_, dependentVar_, indepVarList_] := Module[
    {derivativeList, orderList, i, maxi, partialWRTpartialList, innerSumList,
    dummyderiv, innerSum, outerSumList, integrandResult, k,
    indVar1 = Part[indepVarList, 1], printIL1D, printILLP},

  If[debugIntegrandLeftDerList, printIL1D = Print, Clear[printIL1D],
      Clear[printIL1D]];
  If[debugIntegrandLoop, printILLP = Print, Clear[printILLP],
      Clear[printILLP]];
  printIL1D["debug code: IL1D, Function: integrandLeftDerList1D,"<>
      " File: integrandLeftDerList1D.m"];

  printIL1D["Computing the integrand for the 1-D homotopy operator."];
  printIL1D["At IL1D IN, the expression given is:"];
  printIL1D[expression];
  printIL1D["At IL1D IN, the current dependent variable is"<>
      " u[", dependentVar, "]."];
  printIL1D["At IL1D IN, the independent variable is ", indepVarList];

  timeKeeper = TimeUsed[];

  If[globalHomotopyVerbose,
    printOUT["Beginning the computation of the homotopy integrand for"<>
        " dependent variable "<> ToString[dependentVar] <>".  The algorithm"<>
        " first finds partial derivatives with respect to the jet space"<>
        " coordinates found in the expression, then performs a series of"<>
        " total derivatives on the partial derivatives.  The results are"<>
        " then combined into a sum."]
  ]; (* end If globalHomotopyVerbose *)

  (* Create a list of patial derivatives existing in the expression.    *)
  derivativeList = Union[Cases[expression,
      Derivative[___][u[dependentVar]][indVar1], {0, Infinity}]];
  printIL1D["At IL1D 1, derivatives contained in the given expression are"<>
      " shown in the following list:"];
  printIL1D[derivativeList];

  (* Create a list of the orders of all derivatives that exist in the   *)
  (* expression.                                                        *)
  orderList = Flatten[Map[Head, Map[Head, derivativeList]] /.
      Derivative -> List];
  printIL1D["At IL1D 2, a list of orders for all derivatives is:"];
  printIL1D[orderList];

  maxi = Max[orderList, 0];
  printIL1D["At IL1D 3, the order for the expression is ", maxi];

  (* Create a list of the partial derivatives of the expression with    *)
  (* respect to the partial derivative list created above.              *)
  partialWRTpartialList = Map[D[expression, #] &, derivativeList];
  printIL1D["At IL1D 4, the partial derivatives of the expression with"<>
      " respect to the terms given in the derivative list above are given"<>
      " in the following list:"];
  printIL1D[partialWRTpartialList];

  (* Create a list of terms for the inner sum of the formula using      *)
  (* the partial derivative list created in the last step.              *)
  innerSumList = Table[Map[(-1)^(# - i - 1) &,
      Select[orderList, # > i &]]*
      Thread[dummyderiv[Extract[partialWRTpartialList,
      Position[orderList, _?(# > i &)]],
      Map[{indVar1, # - i - 1} &, Select[orderList, # > i &]]]]
      /. dummyderiv -> D, {i, 0, maxi - 1}];
  printIL1D["At IL1D 5, a list of inner sum terms for"<>
      " variable ", indVar1, " is:"];
  printIL1D[innerSumList];

  If[globalHomotopyVerbose && TimeUsed[] - timeKeeper > 2.0,
    printOUT["Beginning to sum all total derivatives."]
  ]; (* end If globalHomotopyVerbose *)

  (* Add the terms in the sum list.                                     *)
  innerSum = Expand[Apply[Plus, innerSumList, 1]];
  printIL1D["At IL1D 6, after adding all inner sum terms for"<>
      " variable ", indVar1, ":"];
  printIL1D[innerSum];

  (* Multiply the sum by the derivatives of the dependent varialble     *)
  (* with respect to x according to the outer sum.  The result is a     *)
  (* list of terms.                                                     *)
  outerSumList = Expand[Table[Part[innerSum, i + 1]*
      Derivative[i][u[dependentVar]][indVar1], {i, 0, maxi - 1}]];
  printIL1D["At IL1D 7, a list of outer sum terms for"<>
      " variable ", indVar1, " is:"];
  printIL1D[outerSumList];

  (* Add the terms from the list in the last step to get the result of  *)
  (* applying the higher Euler operator.                                *)
  integrandResult = Expand[Apply[Plus, outerSumList]];
  printIL1D["At IL1D 8, the integrand for the 1-D homotopy operator is as"<>
      " follows:"];
  printIL1D[integrandResult];

  If[globalHomotopyVerbose && TimeUsed[] - timeKeeper > 3.0,
    printOUT["All sums for dependent variable "<> ToString[dependentVar] <>
        " have been compiled."]
  ]; (* end If globalHomotopyVerbose *)

  (* Clear all variables not being passed back.                         *)
  Clear[derivativeList, orderList, i, k, maxi, partialWRTpartialList,
      innerSumList, dummyderiv, innerSum, outerSumList, indVar1, printIL1D,
      printILLP];

  Return[integrandResult]
]; (* end Module integrandLeftDerList1D *)

(* ##########        Function: integrandLeftDerList2D        ########## *)

(************************************************************************)
(* integrandLeftDerList2D[expression, dependentVar, indepVarList]       *)
(* Purpose:  To find the integrand for the 2-D homotopy operator        *)
(*           using the form where the derivatives of the dependent      *)
(*           varialble have been propagated outside the inner most sum  *)
(*           (after Kruskal et. al.) In this version, operators are     *)
(*           determined by manipulating  lists                          *)
(* Input:    An expression                                              *)
(*           A number for the dependent variable                        *)
(*           A list of independent variables                            *)
(* Output:   The result of applying the higher Euler operators to the   *)
(*           expression.                                                *)
(* Created:  27 February, 2006 by DP                                    *)
(* Code is in File:  newhomotopyoper/hoingr2d.m                         *)
(* Last Modified: 9 April, 2008, 12:42 by DP at CSM                     *)
(************************************************************************)

integrandLeftDerList2D[expression_, dependentVar_, indepVarList_] :=
Module[{derivativeList, ix, iy, kx, ky, dummyderiv, orderList,
    maxix, maxiy, partialWRTpartialList, innerSumListx, innerSumListy,
    innerSumx, innerSumy, positionNonZeroTermsx, positionNonZeroTermsy,
    listWOzerosx, listWOzerosy, indVar1 = Part[indepVarList, 1],
    indVar2 = Part[indepVarList, 2], integrandResultx, integrandResulty,
    outerSumListx, outerSumListy, printIL2D, printILLP},

  If[debugIntegrandLeftDerList, printIL2D = Print, Clear[printIL2D],
      Clear[printIL2D]];
  If[debugIntegrandLoop, printILLP = Print, Clear[printILLP],
      Clear[printILLP]];
  printIL2D["debug code: IL2D, Function: integrandLeftDerList2D,"<>
      " File: integrandLeftDerList2D.m"];

  printIL2D["Computing the integrand for the 2-D homotopy operator."];
  printIL2D["At IL2D IN, the expression given is:"];
  printIL2D[expression];
  printIL2D["At IL2D IN, the current dependent variable is"<>
      " u[", dependentVar, "]."];
  printIL2D["At IL2D IN, the independent variables are ", indepVarList];

  timeKeeper = TimeUsed[];

  If[globalHomotopyVerbose,
    printOUT["Beginning the computation of the homotopy integrand for"<>
        " dependent variable "<> ToString[dependentVar] <>".  The algorithm"<>
        " first finds partial derivatives with respect to the jet space"<>
        " coordinates found in the expression, then performs a series of"<>
        " total derivatives on the partial derivatives.  The results are"<>
        " then combined into a sum."]
  ]; (* end If globalHomotopyVerbose *)

  (* Create a list of patial derivatives existing in the expression.    *)
  derivativeList = Union[Cases[expression,
      Derivative[___][u[dependentVar]][indVar1, indVar2], {0, Infinity}]];
  printIL2D["At IL2D 1, the following list contains derivative terms from"<>
      " the given expression:"];
  printIL2D[derivativeList];

  (* Create lists of the orders of all derivatives that exist in the    *)
  (* expression. *)
  orderList = Map[Head, Map[Head, derivativeList]] /. Derivative -> List;
  printIL2D["At IL2D 2, a list of orders for the derivative terms shown"<>
      " above with orders for each variable listed in the form"<>
      " {", indVar1, ", ", indVar2, "} is:"];
  printIL2D[orderList];

  (* Determine the highest order derivative for each variable.          *)
  maxix = Max[Map[Part[#, 1] &, orderList], 0];
  maxiy = Max[Map[Part[#, 2] &, orderList], 0];
  printIL2D["At IL2D 3, the order for the expression in"<>
      " variable ", indVar1, " is ", maxix, " and the order for the"<>
      " expression in variable ", indVar2, " is ", maxiy, "."];

  (* Create lists for each variable of the partial derivatives of the   *)
  (* expression with respect to the partial derivative list created     *)
  (* above.                                                             *)
  partialWRTpartialList = Map[D[expression, #] &, derivativeList];
  printIL2D["At IL2D 4, the partial derivatives of the expression with"<>
      " respect to the terms given in the derivative list above are:"];
  printIL2D[partialWRTpartialList];

  (* Create lists of terms for the inner sum of the new formula using   *)
  (* the partial derivative list created in the last step.              *)
  innerSumListx = Table[Map[(-1)^(Part[#, 1] + Part[#, 2] - ix - iy - 1)*
      Binomial[Part[#, 1] + Part[#, 2] - ix - iy - 1, Part[#1, 1] - ix - 1]/
      Binomial[Part[#, 1] + Part[#, 2], Part[#, 1]] &,
      Select[Select[orderList, Part[#, 1] > ix &], Part[#, 2] >= iy &]]*
      Thread[dummyderiv[Extract[partialWRTpartialList,
      Position[orderList, {_, _}?(Part[#, 1] > ix && Part[#, 2] >= iy &)]],
      Map[{indVar1, Part[#, 1] - ix - 1} &,
      Select[Select[orderList, Part[#, 1] > ix &], Part[#, 2] >= iy &]],
      Map[{indVar2, Part[#, 2] - iy} &,
      Select[Select[orderList, Part[#, 1] > ix &], Part[#, 2] >= iy &]]]]
      /. dummyderiv -> D, {iy, 0, maxiy}, {ix, 0, maxix - 1}];
  printIL2D["At IL2D 5x, a list of inner sum terms for"<>
      " variable ", indVar1, " is:"];
  printIL2D[innerSumListx];

  innerSumListy = Table[Map[(-1)^(Part[#, 1] + Part[#, 2] - ix - iy - 1)*
      Binomial[Part[#, 1] + Part[#, 2] - ix - iy - 1, Part[#, 2] - iy - 1]/
      Binomial[Part[#, 1] + Part[#, 2], Part[#, 2]] &,
      Select[Select[orderList, Part[#, 1] >= ix &], Part[#, 2] > iy &]]*
      Thread[dummyderiv[Extract[partialWRTpartialList,
      Position[orderList, {_, _}?(Part[#, 1] >= ix && Part[#, 2] > iy &)]],
      Map[{indVar1, Part[#, 1] - ix} &,
      Select[Select[orderList, Part[#, 1] >= ix &], Part[#, 2] > iy &]],
      Map[{indVar2, Part[#, 2] - iy - 1} &,
      Select[Select[orderList, Part[#, 1] >= ix &], Part[#, 2] > iy &]]]]
      /. dummyderiv -> D, {iy, 0, maxiy - 1}, {ix, 0, maxix}];
  printIL2D["At IL2D 5y, a list of inner sum terms for"<>
      " variable ", indVar2, " is:"];
  printIL2D[innerSumListy];

  If[globalHomotopyVerbose && TimeUsed[] - timeKeeper > 2.0,
    printOUT["Beginning to sum all total derivatives."]
  ]; (* end If globalHomotopyVerbose *)

  (* Add the terms in the sum lists to produce lists of inner sums      *)
  (* corresponding to each outer index for ix and iy. *)
  innerSumx = Expand[Apply[Plus, innerSumListx, {2, 2}]];
  printIL2D["At IL2D 6x, after adding all inner sum terms for"<>
      " variable ", indVar1, ":"];
  printIL2D[innerSumx];

  innerSumy = Expand[Apply[Plus, innerSumListy, {2, 2}]];
  printIL2D["At IL2D 6y, after adding all inner sum terms for"<>
      " variable ", indVar2, ":"];
  printIL2D[innerSumy];

  (* Remove all zeros from the lists of sums above and keep track of    *)
  (* the positions of the remaining terms.  The positions give the      *)
  (* outer index, ix and iy, for each term (in reverse order).          *)
  positionNonZeroTermsx = Position[innerSumx,
      _?(# =!= 0 && # =!= List &), {2, 2}];
  printIL2D["At IL2D 7x, the positions of all nonzero terms in the inner"<>
      " sum list for ", indVar1, " shifted by {1, 1, 1}:"];
  printIL2D[positionNonZeroTermsx];

  positionNonZeroTermsy = Position[innerSumy,
      _?(# =!= 0 && # =!= List &), {2, 2}];
  printIL2D["At IL2D 7y, the positions of all nonzero terms in the inner"<>
      " sum list for ", indVar2, " shifted by {1, 1, 1}:"];
  printIL2D[positionNonZeroTermsy];

  listWOzerosx = Flatten[DeleteCases[innerSumx, 0, 2]];
  printIL2D["At IL2D 8x, the inner sum list for ", indVar1, " after"<>
      " removing zero terms:"];
  printIL2D[listWOzerosx];

  listWOzerosy = Flatten[DeleteCases[innerSumy, 0, 2]];
  printIL2D["At IL2D 8y, the inner sum list for " indVar2, " after"<>
      " removing zero terms:"];
  printIL2D[listWOzerosx];

  (* Multiply the sum by the partial derivative of the dependent        *)
  (* variable with respect to the outer sum indices (ix and iy) and     *)
  (* constant determined by the outer sum indicies.                     *)
  outerSumListx = Thread[Times[listWOzerosx,
      Map[Binomial[Part[#, 2] + Part[#, 1] - 2, Part[#, 2] - 1]*
      Derivative[Part[#, 2] - 1, 
                 Part[#, 1] - 1][u[dependentVar]][indVar1, indVar2] &,
      positionNonZeroTermsx]]];
  printIL2D["At IL2D 9x, a list of outer sum terms for"<>
      " variable ", indVar1, " is:"];
  printIL2D[outerSumListx];

  outerSumListy = Thread[Times[listWOzerosy,
      Map[Binomial[Part[#, 2] + Part[#, 1] - 2, Part[#, 1] - 1]*
      Derivative[Part[#, 2] - 1, 
                 Part[#, 1] - 1][u[dependentVar]][indVar1, indVar2] &,
      positionNonZeroTermsy]]];
  printIL2D["At IL2D 9y, a list of outer sum terms for"<>
      " variable ",indVar2, " is:"];
  printIL2D[outerSumListy];

  (* Add the terms from the lists in the last step to get the result of *)
  (* applying the higher Euler operator for each variable.              *)
  integrandResultx = Expand[Apply[Plus, outerSumListx]];
  integrandResulty = Expand[Apply[Plus, outerSumListy]];
  printIL2D["At IL2D 10, the integrands for the 2-D homotopy operator for"<>
      " dependent variable u[", dependentVar, "] are as follows:"];
  printIL2D["For variable ", indVar1, ": ", integrandResultx];
  printIL2D["For variable ", indVar2, ": ", integrandResulty];

  If[globalHomotopyVerbose && TimeUsed[] - timeKeeper > 3.0,
    printOUT["All sums for dependent variable "<> ToString[dependentVar] <>
        " have been compiled."]
  ]; (* end If globalHomotopyVerbose *)

  (* Clear all variables not being passed back.                         *)
  Clear[derivativeList, ix, iy, kx, ky, dummyderiv, orderList,
    maxix, maxiy, partialWRTpartialList, innerSumListx, innerSumListy,
    innerSumx, innerSumy, positionNonZeroTermsx, positionNonZeroTermsy,
    listWOzerosx, listWOzerosy, indVar1, indVar2, outerSumListx,
    outerSumListy, printIL2D, printILLP];

  Return[{integrandResultx, integrandResulty}]
]; (* end Module  integrandLeftDerList2D *)

(* ##########        Function: integrandLeftDerList3D        ########## *)

(************************************************************************)
(* integrandLeftDerList3D[expression, dependentVar, indepVarList]       *)
(* Purpose:  To find the integrand for the 3-D homotopy operator        *)
(*           using the form where the derivatives of the dependent      *)
(*           varialble have been propagated outside the inner most sum  *)
(*           (after Kruskal et. al.) In this version, operators are     *)
(*           determined by manipulating  lists                          *)
(* Input:    An expression                                              *)
(*           A number for the dependent variable                        *)
(*           A list of independent variables                            *)
(* Output:   The result of applying the higher Euler operators to the   *)
(*           expression.                                                *)
(* Created:  2 March, 2006 by DP                                        *)
(* Code is in File:  newhomotopyoper/hoingr3d.m                         *)
(* Last Modified: 9 April, 2008, 12:43 by DP at CSM                     *)
(************************************************************************)

integrandLeftDerList3D[expression_, dependentVar_, indepVarList_] :=
Module[{derivativeList, orderList, maxix, maxiy, maxiz, partialWRTpartialList,
    innerSumListx, innerSumListy, innerSumListz, ix, iy, iz, kx, ky, kz,
    innerSumx, innerSumy, innerSumz, positionNonZeroTermsx,
    positionNonZeroTermsy, positionNonZeroTermsz, listWOzerosx,
    listWOzerosy, listWOzerosz, dummyderiv, indVar1 = Part[indepVarList, 1],
    indVar2 = Part[indepVarList, 2], indVar3 = Part[indepVarList, 3],
    outerSumListx, outerSumListy, outerSumListz, integrandResultx,
    integrandResulty, integrandResultz, printIL3D, printILLP},

  If[debugIntegrandLeftDerList, printIL3D = Print, Clear[printIL3D],
      Clear[printIL3D]];
  If[debugIntegrandLoop, printILLP = Print, Clear[printILLP],
      Clear[printILLP]];
  printIL3D["debug code: IL3D, Function: integrandLeftDerList3D,"<>
      " File: integrandLeftDerList3D.m"];

  printIL3D["Computing the integrand for the 3-D homotopy operator."];
  printIL3D["At IL3D IN, the expression given is:"];
  printIL3D[expression];
  printIL3D["At IL3D IN, the current dependent variable is"<>
      " u[", dependentVar, "]."];
  printIL3D["At IL3D IN, the independent variables are ", indepVarList];

  timeKeeper = TimeUsed[];

  If[globalHomotopyVerbose,
    printOUT["Beginning the computation of the homotopy integrand for"<>
        " dependent variable "<> ToString[dependentVar] <>".  The algorithm"<>
        " first finds partial derivatives with respect to the jet space"<>
        " coordinates found in the expression, then performs a series of"<>
        " total derivatives on the partial derivatives.  The results are"<>
        " then combined into a sum."]
  ]; (* end If globalHomotopyVerbose *)

  (* Create a list of patial derivatives existing in the expression.    *)
  derivativeList = Union[Cases[expression,
      Derivative[___][u[dependentVar]][indVar1, indVar2, indVar3],
      {0, Infinity}]];
  printIL3D["At IL3D 1, the following list contains derivative terms from"<>
      " the given expression:"];
  printIL3D[derivativeList];

  (* Create lists of the orders of all derivatives that exist in the    *)
  (* expression.                                                        *)
  orderList = Map[Head, Map[Head, derivativeList]] /. Derivative -> List;
  printIL3D["At IL3D 2, a list of orders for the derivative terms shown"<>
      " above with orders for each variable listed in the form"<>
      " {", indVar1, ", ", indVar2, ", ", indVar3, "} is:"];
  printIL3D[orderList];

  (* Determine the highest order derivative for each variable.          *)
  maxix = Max[Map[Part[#, 1] &, orderList], 0];
  maxiy = Max[Map[Part[#, 2] &, orderList], 0];
  maxiz = Max[Map[Part[#, 3] &, orderList], 0];
  printIL3D["At IL3D 3, the order for the expression in"<>
      " variable ", indVar1, " is ", maxix, ",\nthe order for the"<>
      " expression in variable ", indVar2, " is ", maxiy, ",\nand the"<>
      " order for the expression in variable ", indVar3, " is ", maxiz, "."];

  (* Create lists for each variable of the partial derivatives of the   *)
  (* expression with respect to the partial derivative list created     *)
  (* above.                                                             *)
  partialWRTpartialList = Map[D[expression, #] &, derivativeList];
  printIL3D["At IL3D 4, the partial derivatives of the expression with"<>
      " respect to the terms given in the derivative list above are:"];
  printIL3D[partialWRTpartialList];

  (* Create lists of terms for the inner sum of the higher Euler        *)
  (* formula using the partial derivative list created in the last step.*)
  innerSumListx = Table[Map[(-1)^(Part[#, 1] + Part[#, 2] +
      Part[#, 3] - ix - iy - iz - 1)*
      Binomial[Part[#, 1] + Part[#, 2] + Part[#, 3] - ix - iy - iz - 1,
          Part[#, 1] - ix - 1]*
      Binomial[Part[#, 2] + Part[#, 3] - iy - iz,  Part[#, 2] - iy]/
      (Binomial[Part[#, 1] + Part[#, 2] + Part[#, 3], Part[#, 1]]*
      Binomial[Part[#, 2] + Part[#, 3], Part[#, 3]]) &,
      Select[Select[ Select[orderList,
          Part[#, 1] > ix &], Part[#, 2] >= iy &], Part[#, 3] >= iz &]]*
      Thread[dummyderiv[Extract[partialWRTpartialList,
      Position[orderList, {_, _, _}?(Part[#1, 1] > ix && Part[#1, 2] >= iy
          && Part[#, 3] >= iz &)]],
      Map[{indVar1, Part[#, 1] - ix - 1} &, Select[Select[Select[orderList,
          Part[#, 1] > ix &], Part[#, 2] >= iy &], Part[#, 3] >= iz &]],
      Map[{indVar2, Part[#, 2] - iy} &, Select[Select[Select[orderList,
          Part[#, 1] > ix &], Part[#, 2] >= iy &], Part[#, 3] >= iz &]],
      Map[{indVar3, Part[#, 3] - iz} &, Select[Select[Select[orderList,
          Part[#, 1] > ix &], Part[#, 2] >= iy &], Part[#, 3] >= iz &]]]] /.
      dummyderiv -> D, {iz, 0, maxiz}, {iy, 0, maxiy}, {ix, 0, maxix - 1}];
  printIL3D["At IL3D 5x, a list of terms for the inner sum for"<>
      " variable ", indVar1, " is:"];
  printIL3D[innerSumListx];

  innerSumListy = Table[Map[(-1)^(Part[#, 1] + Part[#, 2] +
      Part[#, 3] - ix - iy - iz - 1)*
      Binomial[Part[#, 1] + Part[#, 2] + Part[#, 3] - ix - iy - iz - 1,
          Part[#1, 2] - iy - 1]*
      Binomial[Part[#, 1] + Part[#, 3] - ix - iz,  Part[#, 3] - iz]/
      (Binomial[Part[#, 1] + Part[#, 2] + Part[#, 3], Part[#, 2]]*
      Binomial[Part[#, 1] + Part[#, 3], Part[#, 3]]) &,
      Select[Select[Select[orderList, Part[#, 1] >= ix &],
          Part[#, 2] > iy &], Part[#, 3] >= iz &]]*
      Thread[dummyderiv[Extract[partialWRTpartialList,
      Position[orderList, {_, _, _}?(Part[#, 1] >= ix && Part[#, 2] > iy
          && Part[#, 3] >= iz &)]],
      Map[{indVar1, Part[#, 1] - ix} &, Select[Select[Select[orderList,
          Part[#, 1] >= ix &], Part[#, 2] > iy &], Part[#, 3] >= iz &]],
      Map[{indVar2, Part[#, 2] - iy - 1} &, Select[Select[Select[orderList,
          Part[#, 1] >= ix &], Part[#, 2] > iy &], Part[#, 3] >= iz &]],
      Map[{indVar3, Part[#, 3] - iz} &, Select[Select[Select[orderList,
          Part[#, 1] >= ix &], Part[#, 2] > iy &], Part[#, 3] >= iz &]]]] /.
      dummyderiv -> D, {iz, 0, maxiz}, {iy, 0, maxiy - 1}, {ix, 0, maxix}];
  printIL3D["At IL3D 5y, a list of terms for the inner sum for"<>
      " variable ", indVar2, " is:"];
  printIL3D[innerSumListy];

  innerSumListz = Table[Map[(-1)^(Part[#, 1] + Part[#, 2] +
      Part[#, 3] - ix - iy - iz - 1)*
      Binomial[Part[#1, 1] + Part[#, 2] + Part[#, 3] - ix - iy - iz - 1,
          Part[#, 3] - iz - 1]*
      Binomial[Part[#, 1] + Part[#, 2] - ix - iy,  Part[#, 1] - ix]/
      (Binomial[Part[#, 1] + Part[#, 2] + Part[#, 3], Part[#, 3]]*
      Binomial[Part[#, 1] + Part[#, 2], Part[#, 1]]) &,
      Select[Select[Select[orderList, Part[#, 1] >= ix &],
          Part[#, 2] >= iy &], Part[#, 3] > iz &]]*
      Thread[dummyderiv[Extract[partialWRTpartialList,
      Position[orderList, {_, _, _}?(Part[#, 1] >= ix && Part[#, 2] >= iy
           && Part[#, 3] > iz &)]],
      Map[{indVar1, Part[#, 1] - ix} &, Select[Select[Select[orderList,
          Part[#, 1] >= ix &], Part[#, 2] >= iy &], Part[#, 3] > iz &]],
      Map[{indVar2, Part[#, 2] - iy} &, Select[Select[Select[orderList,
          Part[#, 1] >= ix &], Part[#, 2] >= iy &], Part[#, 3] > iz &]],
      Map[{indVar3, Part[#, 3] - iz - 1} &, Select[Select[Select[orderList,
          Part[#, 1] >= ix &], Part[#, 2] >= iy &], Part[#, 3] > iz &]]]] /.
      dummyderiv -> D, {iz, 0, maxiz - 1}, {iy, 0, maxiy}, {ix, 0, maxix}];
  printIL3D["At IL3D 5z, a list of terms for the inner sum for"<>
      " variable ", indVar3, " is:"];
  printIL3D[innerSumListz];

  If[globalHomotopyVerbose && TimeUsed[] - timeKeeper > 2.0,
    printOUT["Beginning to sum all total derivatives."]
  ]; (* end If globalHomotopyVerbose *)

  (* Add the terms in the sum lists to produce lists of inner sums      *)
  (* corresponding to each outer index for ix and iy.                   *)
  innerSumx = Expand[Apply[Plus, innerSumListx, {3, 3}]];
  printIL3D["At IL3D 6x, after adding all inner sum terms for"<>
      " variable ", indVar1, ":"];
  printIL3D[innerSumx];

  innerSumy = Expand[Apply[Plus, innerSumListy, {3, 3}]];
  printIL3D["At IL3D 6y, after adding all inner sum terms for"<>
      " variable ", indVar2, ":"];
  printIL3D[innerSumy];

  innerSumz = Expand[Apply[Plus, innerSumListz, {3, 3}]];
  printIL3D["At IL3D 6z, after adding all inner sum terms for"<>
      " variable ", indVar3, ":"];
  printIL3D[innerSumz];

  (* Remove all zeros from the lists of sums above and keep track of    *)
  (* the positions of the remaining terms.  The positions give the      *)
  (* outer index, ix, iy, and iz, for each term (The positions in the   *)
  (* list are given in the order z, y, x).                              *)
  positionNonZeroTermsx = Map[Part[#] - {1, 1, 1} &,
      Position[innerSumx, _?(# =!= 0 && # =!= List &), {3, 3}]];
  printIL3D["At IL3D 7x, the positions of all nonzero terms in the inner"<>
      " sum list for ", indVar1, " shifted by {1, 1, 1}:"];
  printIL3D[positionNonZeroTermsx];

  positionNonZeroTermsy = Map[Part[#] - {1, 1, 1} &,
      Position[innerSumy, _?(# =!= 0 && # =!= List &), {3, 3}]];
  printIL3D["At IL3D 7y, the positions of all nonzero terms in the inner sum"<>
      " sum list for ", indVar2, " shifted by {1, 1, 1}:"];
  printIL3D[positionNonZeroTermsy];

  positionNonZeroTermsz = Map[Part[#] - {1, 1, 1} &,
      Position[innerSumz, _?(# =!= 0 && # =!= List &), {3, 3}]];
  printIL3D["At IL3D 7z, the positions of all nonzero terms in the inner sum"<>
      " sum list for ", indVar3, " shifted by {1, 1, 1}:"];
  printIL3D[positionNonZeroTermsz];

  listWOzerosx = Flatten[DeleteCases[innerSumx, 0, 3]];
  printIL3D["At IL3D 8x, the inner sum list for ", indVar1, " after"<>
      " removing zero terms:"];
  printIL3D[listWOzerosx];

  listWOzerosy = Flatten[DeleteCases[innerSumy, 0, 3]];
  printIL3D["At IL3D 8y, the inner sum list for ", indVar2, " after"<>
      " removing zero terms:"];
  printIL3D[listWOzerosy];

  listWOzerosz = Flatten[DeleteCases[innerSumz, 0, 3]];
  printIL3D["At IL3D 8z, the inner sum list for ", indVar3, " after"<>
      " removing zero terms:"];
  printIL3D[listWOzerosz];

  (* Multiply the sum by the dependent variable and conduct derivatives *)
  (* with respect to x and y according to the outer sum.  The result is *)
  (* a list of terms.                                                   *)
  outerSumListx = Thread[Times[Map[
      Binomial[Part[#, 3] + Part[#, 2] + Part[#, 1], Part[#, 3]]*
      Binomial[Part[#, 2] + Part[#, 1], Part[#, 2]] &,
      positionNonZeroTermsx]*Map[Derivative[
      Part[#, 3], Part[#, 2], Part[#, 1]][u[dependentVar]]
          [indVar1, indVar2, indVar3] &,
      positionNonZeroTermsx], listWOzerosx]];
  printIL3D["At IL3D 9x, a list of outer sum terms for"<>
      " variable ", indVar1, " is:"];
  printIL3D[outerSumListx];

  outerSumListy = Thread[Times[Map[
      Binomial[Part[#, 3] + Part[#, 2] + Part[#, 1], Part[#, 2]]*
      Binomial[Part[#, 3] + Part[#, 1], Part[#, 1]] &,
      positionNonZeroTermsy]*Map[Derivative[
      Part[#, 3], Part[#, 2], Part[#, 1]][u[dependentVar]]
          [indVar1, indVar2, indVar3] &,
      positionNonZeroTermsy], listWOzerosy]];
  printIL3D["At IL3D 9y, a list of outer sum terms for"<>
      " variable ", indVar2, " is:"];
  printIL3D[outerSumListy];

  outerSumListz = Thread[Times[Map[
      Binomial[Part[#, 3] + Part[#, 2] + Part[#, 1], Part[#, 1]]*
      Binomial[Part[#, 3] + Part[#, 2], Part[#, 3]] &,
      positionNonZeroTermsz]*Map[Derivative[
      Part[#, 3], Part[#, 2], Part[#, 1]][u[dependentVar]]
          [indVar1, indVar2, indVar3] &,
      positionNonZeroTermsz], listWOzerosz]];
  printIL3D["At IL3D 9z, a list of outer sum terms for"<>
      " variable ", indVar3, " is:"];
  printIL3D[outerSumListz];

  (* Add the terms from the lists in the last step to get the result of *)
  (* applying the higher Euler operator for each variable.              *)
  integrandResultx = Expand[Apply[Plus, outerSumListx]];
  integrandResulty = Expand[Apply[Plus, outerSumListy]];
  integrandResultz = Expand[Apply[Plus, outerSumListz]];
  printIL3D["At IL3D 11, the integrands for the 3-D homotopy operator for"<>
      " dependent variable u[", dependentVar, "] are as follows:"];
  printIL3D["For variable ", indVar1, ": ", integrandResultx];
  printIL3D["For variable ", indVar2, ": ", integrandResulty];
  printIL3D["For variable ", indVar3, ": ", integrandResultz];

  If[globalHomotopyVerbose && TimeUsed[] - timeKeeper > 3.0,
    printOUT["All sums for dependent variable "<> ToString[dependentVar] <>
        " have been compiled."]
  ]; (* end If globalHomotopyVerbose *)

  (* Clear all variables not being passed back.                         *)
  Clear[derivativeList, orderList, maxix, maxiy, maxiz, partialWRTpartialList,
    innerSumListx, innerSumListy, innerSumListz, ix, iy, iz, kx, ky, kz,
    innerSumx, innerSumy, innerSumz, positionNonZeroTermsx, indVar2, indVar3,
    positionNonZeroTermsy, positionNonZeroTermsz, listWOzerosx,
    listWOzerosy, listWOzerosz, dummyderiv, indVar1, outerSumListx,
    outerSumListy, outerSumListz, printIL3D, printILLP];

  Return[{integrandResultx, integrandResulty, integrandResultz}]
]; (* end Module integrandLeftDerList3D *)

(**********************************************************************)
(*                                                                    *)
(*    BEGINNING THE REMOVAL OF CURL TERMS FROM THE HOMOTOPY RESULT    *)
(*                                                                    *)
(**********************************************************************)

(* ##########          Function: removeCurlTerms           ########## *)

(**********************************************************************)
(* removeCurlTerms[resultFromHomotopy, independVarlist]               *)
(* Purpose: To reduce the number of terms in the final homotopy       *)
(*          operator result to least amount possible.  This is a      *)
(*          modified form of extractSystemForCoefficient              *)
(* Input:   The result from the homotopy operator                     *)
(*          A list of independent variables                           *)
(* Output:  A simplified result                                       *)
(* Created: 20 September 2006 by DP at home                           *)
(* Code is in File:  newhomotopyoper/horecute.m                       *)
(* Major modifications to original file hosimhom.m                    *)
(*          made 18 May, 2007 by DP                                   *)
(* Last Modified:  8 June, 2009, 11:23 by DP at CSM                   *)
(**********************************************************************)

removeCurlTerms[resultFromHomotopy_, allIndependVarList_] :=
Module[{lengthIndVar, lengthComp, unknownConstants,  independVarList,
    homotopyExpnWithConstants, divergenceHomotopyWithK, divCheckDiv,
    divergenceHomotopyNoK, system, zeroEquations, zeroConstants,
    j1 = 1, j2, nonzeroConstants, solveConstantList, systemSoln,
    finalHomotopy, keepHomotopy, divCheck, ivar, i, i1 = 1, printRCT,
    homotopyResultPart, commonDenominator},

  If[debugRemoveCurlTerms, printRCT = Print, Clear[printRCT], Clear[printRCT]];
  printRCT["debug code: RCT, Function: removeCurlTerms File: horecute.m"];

  printRCT["At RCT IN, the vector from homotopy operator:"];
  printRCT[resultFromHomotopy];
  printRCT["At RCT IN, the list of independent variables is ",
      allIndependVarList];

  (* If no result is given, return immediately.                         *)
  If[resultFromHomotopy === {}, Return[resultFromHomotopy]];

  Clear[k];

  printOUT["The homotopy operator has finished inverting the divergence. "<>
      " The program is beginning the removal of curl terms."];

  (* Determine the number of independent variables.                     *)
  independVarList = Part[allIndependVarList, 1];
  lengthIndVar = Length[resultFromHomotopy];

  (* Assign independent variables to consistent terms.                  *)
  ivar[1] = Part[independVarList, 1];
  If[lengthIndVar > 1,  ivar[2] = Part[independVarList, 2]];
  If[lengthIndVar > 2,  ivar[3] = Part[independVarList, 3]];

  (* Find the number of terms in each component of the homotopy result. *)
  lengthComp[1] = If[Head[Part[resultFromHomotopy, 1]] === Plus,
      Length[Part[resultFromHomotopy, 1]], 1];
  lengthComp[2] = If[Head[Part[resultFromHomotopy, 2]] === Plus,
      Length[Part[resultFromHomotopy, 2]], 1];
  lengthComp[3] = If[lengthIndVar === 3,
      If[Head[Part[resultFromHomotopy, 3]] === Plus,
      Length[Part[resultFromHomotopy, 3]], 1], 0];

  (* If only one term exists in each component, there is no need to     *)
  (* go through the reduction process.                                  *)
  If[lengthComp[1] <= 1 && lengthComp[2] <= 1 && lengthComp[3] <= 1,
    Return[resultFromHomotopy]
  ]; (* end If lengthPart *)

  If[printRCT === Print, globalHomotopyVerbose = True];
  If[globalHomotopyVerbose,
    printOUT["The inverted divergence vector given by the homotopy operator"<>
        " has ", lengthIndVar, " components.  The first component has"<>
        " ", lengthComp[1], If[lengthComp[1] === 1, " term", " terms"],
        If[lengthIndVar === 3, ", "," and "],
        "the second component has ", lengthComp[2],
        If[lengthComp[2] === 1, " term", " terms"],
        If[lengthIndVar === 3, ", and the third component has "<>
        ToString[lengthComp[3]]<>
        If[lengthComp[3] === 1, " term.", " terms."], "."]]
  ]; (* end If globalHomotopyVerbose *)
  printOUT["If the number of terms in the components is large, say several"<>
      " hundred, the removal of curl terms may take a long time."];
  printOUT["------------------------------------------"];

  (* Distribute constants through the terms of the homotopy result.     *)
  {homotopyExpnWithConstants[1], homotopyExpnWithConstants[2],
      homotopyExpnWithConstants[3], unknownConstants} =
      applyConstantsToVector[resultFromHomotopy, lengthIndVar];
  printRCT["At RCT 1, the result from the homotopy operator with constants"<>
      " added to each term:"];
  printRCT[{homotopyExpnWithConstants[1], homotopyExpnWithConstants[2],
      homotopyExpnWithConstants[3]}];

  (* Take the divergence on the original homotopy result and the        *)
  (* homotopy result with constants.                                    *)
  If[lengthIndVar === 2,
    divergenceHomotopyWithK = divergence[{homotopyExpnWithConstants[1],
        homotopyExpnWithConstants[2], 0}, ivar[1], ivar[2], ivar[3]];
    divergenceHomotopyNoK = divergence[Append[resultFromHomotopy, 0],
        ivar[1], ivar[2], ivar[3]]
  ]; (* end If lengthIndVar === 2 *)
  If[lengthIndVar === 3,
    divergenceHomotopyWithK = divergence[{homotopyExpnWithConstants[1],
        homotopyExpnWithConstants[2], homotopyExpnWithConstants[3]},
        ivar[1], ivar[2], ivar[3]];
    divergenceHomotopyNoK = divergence[resultFromHomotopy,
        ivar[1], ivar[2], ivar[3]]
  ]; (* end If lengthIndVar === 3 *)
  divergenceHomotopyNoK = Expand[Factor[divergenceHomotopyNoK]];
  divCheckDiv = divergenceHomotopyNoK;
  printRCT["At RCT 2, after divergence is applied to the homotopy result"<>
      " with constants:"];
  printRCT[divergenceHomotopyWithK];
  printRCT["At RCT 3, after divergence is applied to the homotopy result"<>
      " without constants:"];
  printRCT[divergenceHomotopyNoK];

  (* Check for rational expressions.                                    *)
  denom = MemberQ[resultFromHomotopy, Power[n2__, n1_] /; n1 < 0 /;
      MemberQ[n2, u[_][__], {0, Infinity}] ||
      MemberQ[n2, Derivative[__][u[_]][__], {0, Infinity}], {0, Infinity}];
  printRCT["At RCT 4, the homotopy result is rational:  ", denom];

  If[denom,
    denominatorExpnNoK = Union[Map[Denominator,
        Apply[List, divergenceHomotopyNoK]]];
    denominatorExpnNoK = Apply[PolynomialLCM, denominatorExpnNoK];
    printRCT["At RCT 5, the denominator for the divergence without"<>
        " constants:"];
    printRCT[denominatorExpnNoK];

    divergenceHomotopyWithK = Map[Times[#, denominatorExpnNoK] &,
        divergenceHomotopyWithK];
    divergenceHomotopyWithK = Map[Factor, divergenceHomotopyWithK];
    printRCT["At RCT 6, the denominator from the divergence with no"<>
        " constants is multiplied to the divergence with constants:"];
    printRCT[divergenceHomotopyWithK];

    denominatorExpnWithK = Union[Map[Denominator,
        Apply[List, divergenceHomotopyWithK]]];
    denominatorExpnWithK = Apply[PolynomialLCM, denominatorExpnWithK];
    printRCT["At RCT 7, the denominator for the divergence with constants:"];
    printRCT[denominatorExpnWithK];

    divergenceHomotopyWithK = Map[Times[#, denominatorExpnWithK] &,
        divergenceHomotopyWithK];
    divergenceHomotopyWithK = Map[Factor, divergenceHomotopyWithK];
    printRCT["At RCT 8, the divergence with constants with the denominator"<>
        " removed:"];
    printRCT[divergenceHomotopyWithK];

    divergenceHomotopyNoK = Map[Times[#, denominatorExpnNoK*
        denominatorExpnWithK] &, divergenceHomotopyNoK];
    divergenceHomotopyNoK = Map[Factor, divergenceHomotopyNoK];
    printRCT["At RCT 9, the divergence without constants with the"<>
        " denominator removed:"];
    printRCT[divergenceHomotopyNoK]
  ]; (* end If denom *)

  (* Make a system of equations on the unknown constants.                *)
  (* The first case will try to construct a small system, assuming that  *)
  (* most of the k[i] are zero.  This may lead to an inconsistent system *)
  (* in some situations.  If so, then case two will be applied.          *)
  system =
      makeCoefficientEquations[divergenceHomotopyWithK, unknownConstants,
      divergenceHomotopyNoK, Flatten[allIndependVarList], "caseONE"];
  printRCT["At RCT 10, the system of unknown constants to be solved:"];
  printRCT[system];

  If[Not[MemberQ[system, False, {0, Infinity}]],
    solveConstantList = Cases[system, k[_], {0, Infinity}];
    printRCT["At RCT 11, a list of all constants to solve for. Order is"<>
        " not important in this case:"];
    printRCT[solveConstantList];

    zeroConstants = Complement[unknownConstants, solveConstantList];
    zeroConstants = Map[Rule[#, 0] &, zeroConstants];
    printRCT["At RCT 12, a list of constants taken to be zero:"];
    printRCT[zeroConstants];

    (* Find two sets of solutions by reversing the order of the solve   *)
    (* list.  These will be compared for balance.                       *)
    systemSoln[1] = solveCoefficientSystem[system, solveConstantList,
        zeroConstants];
    printRCT["At RCT 13, the first solutions from the coefficient system"<>
        " solver:"];
    printRCT[systemSoln[1]];

    systemSoln[2] = solveCoefficientSystem[system, Reverse[solveConstantList],
        zeroConstants];
    printRCT["At RCT 14, the alternate solutions from the coefficient system"<>
       " solver:"];
    printRCT[systemSoln[2]];

    If[systemSoln[1] === systemSoln[2], j2 = 1, j2 = 2];
    printRCT["At RCT 15, j2 = ", j2, ". If j2 = 2, the two solutions are"<>
        " different and a comparison will be made."];

    While[j1 <= j2,
      (* Replace constants in the homotopy operator with their solutions. *)
      finalHomotopy[j1] = Flatten[{homotopyExpnWithConstants[1],
          homotopyExpnWithConstants[2],
          homotopyExpnWithConstants[3]}] /. systemSoln[j1];
      printRCT["At RCT 16, the no. ",j1, " homotopy result after constants"<>
          " have been replaced using the system solution:"];
      printRCT[finalHomotopy[j1]];
      j1++
    ]; (* end While *)

    If[j2 === 2,
      keepHomotopy = compareCurlFreeSolutions[finalHomotopy[1],
          finalHomotopy[2]],
    (* else *)
      keepHomotopy = finalHomotopy[1]
    ], (* end If j2 *)
  (* else *)
    printRCT["At RCT 17, the system is inconsistent and cannot be solved."<>
      " A new system must be formed."];

    system = makeCoefficientEquations[divergenceHomotopyWithK,
        unknownConstants, divergenceHomotopyNoK, Flatten[allIndependVarList],
        "caseTWO"];
    printRCT["At RCT 18, the system of unknown constants to be solved:"];
    printRCT[system];

    (* Identify all equations in the system where the right hand side   *)
    (* is equal to zero.                                                *)
    zeroEquations = Cases[system,Equal[_,0], {0,2}];
    printRCT["At RCT 19, a list of equations where the right hand side is 0:"];
    printRCT[zeroEquations];

    (* Make a list of all constants, k[i], in zeroEquations.            *)
    zeroConstants = Union[Cases[zeroEquations, k[_], {0, Infinity}]];
    printRCT["At RCT 20, a list of all constants that are in the equations"<>
        " where the right hand side is 0:"];
    printRCT[zeroConstants];

    (* Make a list of all constants that don't show up in zeroEquations *)
    nonzeroConstants = Complement[unknownConstants, zeroConstants];
    printRCT["At RCT 21, a list of all constants that are NOT in the"<>
        " equations where the right hand side is 0:"];
    printRCT[nonzeroConstants];

    (* Solve the system, solving for the constants that don't appear in *)
    (* the zeroEquations first, then solving for all other constants.   *)
    (* Two different solve orders will be applied to give choices for   *)
    (* the solution.                                                    *)
    solveConstantList = Flatten[Append[nonzeroConstants, zeroConstants]];
    printRCT["At RCT 22, a list of all constants in the order they need"<>
      " to be solved for:"];
    printRCT[solveConstantList];

    systemSoln[1] = solveCoefficientSystem[system, solveConstantList, {}];
    printRCT["At RCT 23, the initial solution to the system:"];
    printRCT[systemSoln[1]];

    (* Reverse the order of the constants from the nonzero equations    *)
    (* and solve the system again.  This gives the alternative.         *)
    solveConstantList = Flatten[Append[Reverse[nonzeroConstants],
        zeroConstants]];
    printRCT["At RCT 24, a list of all constants in the order they need"<>
        " to be solved for with the leading constants reversed in order:"];
    printRCT[solveConstantList];

    systemSoln[2] = solveCoefficientSystem[system, solveConstantList, {}];
    printRCT["At RCT 25, the alternate solution to the system:"];
    printRCT[systemSoln[2]];

    If[systemSoln[1] === systemSoln[2], j2 = 1, j2 = 2];
    printRCT["At RCT 26, j2 = ", j2, ". If j2 = 2, the two solutions are"<>
        " different and a comparison will be made."];

    While[j1 <= j2,
      (* Replace constants in the homotopy operator with their          *)
      (* solutions.                                                     *)
      finalHomotopy[j1] = Flatten[{homotopyExpnWithConstants[1],
          homotopyExpnWithConstants[2],
          homotopyExpnWithConstants[3]}] /. systemSoln[j1];
      printRCT["At RCT 27, the no. ",j1, " homotopy result after constants"<>
            " have been replaced using the system solution:"];
      printRCT[finalHomotopy[j1]];
      j1++
    ]; (* end While *)

    (* Compare simplified results. Report the result with the most      *)
    (* balanced set of terms across the components.                     *)
    If[j2 === 2,
      keepHomotopy = compareCurlFreeSolutions[finalHomotopy[1],
          finalHomotopy[2]],
    (* else *)
      keepHomotopy = finalHomotopy[1]
    ] (* end If j2 *)
  ]; (* end If MemberQ *)

  printRCT["At RCT 28, the final homotopy result is:"];
  printRCT[keepHomotopy];

  (* Do a divergence check on the result after curl terms have been     *)
  (* removed.  If the test fails, return the result given by the        *)
  (* homotopy operator.                                                 *)
  divCheck = divergenceCheck[divCheckDiv, keepHomotopy, independVarList];
  printRCT["At RCT 29, Div[homotopyResult] - Div[curlTermsRemoved] = ",
      divCheck];

  If[divCheck =!= 0 || keepHomotopy === {},
    printOUT["The program was unable to remove curl terms from the vector"<>
        " provided by the homotopy operator.  The result with curl terms"<>
        " is begin returned."];
    keepHomotopy = resultFromHomotopy
   ];(* end divCheck *)

  Clear[lengthIndVar, lengthComp, homotopyResultList, unknownConstants,
      termsWithConstants, homotopyExpnWithConstants, divergenceHomotopyWithK,
      divergenceHomotopyNoK, system, zeroEquations, zeroConstants, j1, j2,
      nonzeroConstants, solveConstantList, systemSoln, finalHomotopy,
      divCheck, commonDenominator, ivar, i, divCheckDiv, printRCT];

  Return[keepHomotopy]
]; (* end Module removeCurlTerms *)

(* ##########     Function: makeCoefficientEquations       ########## *)

(**********************************************************************)
(* makeCoefficientEquations[expressionWithC, unknownCoefList,         *)
(*     expressionNoC, indVarList, case]                               *)
(*          modified form of extractSystemForCoefficients             *)
(* Purpose: To extract a linear system to solve for the coefficients  *)
(*          and match the coefficients to coefficients in an          *)
(*          equivalent system.                                        *)
(* Input:   The divergence of the homotopy expression after           *)
(*              coefficients, c[i], have been added                   *)
(*          The divergence of the homotopy expression with no         *)
(*              coefficients                                          *)
(*          List of unknown coefficients                              *)
(*          A list of independent variables                           *)
(*          A string telling whether the long or short method is      *)
(*              being used                                            *)
(* Output:  A system of coefficient equations                         *)
(* Created: 17 September 2006 by DP at CSM                            *)
(* Code is in File:  newhomotopyoper/homacoeq.m                       *)
(* Last Modified: 10 April, 2008, 9:06 by DP at CSM                   *)
(**********************************************************************)

makeCoefficientEquations[expressionWithK_, unknownCoefList_, expressionNoK_,
    indVarList_, case_String] :=
Module[{expandedExpressionWithK, expandedExpressionNoK, termsToRemove, dco,
    commonTermsListWithK, commonTermsListNoK, coefficientEqns, zeroConstants,
    n1, n2, n3, indVarToZeroRules, excludeList, printMCE},

  If[debugMakeCoefficientEquations, printMCE = Print, Clear[printMCE],
      Clear[printMCE]];
  printMCE["debug code: MCE, Function: makeCoefficientEquations,"<>
      " File: homacoeq.m"];

  expandedExpressionWithK = Expand[expressionWithK];
  expandedExpressionNoK = Expand[expressionNoK];

  printMCE["At MCE IN, the expression with k[i] coming into"<>
      " makeCoefficientsEquations: "];
  printMCE[expandedExpressionWithK];
  printMCE["At MCE IN, unknown coefficients to find: ", unknownCoefList];
  printMCE["At MCE IN, the expression without k[i]:"];
  printMCE[expandedExpressionNoK];
  printMCE["At MCE IN, the independent variable list contains ", indVarList];
  printMCE["At MCE IN, the case is ", case];

  (* Make a list of terms needed to identify common coefficients.        *)
  (* The first case will form a system for the constants where the rhs   *)
  (* is never zero.  This may or may not form a consistent system.       *)
  (* The second case sets up the whole coefficient system and will be    *)
  (* used when the first case fails.                                     *)
  Switch[case,
    "caseONE", termsToRemove =
        Reverse[Union[homotopyStripper[expandedExpressionNoK, indVarList]]],
    "caseTWO", termsToRemove =
        Reverse[Union[homotopyStripper[expandedExpressionWithK, indVarList]]]
  ]; (* end Switch *)
  printMCE["AT MCE 1, terms used as a basis for forming coefficient"<>
      " equations:"];
  printMCE[termsToRemove];

  excludeList = Flatten[{u[_][__], Derivative[__][u[_]][__], indVarList,
      Map[Power[#, _] &, indVarList]}];

  (* Match all terms in the expression with k[i] according to the terms  *)
  (* given in the previous list.                                         *)
  rememberCommonTerms =
  commonTermsListWithK = Apply[Plus, Map[Cases[expandedExpressionWithK,
      Times[___?(Not[MemberQ[#, Apply[Alternatives, excludeList],
      {0, Infinity}]] &), #], {0, 2}] &, termsToRemove], {1}];
  printMCE["AT MCE 2, a list of expressions taken from the divergence with"<>
      " k[i] corresponding to each term in the basis list:"];
  printMCE[commonTermsListWithK];

  (* Match all terms in the expression without c[i] according to the     *)
  (* terms given in the previous list.                                   *)
  expandedExpressionNoK = Map[Times[#, dco] &, expandedExpressionNoK];
  commonTermsListNoK = Apply[Plus, Map[Cases[expandedExpressionNoK,
     Times[___?(Not[MemberQ[#, Apply[Alternatives, excludeList],
      {0, Infinity}]] &), #], {0, 2}] &, termsToRemove], {1}];
  printMCE["AT MCE 3, a list of expressions taken from the divergence that"<>
      " has no k[i] corresponding to each term in the basis list:"];
  printMCE[commonTermsListNoK];

  commonTermsListWithK = Expand[Thread[Times[commonTermsListWithK,
      Map[Power[#, -1] &, termsToRemove]]]];
  printMCE["AT MCE 4, all dependent variables and their derivatives have"<>
      " been removed from the expressions with k[i]:"];
  printMCE[commonTermsListWithK];

  commonTermsListNoK = Expand[Thread[Times[commonTermsListNoK,
      Map[Power[# * dco, -1] &, termsToRemove]]]];
  printMCE["AT MCE 5, all dependent variables and their derivatives have"<>
      " been removed from the expressions with no k[i]:"];
  printMCE[commonTermsListNoK];

  (* Set the corresponding expressions from each list equal to each     *)
  (* other, forming coefficient equations.                              *)
  coefficientEqns = Thread[Equal[commonTermsListWithK, commonTermsListNoK]];
  printMCE["AT MCE 6, a system of equations formed from the coefficient"<>
      " expressions:"];
  printMCE[coefficientEqns];

  If[case === "caseONE",
    remainingTerms = Complement[Apply[List, expandedExpressionWithK],
        rememberCommonTerms /. Plus -> Sequence];
    printMCE["At MCE 7, all terms from the divergence with k[i] that do"<>
        " not form coefficient expressions:"];
    printMCE[remainingTerms];

    zeroConstants = Union[Cases[remainingTerms, k[_], {0, Infinity}]];
    zeroConstants = Map[Rule[#, 0] &, zeroConstants];
    printMCE["AT MCE 8, divergence constants which are set automatically"<>
        " to zero (Note: This may form an inconsistent system.):"];
    printMCE[zeroConstants];

    coefficientEqns = coefficientEqns /. zeroConstants;
    printMCE["AT MCE 9, a shortened system of equations with both sides"<>
        " reduced to coefficients only:"];
    printMCE[coefficientEqns];

    coefficientEqns = coefficientEqns /. Equal[0, _] -> False;
    printMCE["AT MCE 10, 0 = parameter must be changed to false:"];
    printMCE[coefficientEqns]
  ]; (* end If case *)

  Clear[expandedExpressionWithK, expandedExpressionNoK, termsToRemove, r,
      commonTermsListWithK, commonTermsListNoK, zeroConstants, n1, n2, n3,
      tim, printMCE];

  Return[coefficientEqns]
]; (* end Module makeCoefficientEquations *)

(* ##########       Function: solveCoefficientSystem       ########## *)

(**********************************************************************)
(* solveCoefficientSystem[coefSystem, coefOrder, knownZeros]          *)
(* Purpose: To solve a linear system of constants.  In this case, all *)
(*          free constants in the solution will be set to zero        *)
(*          automatically. Zero solutions will be listed.             *)
(* Input:   The system to be solved                                   *)
(*          A list of constants to be solved for given in the order   *)
(*              that they should be solved for                        *)
(*          A list of solutions known to be zero                      *)
(* Output:  The solution to the system with every constant accounted  *)
(*          for.                                                      *)
(* Created: 21 May 2007 by DP at CSM                                  *)
(* Code is in File:  newhomotopyoper/hosocosy.m                       *)
(* Last Modified:  10 March, 2008, 16:05 by DP at CSM                 *)
(**********************************************************************)

solveCoefficientSystem[coefSystem_, coefOrder_List, knownZeros_List] :=
Module[{systemSoln, undeterminedConstants, n1, printSCS},

  If[debugSolveCoefficientSystem, printSCS = Print, Clear[printSCS],
      Clear[printSCS]];
  printSCS["debug code: SCS, Function: solveCoefficientSystem,
      File: hosocosy.m"];

  printSCS["At SCS IN, the given system of constants is:"];
  printSCS[coefSystem];
  printSCS["At SCS IN, the constants to be solved for in the order that"<>
      " they are to be solved:"];
  printSCS[coefOrder];
  printSCS["At SCS IN, the solutions known to be zero:"];
  printSCS[knownZeros];

  (* Solve the system given to the function.                            *)
  systemSoln = Flatten[Solve[coefSystem, coefOrder]];
  printSCS["At SCS 1, the initial solution to the system using the"<>
        " reversed constants:"];
  printSCS[systemSoln];

  (* Collect all undetermined constants.                                *)
  undeterminedConstants = Union[Cases[systemSoln /.
      Rule[_, n1_] -> n1, k[_], {0, Infinity}]];
  printSCS["At SCS 2, all undetermined constants in the system solution:"];
  printSCS[undeterminedConstants];

  (* Formulate rules setting all undetermined constants to zero.        *)
  undeterminedConstants = Map[Rule[#, 0] &, undeterminedConstants];
  printSCS["At SCS 3, all undetermined constants set to zero:"];
  printSCS[undeterminedConstants];

  (* Apply the undetermined constants rules to the solution.            *)
  systemSoln = systemSoln /. undeterminedConstants;
  printSCS["At SCS 4, the solution to the system after undetermined"<>
      " constants have been replaced by zeros:"];
  printSCS[systemSoln];

  (* Make a list showing the solutions to all k[i].                   *)
  systemSoln = Union[systemSoln, undeterminedConstants, knownZeros];
  printSCS["At SCS 5, all solutions listed together:"];
  printSCS[systemSoln];

  Clear[undeterminedConstants, n1, printSCS];

  Return[systemSoln]
]; (* end Module solveCoefficientSystem *)

(* ##########      Function: compareCurlFreeSolutions      ########## *)

(**********************************************************************)
(* compareCurlFreeSolutions[invertedDiv1, invertedDiv2]               *)
(* Purpose: To compare two possibilities for the final form for an    *)
(*          inverted divergence and recommend the one that has the    *)
(*          most balanced solution.                                   *)
(* Input:   A pair of inverted divergences                            *)
(* Output:  The inverted divergence with the most balanced form.      *)
(* Created: 21 May 2007 by DP at CSM                                  *)
(* Code is in File:  newhomotopyoper/hoccfsol.m                       *)
(* Last Modified:  10 March, 2008, 17:13 by DP at CSM                 *)
(**********************************************************************)

compareCurlFreeSolutions[invertedDiv1_, invertedDiv2_] :=
Module[{mylength, lengthCheck, findDiffInLength, lengthDiffs, i1, i2,
    keepHomotopy, printCCFS},

  If[debugCompareCurlFreeSolutions, printCCFS = Print, Clear[printCCFS],
      Clear[printCCFS]];
  printCCFS["debug code: CCFS, Function: compareCurlFreeSolutions,
      File: hoccfsol.m"];

  printCCFS["At CFFS IN, the first inverted divergence is:"];
  printCCFS[invertedDiv1];
  printCCFS["At CFFS IN, the first inverted divergence is:"];
  printCCFS[invertedDiv2];

  (* Compare simplified results. Report the result with the most        *)
  (* balanced set of terms across the components.                       *)
  mylength := If[Head[#] === Plus, Length[#], 1] &;
  lengthCheck[1] = Map[mylength, invertedDiv1];
  lengthCheck[2] = Map[mylength, invertedDiv2];
  printCCFS["At CCFS 1, the first simplification has lengths ", lengthCheck[1],
      " and the second simplification has lengths ", lengthCheck[2]];

  findDiffInLength := Abs[#1 - #2] &;
  lengthDiffs[1] = Table[Part[Append[
      lengthCheck[1], Part[lengthCheck[1], 1]], i1],
      {i2, 1, Length[lengthCheck[1]] - 1}, {i1, i2, i2 + 1}];
  lengthDiffs[2] = Table[Part[Append[
      lengthCheck[2], Part[lengthCheck[2], 1]], i1],
      {i2, 1, Length[lengthCheck[2]] - 1}, {i1, i2, i2 + 1}];
  printCCFS["At CCFS 2, lengths of components paired to find differences"<>
      " lengths for the first simplification are ", lengthDiffs[1],
      " and for the second simplification are ", lengthDiffs[2]];

  lengthDiffs[1] = Max[Apply[findDiffInLength, lengthDiffs[1], {1}]];
  lengthDiffs[2] = Max[Apply[findDiffInLength, lengthDiffs[2], {1}]];
  printCCFS["At CCFS 3, the maximum difference in length for the first"<>
      " simplification is ", lengthDiffs[1], " and for the second"<>
      " simplification is ", lengthDiffs[2]];

  keepHomotopy = If[lengthDiffs[1] <= lengthDiffs[2],
      invertedDiv1, invertedDiv2];
  printCCFS["At CCFS 4, the choice for the final homotopy result:"];
  printCCFS[keepHomotopy];

  Clear[mylength, lengthCheck, findDiffInLength, lengthDiffs, i1, i2,
      printCCFS];

  Return[keepHomotopy]
]; (* end Module compareCurlFreeSolutions *)

(* ##########              Function: divergence              ########## *)

divergence[expr_List, iv1_, iv2_:0, iv3_:0] :=
    D[Part[expr, 1], iv1] + If[iv2 === 0, 0, D[Part[expr, 2], iv2]] +
    If[iv3 === 0, 0, D[Part[expr, 3], iv3]];

(* ##########           Function: divergenceCheck            ########## *)

(************************************************************************)
(* divergenceCheck[expression, homotopyOpResult, indepVarList]          *)
(* Purpose:  To take the derivative of divergence on the result         *)
(*           provided by the homotopy operator and check to see         *)
(*           if it is the same as the given expression.                 *)
(* Input:    The original given expression                              *)
(*           The calculated result from the homotopy operator           *)
(*           A list of independent variables                            *)
(* Output:   Yes or No.                                                 *)
(* Created:  30 August, 2007 by DP                                      *)
(* Code is in File:  newhomotopyoper/hodivchk.m                         *)
(* Last Modified: 20 March, 2008, 11:30 by DP at CSM                    *)
(************************************************************************)

divergenceCheck[expression_, homotopyOpResult_, indepVarList_] :=
Module[{numIndVars = Length[indepVarList], newDiv, divCheck},

  If[debugDivergenceCheck, printDC = Print, Clear[printDC], Clear[printDC]];
  printDC["debug code: DC, Function: divergenceCheck,"<>
      " File: hodivdhk.m"];

  printDC["At DC IN, the given expression is:"];
  printDC[expression];
  printDC["At DC IN, the homotopy result is:"];
  printDC[homotopyOpResult];
  printDC["At DC IN, the list of independent variables given is ",
      indepVarList];

    newDiv = divergence[homotopyOpResult, Apply[Sequence, indepVarList]];
    printDC["At DC 1, the "<>If[numIndVars === 1, "derivative", "divergence"]<>
        " of the result produced by the homotopy operator:"];
    printDC[newDiv];

    printDC["Checking the "<>If[numIndVars === 1, "derivative", "divergence"]<>
        " against the given expression:"];
    divCheck = Expand[expression - newDiv];
    printDC["At DC 2, applying Expand to the difference:"];
    printDC[divCheck];

    If[divCheck =!= 0,
      divCheck = Factor[divCheck];
      printDC["At DC 3, applying Factor to the difference:"];
      printDC[divCheck]
    ]; (* end If divCheck *)

    If[divCheck =!= 0,
      divCheck = Expand[TrigToExp[divCheck]];
      printDC["At DC 4, applying TrigToExp to the difference:"];
      printDC[divCheck]
    ]; (* end If divCheck *)

    If[divCheck =!= 0,
      divCheck = FullSimplify[divCheck, TimeConstraint -> 100];
      printDC["At DC 5, applying FullSimplify to the difference:"];
      printDC[divCheck]
    ]; (* end If divCheck *)

  Clear[numIndVars, newDiv];

  If[divCheck === 0,
    Return[0],
    Return[1]]
]; (* end Module divergenceCheck *)

printLF["divergenceCheck[expression, homotopyOpResult, indepVarList]"<>
    " successfully loaded."];

(* ##########              Function: pdeform               ########## *)

(**********************************************************************)
(* Purpose: Changes an expression: removes functional dependence, and *)
(*          prints derivatives as subscripts.                         *)
(* Authors: Tracy Otto, Tony Miller, and U. Goktas                    *)
(* Input:   An expression                                             *)
(* Output:  Same expression with all of the changes made              *)
(* Created: 1995 at CSM                                               *)
(* Adapted From:  data/newzeala/reu2004/reu2004-0811/software/        *)
(*                whlepedf.m, Created 1995 at CSM                     *)
(* Code is in File:  conservationlaws/dplepdef.m                      *)
(* Last Modified:  19 December, 2006, 11:40 by DP at CSM              *)
(**********************************************************************)

pdeform[expres_] := expres /. {
        Derivative[n__][u[k_]][x__] :>
                SequenceForm[u, Subscript[k],Subscript[","],
                        Subscript[
                         SequenceForm @@ Flatten[Table[#[[1]], {#[[2]]}]& /@
                                        Transpose[{{x}, {n}}]]]],
        u[n_][x__] :> SequenceForm[u,Subscript[n]] }; (* end pdeform *)

(* ##########       Function: printSystemInPdeForm         ########## *)

(**********************************************************************)
(* printSystemInPdeForm[]:                                            *)
(* Author: Mike Colagrosso                                            *)
(* Adapted From:  data/newzeala/reu2004/reu2004-0811/software/        *)
(*                3DCO0510.m                                          *)
(* Major Changes Made: 15 May, 2007, 12:03 by DP at CSM               *)
(**********************************************************************)

(* Major changes to function made since differential expressions and     *)
(* independent variables are passed by the function now.                 *)
printSystemInPdeForm[diffFunctionList_, name_, note_, numDepVars_] := Module[
    {inPDEform, eqnIter},
  Print[If[Length[diffFunctionList] > 1, "System given for ", ""],
      name, "."];
  If[Head[note] === String, Print["NOTE: ", note]];
  inPDEform = Map[pdeform, diffFunctionList];
  For[eqnIter = 1, eqnIter <= numDepVars, eqnIter++,
    Print["\nEquation ",eqnIter," of the system with ",numDepVars," equation",
        If[Length[diffFunctionList] > 1, "s", ""], ":"];
    Print["\n  ", Part[inPDEform, eqnIter]," = 0\n\n"];
  ]; (* end For *)
]; (* end Module printSystemInPdeForm *)

printLF["printSystemInPdeForm: loaded successfully"];

(* ##########           Function: homotopyDriver             ########## *)

(************************************************************************)
(* homotopyDriver[expression, indepVariableList, (optional) parameters] *)
(* Purpose:  To integrate a 1D expression or to invert the divergence   *)
(*           operator on a 2D or 3D expression by applying a homotopy   *)
(*           operator.                                                  *)
(* Input:    An expression                                              *)
(*           A list of independent variables                            *)
(*           A list of parameters (can be given as an option)           *)
(* Output:   The result of applying the higher Euler operators to the   *)
(*           expression.                                                *)
(* Created:  23 March, 2006 by DP                                       *)
(* Code is in File:  newhomotopyoper/hodriver.m                         *)
(* Last Modified: 2 October, 2008, 11:20 by DP at CSM                   *)
(************************************************************************)

homotopyDriver[expression_, indepVariableList_, parameters_:{}] :=
Module[{expressionAltered, integrability, interimResult, indepVariableList2,
    integrableExpression, finalResult, lengthIndVarList, expressionDenom,
    finalCheck, ddd, j, newIndVarList, printHO},

  If[debugHomotopyDriver, printHO = Print, Clear[printHO], Clear[printHO]];
  printHO["debug code: HO, Function: homotopyDriver, File: hodriver.m"];

  (* TO DO needed for test program - to be removed later                 *)
  paramrules = {};

  printHO["At HO IN, the given expression is:"];
  printHO[expression];
  printHO["At HO IN, the independent variable list is ", indepVariableList];
  printHO["AT HO IN, the parameters have been given as ", parameters];

  If[globalHomotopyVerbose,
    printOUT["The expression given to the program is"];
    printOUT[expression];
    printOUT["------------------------------------------"]
  ]; (* end globalHomotopyVerbose *)

  (* Check to see that the independent variables are in a list.         *)
  indepVariableList2 = If[Head[indepVariableList] === Symbol,
      {indepVariableList}, indepVariableList, indepVariableList];
  printHO["At HO 1, the independent variable list in list form is ",
       indepVariableList2];
  lengthIndVarList = Length[indepVariableList2];

  (* If the expression given to the homotopy operator is 0, Return 0.   *)
  If[expression === 0,
    Return[Switch[lengthIndVarList, 1, 0, 2, {0, 0}, 3, {0, 0, 0}]]
    ];

  (* Reformat the given expression so that all dependent variables have *)
  (* the form u[_] and and any auxillary independent variables are      *)
  (* renamed so that they appear to be removed from the expression.     *)
  expressionAltered = expressionVarReassign[expression,
      indepVariableList2];
  printHO["At HO 4, the expression with variables changed is:"];
  printHO[Part[expressionAltered, 1]];
  printHO["At HO 5, the number of dependent variables is:"];
  printHO[Part[expressionAltered, 3]];
  printHO["At HO 6, the reverse change rule for dependent variables is "<>
      "", Part[expressionAltered, 2]];
  printHO["At HO 7, the full list of independent variables is "<>
      "", Part[expressionAltered, 4]];
  printHO["At HO 8, the list of integrating variables will be carried"<>
      " forward as:"];
  printHO[Part[expressionAltered, 5]];
  printHO["At HO 9, the reverse rule(s) for reinserting terms with"<>
      " derivatives on auxillary variables:"];
  printHO[Part[expressionAltered, 6]];

  (* Determine if the expression is integrable (1D) or if it is a       *)
  (* divergence (2D and 3D).                                            *)
  integrability = integrabilityTest[Part[expressionAltered, 1],
      Part[expressionAltered, 3], Part[expressionAltered, 5],
      parameters];
  printHO["At HO 10, the function is ", Part[integrability, 1]];

  (* If the expression is not integrable or not a divergence, stop      *)
  (* the program.                                                       *)
  If[Part[integrability, 1] =!= "Integrable", Abort[]];

  (* If the expression is integrable or is a divergence, pass it to the *)
  (* homotopy operator to integrate it or invert the divergence.        *)
  integrableExpression = If[Part[integrability, 2] =!= {},
      Part[expressionAltered, 1] /. Part[integrability, 2],
      Part[expressionAltered, 1], Part[expressionAltered, 1]];
  interimResult = rationalExpressionCheck[integrableExpression,
      Part[expressionAltered, 5], Part[expressionAltered, 3]];
  printHO["At HO 11, the interim result (before curl terms are removed):"];
  printHO[interimResult];

  (* If homotopy operator did not work, return empty list.              *)
  If[interimResult === {},
    Return[{}]
  ]; (* end If interimResult *)

  (* If the interim result is an inverted divergence, remove the curl   *)
  (* terms.                                                             *)
  If[lengthIndVarList > 1,
    newIndVarList = {Part[expressionAltered, 5],
        Complement[Part[expressionAltered, 4], Part[expressionAltered, 5]]};
    finalResult = removeCurlTerms[ExpandAll[interimResult], newIndVarList],
  (* else *)
    finalResult = interimResult;
  ]; (* end If lengthIndVarList *)
  printHO["At HO 12, the result returned after all curl terms have been"<>
      " removed:"];
  printHO[finalResult];

  (* Revert the expression back to the variables given in the problem.  *)
  (* Part[expressionAltered, 2] is the inversion rule for the           *)
  (* dependent variables. Part[expressionAltered, 6] returns terms with *)
  (* derivatives on auxillary variables.                                *)
  If[Part[expressionAltered, 6] =!= {},
    finalResult = returnAuxIndependentVars[finalResult,
        Part[expressionAltered, 4], Part[expressionAltered, 5],
        Part[expressionAltered, 6]];
    printHO["At HO 13, the final result after the auxillary independent"<>
        " variables have been returned:"];
    printHO[finalResult]
  ]; (* end If Part[expressionAltered, 6] *)

  printHO["At HO 14, the rule to change the dependent variables to their"<>
      " original form is given as ", Part[expressionAltered, 2]];

  (* Return the final result in factored form.                          *)
  (* BEGIN TAKE OUT
  finalResult = Map[Factor, finalResult];
  printHO["At HO 15, after factoring, the final result is:"];
  printHO[finalResult];
  END TAKE OUT *)

  finalResult = MapAll[Apart, finalResult];
  printHO["At HO 15, after factoring, the final result is:"];
  printHO[finalResult];

  finalResult = finalResult /. Part[expressionAltered, 2];
  printHO["At HO 16, the final result has now been changed to:"];
  printHO[finalResult];

  If[Part[integrability, 2] =!= {},
    Print["In the result that follows, parameters have been reassigned"<>
        " under the rules given by:"];
    Print[Part[integrability, 2]];
    paramrules = Part[integrability, 2]
  ]; (* end If Part[integrability, 2] *)
  printHO["At HO 17, replacing parameters with their restricting values:"];
  printHO[finalResult];

  (* Check the final result by comparing its derivative/divergence to   *)
  (* the expression given to the program.  They should match.           *)
  finalCheck = divergenceCheck[expression /. paramrules,
      finalResult, Part[expressionAltered, 5]];

  printHO["At HO 18, the final check on the result, comparing the derivative"<>
      " or divergence of the result to the given expression produces:"];
  printHO[finalCheck];

  If[finalCheck === 1,
    printOUT["The program has been unable to verify that the"<>
        If[lengthIndVarList === 1, " derivative", " divergence"]<>
        " of the result determined by the homotopy operator matches"<>
        " that of the given expression.  The result may not be accurate"<>
        " or the program may not be able to correctly simplify the"<>
        " expression calculated in the final check."],
  (* else *)
    printOUT["Check: Given Expression - "<> If[lengthIndVarList === 1,
        "D[", "Div["] <> "Homotopy Result] = ", finalCheck ]
  ]; (* end If finalCheck *)

  finalResult = If[lengthIndVarList === 1, Part[finalResult, 1],
      finalResult, finalResult];

  printio["Leaving homotopyDriver with final result:"];
  printio[finalResult];

  Clear[expressionAltered, integrability, interimResult, indepVariableList2,
      integrableExpression, finalCheck, lengthIndVarList, expressionDenom,
      ddd, j];

  Return[finalResult]
]; (* end Module homotopyDriver *)

(**********************************************************************)
(* The list of Options controls print  and debug statements.          *)
(* The option is given as a rule in the opening function of the       *)
(* program.                                                           *)
(**********************************************************************)

Options[HomotopyIntegrator] = {Verbose -> False, printHomotopyDriver -> False,
  printStandAloneOutput -> True, printExpressionVarReassign -> False,
  printChangeIndpendentParameters -> False, printIntegrabilityTest -> False,
  printReturnAuxIndependentVars -> False, printIntegrandLoop -> False,
  printHomotopyFinalIntegration -> False, printParameterTest -> False,
  printIntegrandLeftDerList -> False, printApplyConstantsToVector -> False,
  printDegreeOperator -> False, printVariationalDerivative -> False,
  printDivergenceCheck -> False, printRemoveCurlTerms -> False,
  printTakeOnlyNumerator -> False, printMakeCoefficientEquations -> False,
  printSolveCoefficientSystem -> False, printHomotopyStripper -> False,
  printCompareCurlFreeSolutions -> False, printShiftDenominatorTerm -> False,
  printRationalExpressionCheck -> False};

(* ##########        Function: HomotopyIntegrator          ########## *)

(**********************************************************************)
(* HomotopyIntegrator[pde, indepVariables, parameters, options]       *)
(*                                                                    *)
(* This function sets all options, then starts the program by calling *)
(* the driver function.                                               *)
(**********************************************************************)

HomotopyIntegrator[pde_, indepVariables_,
    parameters___?(MemberQ[{List, Symbol}, Head[#]] &), opts___?OptionQ] :=
Module[{result},

  printio["Entering HomotopyIntegrator with the given expression as:"];
  printio[pde];
  printio["The independent variable list is ", indepVariables];
  printio["Parameters have been given as ", parameters];
  printio["The options given are:"];
  printio[opts];
  Off[Solve::svars];

  (* Set up switches for turning on print and debug statements.         *)
  (* Verbose -> True prints a running dialog showing the progression    *)
  (* of  a function through the homotopy operator                       *)
  globalHomotopyVerbose = Verbose /. {opts} /. Options[HomotopyIntegrator];

  (* printStandAloneOutput -> False suppresses all output.              *)
  printOUT = printStandAloneOutput /. {opts} /. Options[HomotopyIntegrator];
  If[printOUT, printOUT = Print, Clear[printOUT], Clear[printOUT]];

  debugHomotopyDriver = printHomotopyDriver /. {opts} /.
      Options[HomotopyIntegrator];
  debugExpressionVarReassign = printExpressionVarReassign /. {opts} /.
      Options[HomotopyIntegrator];
  debugChangeIndependentParameters = printChangeIndpendentParameters /.
      {opts} /. Options[HomotopyIntegrator];
  debugReturnAuxIndependentVars = printReturnAuxIndependentVars /.
      {opts} /. Options[HomotopyIntegrator];
  debugHomotopyFinalIntegration = printHomotopyFinalIntegration /.
      {opts} /. Options[HomotopyIntegrator];
  debugShiftDenominatorTerm = printShiftDenominatorTerm /.
      {opts} /. Options[HomotopyIntegrator];
  debugApplyConstantsToVector = printApplyConstantsToVector /.
      {opts} /. Options[HomotopyIntegrator];
  debugIntegrabilityTest = printIntegrabilityTest /. {opts} /.
      Options[HomotopyIntegrator];
  debugIntegrandLeftDerList = printIntegrandLeftDerList /. {opts} /.
      Options[HomotopyIntegrator];
  debugIntegrandLoop = printIntegrandLoop /. {opts} /.
      Options[HomotopyIntegrator];
  debugParameterTest = printParameterTest /. {opts} /.
      Options[HomotopyIntegrator];
  debugDegreeOperator = printDegreeOperator /. {opts} /.
      Options[HomotopyIntegrator];
  debugVariationalDerivative = printVariationalDerivative /. {opts} /.
      Options[HomotopyIntegrator];
  debugDivergenceCheck = printDivergenceCheck /. {opts} /.
      Options[HomotopyIntegrator];
  debugRemoveCurlTerms = printRemoveCurlTerms /. {opts} /.
      Options[HomotopyIntegrator];
  debugTakeOnlyNumerator = printTakeOnlyNumerator /. {opts} /.
      Options[HomotopyIntegrator];
  debugMakeCoefficientEquations = printMakeCoefficientEquations /.
      {opts} /. Options[HomotopyIntegrator];
  debugSolveCoefficientSystem = printSolveCoefficientSystem /. {opts} /.
      Options[HomotopyIntegrator];
  debugCompareCurlFreeSolutions = printCompareCurlFreeSolutions /.
      {opts} /. Options[HomotopyIntegrator];
  debugHomotopyStripper = printHomotopyStripper /. {opts} /.
      Options[HomotopyIntegrator];
  debugRationalExpressionCheck = printRationalExpressionCheck /. {opts} /.
      Options[HomotopyIntegrator];

  result = homotopyDriver[pde, indepVariables, parameters];

  If[printOUT === Print, Off[Solve::svars], On[Solve::svars]];

  Return[result]
]; (* end Module HomotopyIntegrator *)

Print["Code HomotopyIntegrator.m of July 20, 2009 successfully loaded."];

(* ************************** end of all ******************************* *)

