
(* Time-stamp: 14 October, 2010, 15:55 Mountain Daylight Time *)

(*****************************************************************************)
(*                                                                           *)
(*            *** M A T H E M A T I C A   P R O G R A M ***                  *)
(*                                                                           *)
(*      SYMBOLIC COMPUTATION SOLITARY WAVE AND SOLITON SOLUTIONS             *)
(*         OF CERTAIN CLASSES OF NONLINEAR EVOLUTION EQUATIONS               *)
(*                 USING A HOMOGENIZATION METHOD                             *)
(*           A.K.A. A SIMPLIFIED VERSION OF HIROTA'S METHOD                  *)
(*                                                                           *)
(* Program name: Homogenize-And-Solve.m                                      *)
(*                                                                           *)
(* Purpose: Compute solitary wave solutions and 1, 2, or 3 soliton solutions *)
(*          of certain classes of nonlinear evolution PDEs in (1+1)          *)
(*          dimensions by transforming the given PDE into a PDE that is      *)
(*          homogeneous of degree in a new dependent variable.               *)
(*          Only logarithmic derivative transformations for the dependent    *)
(*          are considered.                                                  *)
(*                                                                           *)
(* Input to Homogenize-And-Solve.m:                                          *)
(*          Single nonlinear evolution equation of any order, in any degree, *)
(*          with polynomial terms.                                           *)
(*          Independent variables x and t, dependent variable u(x,t).        *)
(*          Any parameters must be constant.                                 *)
(*          A typical PDE must be of the form                                *)
(*          u_t = f(u, u_x, ... , u_{nx}), where n is the order of the PDE.  *)
(*                                                                           *)
(* Output : Solitary wave solution or soliton solutions of order 1, 2, or 3  *)
(*          if they exist.                                                   *)
(*                                                                           *)
(* Language : Mathematica, version 8                                         *)
(*                                                                           *)
(* Authors: Andrew Cook, Willy Hereman                                       *)
(*                                                                           *)
(*          Department of Applied Mathematics and Statistics                 *)
(*          Colorado School of Mines                                         *)
(*          Golden, CO 80401-1887, U.S.A.                                    *)
(*                                                                           *)
(*          and Unal Goktas                                                  *)
(*                                                                           *)
(*          Department of Computer Engineering                               *)
(*          Turgut Ozal University                                           *)
(*          Kecioren, Ankara 06010, Turkey                                   *)
(*                                                                           *)
(* Version 1: Released on October 14, 2012                                   *)
(*                                                                           *)
(*****************************************************************************)

(* Old notes Andrew Cook *)
(* 
   PROGRAM VERSION: acHo0904.m 

   AUTHORS: Andrew Cook, Willy Hereman, and Unal Goktas

   LAST UPDATED: October 14, 2012 at 4:30p.m. by Hereman in Boulder

   Previously updated: September 04, 2012 at 17:00 by Cook in Golden

   PURPOSE: Given a nonlinear partial differential equation (PDE), 
            the program attempts to apply a logarithmic derivative 
            transformation to the PDE to transform the given PDE into a PDE 
            that is homogenous in degree in the new dependent variable . 
            Subsequently, it attempts to find solitary wave and soliton 
            solutions to the homogenized PDE.

   INPUT:   A PDE in (1+1) dimensions with dependent variable u(x,t)

   OUTPUT:  A menu with options of PDEs, the corresponding homogenized PDE,
            and a soltary wave or soliton solution to the given PDE.

   LANGUAGE: Mathematica version 8

   TO DO: * Adjust 'pdeform' module to work with Mathematica 8
          * Check code to verify that it works for all PDEs
          * Get u[x,t] soln in terms of hyperbolic tangent
          * Make more data files
          * Edit code so it can handle systems of equations
*)

(* Notes about program modification written by Andrew Cook *)

(* Edit 5/28 to include proper naming of variables & include time used. 
      Removed unnecessary 'if' statement that displayed homogenized PDE. 
      Minor changes in print statement wording. Added print statement 
      to show file loaded successfully. Improvements to header display made. 
      Added debug flag. Simplified 'k' print statement *)
(* Edit 5/29 to include menu interface. Added functions printpage and menu, 
      originally created in Hereman's CONDENS.M file. Edited these 
      functions to fit with program better *)
(* Edit 5/30 to remove tabs and improve line length *)
(* Edit 6/2 to copy 'highestpair' and 'length' modules from condens.m *)
(* Edit 6/3 to add more detail to module headers. Modified 
      highestorder module to simplify module and update debug flags *)
(* Edit 6/4 to improve bracket matching *)
(* Edit 6/5 on 'homogenize' module to separate x and t derivatives prior to 
      integration. Added 'welcome' module. *)
(* Edit 6/7 to add data files, clean up homogenize module *)
(* Edit 6/10 to add additional print statements to homogenize module *)
(* Edit 6/12 to add looping logic to integration in homogenize module *)
(* Edit 6/15 to let user pick homogenized PDE *)
(* Edit 6/17 to add framework for solitonsoln module *)
(* Edit 6/18 to add more print statements to solitonsoln module *)
(* Edit 6/20 to improve epsilon1 logic for solitonsoln module *)
(* Edit 6/21 to add more appropriate f[x,t] and additional 
      logic behind that to solitonsoln module *)
(* Edit 6/23 for picking g[2][x,t] for epsilon2 in solitonsoln *)
(* Edit 6/24 to clarify print statements in solitonsoln *)
(* Edit 6/26 to add improved logic to epsilon2 scenario *)
(* Edit 7/1 for epsilon3 logic in solitonsoln *)
(* Edit 7/2 to begin cleaning up solitonsoln *)
(* Edit 7/7 to add epsilon 4 & 5 logic & increase speed for K-K *)
(* Edit 7/9 to work on while loop for solitonsoln *)
(* Edit 7/20 to being to add stopping criterion in solitonsoln *)
(* Edit 7/21 to change return of homogenize, finish stopping criterion, adjust
      looping structure in solitonsoln, plug f into u to give soln *)
(* Edit 7/22 to create solitonsolnOper module, 
      add local vars to solitonsoln *)
(* Edit 7/23 to add basic epsilon2 logic and adjust replacement rules *)
(* Edit 7/24 to attempt to better phiReplacementRules in solitonsolnOper *)
(* Edit 7/25 to clean up code, test u and f solutions *)
(* Edit 7/26 to add operator rules and adjust nonlinear operator *)
(* Edit 7/29 to adjust operator/derivatives logic in solitonsolnOper *)
(* Edit 7/30 to improve variable changing rules for nonlinear operator, 
      add epsilon3 logic, add stopping rules in solitonsolnOper module *)
(* Edit 7/31 to add looping logic for epsilon2+, add f[x,t] and
      other checks in solitonsolnOper *)
(* Edit 8/1 to replace logic for checking u[1][x,t] soln *)
(* Edit 8/4 to display E^Theta instead of Phis in F[x,t] soln *)
(* Edit 8/6 to make code more user friendly *)
(* Edit 8/8 to change logic for dealing with RHS in solitonsolnOper *)
(* Edit 8/9 to continue edits from 8/8 and attempt to improve speed *)
(* Edit 8/10 to delete the zero a[i][j] terms appearing in solitonsolnOper.
      Try to improve speed of monomialList by adding 'stripper' module, but
      will likely delete due to inefficiency. Rules have been added instead *)
(* Edit 8/12 to verify code works for most cases, added comments in code 
      up until epsilon2 while loop in solitonsolnOper *)
(* Edit 8/13 to continue adding comments in code *)
(* Edit 8/16 to continue with comments, also added h_cal.m and h_cal2.m.
      Added extra condition in homogenize when solving so K != 0.
      Changed logic to determine which formula to use in solitonsolnOper *)
(* Edit 8/17 to check to see if given PDE is already homogenized. Changed
      yesterday's K != 0 to something more logical  *)
(* Edit 8/19 to add logic in case epsilon1 was zero.  *)
(* Edit 8/22 to add comments, better error trapping, and improved print 
      statements *)
(* Edit 8/29 to create CodeOutline.txt and eliminate unnecessary prints *)
(* Edit 9/2 to update CodeOutline.txt, add local variables *)
(* Edit 9/3 to finish CodeOutline.txt and do a final cleaning of the code *)
(* Edit 9/4 to proofread over the code one final time & remove any tabs.
      Add additional print statements to 'solitonsoln' *)

(****************************************************************************)

(* debugflags *)

debugHighestOrder = False;
debugHomogenizedEqn = False;
debugSolitonSoln = False;
debugSolitonSolnOper = False;
debugCheckSolns = False;  (* User is asked if they want to verify solutions. *)
(* debugStripper = False; *)

(************************************************************************** *)
(* homogenize[_]: Homogenizes the given PDE by trying a solution of the     *)
(*                form u[1][x, t] = K*D[Log[F[x, t]], {x, n}]               *)
(*                                                                          *)
(* Modules called upon - none                                               *)
(************************************************************************** *)

homogenize[givenPDE_]:=
   Module[{time, lambdaEqn, homEqnwithF, termList, tDerivList, xDerivs, 
           tDerivs, numIntegrations=0, integratedEqn, numDeriv=0, simpIntEqn, 
           fPowersIntEqn, noFPowersIntEqn, kValuesList = {}, kRule, 
           kRuleIntEqn, integratedEqnFinal, homLambda, homLambdaList,
           homOrder, homPDEList = {}, userPDEChoice},
Print["***** Begin 'homogenize' module *****"];
time = TimeUsed[];
If[debugHomogenizedEqn,
   Print["At Pt 1 in the Homogenize module, given is the PDE: "];
   Print[givenPDE];
  ];

(* First, checking to make sure given function is not already homogenized.  *)
lambdaEqn = givenPDE /. {u[1][x, t] -> \[Lambda], 
               Derivative[m_Integer, n_Integer][u[1]][x, t] -> \[Lambda]};
If[FreeQ[lambdaEqn, Plus],
(* then1 *)
   Print["The PDE is already homogenized!"];
   homEqnwithF = givenPDE /. 
                    {u[1] -> F, Derivative[m_,n_][u[1]]->Derivative[m,n][F]};
   If[FreeQ[lambdaEqn, Power],
   (* then2 *)
      Return[{homEqnwithF, 1, k->0, F[x,t]}],
   (* else2 *)
      Return[{homEqnwithF, FactorList[lambdaEqn][[2]][[2]], k->0, F[x,t]}];
     ];
  ];

(* Separates x and t derivatives.                                           *)
termList = givenPDE /. Plus -> List;
tDerivList = Cases[termList, Derivative[0, _][_][x, t]];
tDerivs = tDerivList /. List -> Plus;
xDerivs = Complement[termList, tDerivList] /. List -> Plus;

(* To make the homogenization easier, we will attempt to integrate our PDE  *)
(* with respect to x as many times as possible. We'll continue doing this   *)
(* until we cannot integrate a term NOT containing a t derivative.          *)
If[(FreeQ[Integrate[xDerivs,x], Integrate] === True), 
(* then *)
   While[(FreeQ[Integrate[xDerivs,x], Integrate] === True),
      integratedEqn = Integrate[xDerivs,x] + Integrate[tDerivs,x];
      termList = integratedEqn /. Plus -> List;
      tDerivList = Cases[termList, Derivative[0, _][_][x, t]];
      tDerivs = tDerivList /. List -> Plus;
      xDerivs = Complement[termList, tDerivList] /. List -> Plus;
      numIntegrations++
        ];
   Print["The given PDE can be integrated ", numIntegrations, " time(s) "<>
         "with respect to x. After integration, the PDE becomes: "];
   Print[integratedEqn],
(* else *)
   integratedEqn = givenPDE;
   Print["The PDE cannot be integrated. Thus, the PDE is still: "];
   Print[integratedEqn]
  ];

(* The following while loop will run until a suitable solution is found.    *)
(* start of while *)
While[kValuesList === {},
   ++numDeriv;
(* 5 is arbitrarily chosen as the breaking point in the while loop,         *)
(* feel free to adjust as needed.                                           *)
   If[numDeriv === 5,
      Print[
      "The program was unable to find a valid replacement rule after trying "<>
      "4 different transformations. The code stopped searching and aborted."];
      Abort[];
     ];
   u[1][x_, t_] := K*D[Log[F[x, t]], {x, numDeriv}];
   Print["The code attempts to homogenize the PDE with the following " <>
         "logarithmic-derivative transformation: "];
   Print["u[x,t] = ", u[1][x,t]];
   If[debugHomogenizedEqn,
      Print["At Pt 2 in the Homogenize module, after substitution of " <>
            "u[1][x,t], the PDE becomes: "];
      Print[integratedEqn];
     ];
(* We want the numerator without any K's out front. Below is perhaps        *)
(* overcomplicated way of doing so.                                         *)
   simpIntEqn = Part[Part[FactorList[
                       Numerator[Factor[integratedEqn]]], 3],1];
   If[debugHomogenizedEqn,
      Print["At Pt 3 in the Homogenize module, the PDE simplifies to: "];
      Print[simpIntEqn];
     ];
   fPowersIntEqn = Map[Factor,Collect[simpIntEqn, F[x,t]]];
   If[debugHomogenizedEqn,
      Print["At Pt 4 in the Homogenize module, before substitution of " <>
            "K, the homogenized PDE in powers of F[x,t] is: "];
      Print[fPowersIntEqn];
     ];
   noFPowersIntEqn = Factor[fPowersIntEqn /. {F[x,t] -> 0}];
   If[debugHomogenizedEqn,
      Print["At Pt 5 in the Homogenize module, after letting " <>
            "F(x,t) -> 0, the constraint to be solved for K is: "];
      Print[noFPowersIntEqn];
     ];

(* We want to find a nonzero K value. A zero K value would imply that our   *)
(* PDE is homogeneous, which it cannot be since a homogenization check was  *)
(* performed at the beginning of the module.                                *)
   kValuesList = Flatten[Solve[
                         noFPowersIntEqn == 0 && K != 0, K]];
   If[kValuesList === {},
   (* then *)
      Print["No suitable nonzero values for K were found. "<>
            "The code tries another transformation."], 
   (* else *)
      Print["Possible nonzero K values are: "];
      Print[kValuesList];
     ];
  ]; (* end of while *)

(* Loops over all possible values of K and displays the corresponding       *)
(* homogenized equation.                                                    *)
(* start of do loop *)
Do[
   kRule = Part[kValuesList, i];
   Print["------------------------------------"];
   Print["CHOICE ", i, ": Examining the PDE using ", kRule];
   kRuleIntEqn = Factor[fPowersIntEqn /. kRule];
   If[debugHomogenizedEqn,
      Print["At Pt 6 in the Homogenize module, for K = ", K /. kRule, 
            ",\n the homogenized PDE in factored form: "];
      Print[kRuleIntEqn];
     ];
   integratedEqnFinal = Last[FactorList[kRuleIntEqn]][[1]];
   If[debugHomogenizedEqn,
      Print["At Pt 7 in the Homogenize module, for K = ", K /. kRule,
            ",\n the homogenized PDE in factored form, after common " <> 
            "factors were canceled:"];
      Print[integratedEqnFinal];
     ];

(* Checks to see if PDE is homogenized by replacing powers with lambdas.    *)
(* Similar to check done at the beginning of the module.                    *)
   homLambda = Factor[integratedEqnFinal /. {F[x, t] -> \[Lambda]*F[x, t], 
                  Derivative[m_Integer, n_Integer][F][x, t] -> 
                  \[Lambda]*Derivative[m, n][F][x, t]}];
   homLambdaList = FactorList[homLambda];
   homOrder = Cases[homLambdaList, {\[Lambda], p_Integer}][[1]][[2]];
   If[debugHomogenizedEqn,
      Print["At Pt 8 in the Homogenize module, the code verifies that the " <>
            "PDE is homogenized. Each term involving the symbol" <>
            "F[x,t] is multiplied by a Lambda. Hence, there should be " <>
            "a common number of Lambdas out front in the following eqn: "];
      Print[homLambda];
      Print["At Pt 9 in the Homogenize module, the PDE appears to be " <>
            "homogenized to order: ", homOrder];
     ];
   Print["The homogenized PDE is: "];
   Print[integratedEqnFinal, " = 0"];
   Print["The PDE is homogeneous in the dependent variable "<> 
         "with degree: ", homOrder];

(* This will store the homogenized equation, order, K replacement rule,     *)
(* and u[1][x,t] solution as its own list for each possible K value found.  *)
   homPDEList = Append[homPDEList, {integratedEqnFinal,homOrder,kRule,
                         u[1][x,t]}];
   If[debugHomogenizedEqn,
      Print["At Pt 10 in the Homogenize module, the code returns a part of "<>
            "this list. This list contains the homogenized PDE, K value, "<>
            "and the u[x,t] solution for each choice."];
      Print[homPDEList];
     ],
  {i, 1, Length[kValuesList]}
  ]; (* end of do loop *)

Print["Time used by 'homogenize' module: ", TimeUsed[]-time," seconds"];
Print["***** End 'homogenize' module *****"];
userPDEChoice = 
   Input["Enter your choice (a number) of homogenized equation. "<>
          "Scroll down to see the choices. " <> 
          "If you wish to abort, enter 'qq'"];
If[userPDEChoice === qq,
   Abort[];
  ];
While[userPDEChoice > Length[homPDEList] || userPDEChoice <= 0 || 
      !IntegerQ[userPDEChoice],
   userPDEChoice = Input["The choice you selected does not exist! " <>
                         "Please pick from the choices given below."];
  ];
Return[homPDEList[[userPDEChoice]]];
]; (* closes homogenize module *)

(************************************************************************** *)
(* rootMonomialList[_,_] and monomialList[_,_,_]:                           *)
(*                                                                          *)
(* Modules called upon - none                                               *)
(************************************************************************** *)

rootMonomialList[expr_,funcs_]:=
Cases[{expr},q_/;(Not[And@@(FreeQ[q,#]&/@funcs)]&&Head[q]=!=
      Power&&Head[q]=!=Times&&Head[q]=!=Plus),-1]
If[$VersionNumber>4.,
   monomialList[expr_,rootmonomials_,opts___]:=Module[{s},
   If[$VersionNumber>6.,
   s=GroebnerBasis`DistributedTermsList[expr,rootmonomials,opts],
   s=Internal`DistributedTermsList[expr,rootmonomials,opts]];
(Times@@MapThread[Power,{s[[2]],#}]&)/@(#[[1]]&/@s[[1]])],
   monomialList=MonomialList]

(************************************************************************** *)
(* solitonsoln[_,_]: Finds and returns the soliton solution to a given PDE, *)
(*    or returns zero if one cannot be found.                               *)
(*                                                                          *)
(* Modules called upon - monomialList                                       *)
(************************************************************************** *)

solitonsoln[homPDE_, homOrder_]:=
   Module[{time, epsilonPicker, numExponentials, numPieces, 
           firstEpsTerms, omegaEqn, omegaRule, epsilonLevel,
           (* omegaRuleSol, *) firstZeroOccurance, epsilon, epsilonEval,
           epsilonRHS, solnList, epsilonLHSderivs, (* epsilonPhiList, *)
           epsilonLHS, epsilonFinal, epsilonRHSTerms, epsilonLHSTerms, 
           solnListAll = {}, fSolnTheta},
Print["***** Begin 'solitonsoln' module *****"];
Clear[i];
epsilonPicker = {};
numExponentials = Input["Enter the number of exponentials you wish to " <>
   "have in the soliton solution. ENTER 1, 2, or 3"];
While[numExponentials != 1 && numExponentials != 2 && numExponentials != 3,
   numExponentials = Input["It appears you are either trying to quit," <>
                           "mistyped or misunderstood. " <>
                           "Once again, enter either 1, 2, or 3"];
  ];
time = TimeUsed[];
(* numPieces is arbitrarily chosen number. Very few cases seem to have      *)
(* terms past f[5][x,t].                                                    *)
numPieces = 8;

(* Epsilon serves purely as a bookkeeping variable.                         *)
F[x_,t_] := 1 + Sum[\[Epsilon]^j*f[j][x,t], {j, 1, numPieces}];
Print["Substituting F[x,t]: "];
Print["F[x,t] = ", F[x,t]];
If[debugSolitonSoln,
   Print["At Pt 1 in the Solitonsoln module, after substitution of " <>
         "F[x,t], the PDE becomes: "];
   Print[homPDE];
  ];
epsilonPicker = CoefficientList[homPDE, \[Epsilon]];

(* This may take a long time to print (depending on the size of numPieces), *)
(* however it is nice to see when debugging.                                *)
(*
If[debugSolitonSoln,
   Print["At Pt 2 in the Solitonsoln module, epsilonPicker is: "];
   Do[
      Print["At ",\[Epsilon]^(i-1),":", Part[epsilonPicker, i]],
      {i, 1, Length[epsilonPicker]}
     ];
  ];
*)

(* Based on knowledge that Phi[i][x,t] = E ^ (k[i]*x - omega[i]*t + delta)  *)
phiDerivRules = {Derivative[m_, n_][\[Phi][z_]][x_, t_] -> 
                        k[z]^m*(-\[Omega][z])^n*\[Phi][z][x, t]};
If[debugSolitonSoln,
   Print["At Pt 3 in the Solitonsoln module, the derivative rules " <>
         "for \[Phi] are: "];
   Print[phiDerivRules];
  ];

(* We will always find the Omega to k relationship from the first nonzero   *)
(* Epsilon level. Typically, this is Epsilon 1.                             *)
epsilonLevel = 1;
While[Part[epsilonPicker, epsilonLevel+1] === 0,
   epsilonLevel++;
  ];
If[debugSolitonSoln,
   Print["At Pt 4 in the Solitonsoln module, the first nonzero epsilon " <>
         "occurs at order ",\[Epsilon]^epsilonLevel];
  ];

(******* BEGIN FIRST NONZERO EPSILON LOGIC *******)
firstEpsTerms = Part[epsilonPicker, epsilonLevel+1];
If[debugSolitonSoln,
   Print["------------------------------------------"];
   Print["At Pt 5 in the Solitonsoln module, consider the terms at order ",
         \[Epsilon]^epsilonLevel,": "];
   Print[firstEpsTerms];
  ];
f[1][x_,t_] := Sum[\[Phi][i][x, t], {i, 1, numExponentials}];
Print["Because of the user's choice of exponentials, the code takes "<> 
      "f[1][x,t] as follows: "];
Print[f[1][x,t]];
Print["where \[Phi][i][x,t] = Exp[\[Theta][i]] and "<> 
      "\[Theta][i] = k[i]*x - \[Omega][i]*t + \[Delta][i]"];

(* We assume that all Omega[i]'s will be equal to some relationship         *)
(* involving k[i]'s. Thus, we need to turn all Omega's and k's into         *)
(* generic Omega[i] and k[i], find the relationship, and then transfer      *)
(* them back to their original state.                                       *)
omegaEqn = firstEpsTerms /. phiDerivRules
             /. {\[Omega][_] -> \[Omega][i] , k[_] -> k[i]};
If[debugSolitonSoln,
   Print["At Pt 6 in the Solitonsoln module, setting the k's and " <>
         "\[Omega]'s equal to k[i] and \[Omega][i] in ", 
         \[Epsilon]^epsilonLevel," gives: "];
   Print[omegaEqn];
  ];
omegaRuleSol = Flatten[Solve[omegaEqn == 0, \[Omega][i]]];
omegaRule = {\[Omega][i_] -> omegaRuleSol[[1]][[2]]};
Print["The dispersion relation (relation for k and \[Omega]) is: "];
Print[omegaRule];
count = 1;
Do[
   \[Omega][count] = (\[Omega][i] /. omegaRule) /. {k[i] -> k[count]},
   {count, numExponentials}
  ];
If[debugSolitonSoln,
   Print["At Pt 7 in the Solitonsoln module, verifying the dispersion "<> 
         "relation at order ",\[Epsilon]^epsilonLevel," (should give zero): "];
   Print[Factor[omegaEqn /. omegaRule]];
  ];

(******** BEGIN LOGIC FOR REMAINING EPSILONS ********)
(* The value in firstZeroOccurance will be modified in the loop when the    *)
(* first zero shows up on the RHS. An arbitrary value is initially stored   *)
(* in it to prevent endless looping.                                        *)
arbStoppingPt = 6;
firstZeroOccurance = arbStoppingPt;

(* This is just for ease of reading the code. The value 'i' will represent  *)
(* the level (order) of Epsilon that the loop is on.                        *)
i = epsilonLevel + 1;

(* This stopping was determined by looking at operators (which are          *)
(* implemented in 'solitonsolnOper'). This should be the point were all     *)
(* operators reduce down to zero. BEWARE: This is the stopping point if     *)
(* Epsilon1 is nonzero. I'm not certain this is the correct stopping point  *)
(* in the rare case that Epsilon1 is zero.                                  *)
While[i <= homOrder*(firstZeroOccurance-1),
   epsilon[i] = Factor[Part[epsilonPicker,i+1] /. phiDerivRules];
   If[debugSolitonSoln,
      Print["---------------------------------------------------------"];
      Print["At Pt 8 in the Solitonsoln module, at ",\[Epsilon]^i," :"];
      Print[epsilon[i]];
     ];

(* We are choosing to split up our equation into LHS and RHS to create a    *)
(* more operator-like scheme. The LHS will contain the linear operator on   *)
(* f[i][x,t] & the RHS will contain the rest (f[i-1][x,t] terms and lower). *)
   epsilonRHS[i] = epsilon[i] /. 
                     {f[i][x, t] -> 0, Derivative[p_, q_][f[i]][x, t] -> 0};
   If[debugSolitonSoln,
      Print["At Pt 9 in the Solitonsoln module, the RHS at ",
            \[Epsilon]^i, " is:"];
      Print[epsilonRHS[i]];
     ];

   (* begin if *)
   If[epsilonRHS[i] === 0,
   (* then *)
      Print["At ",\[Epsilon]^i," all terms not involving f[", i, "][x,t] "<>
            "are zero, therefore the code sets f[", i, "][x,t] = 0"];
      f[i][x_,t_] = 0;
      If[debugSolitonSoln,
         Print["At Pt 10 in the Solitonsoln module, at ",
               \[Epsilon]^i," the terms reduce to: "];
         Print[epsilon[i]];
        ];

(* This is where we modify firstZeroOccurance, thus changing the stopping   *)
(* criteria for our loop. Note that the firstZeroOccurance may not          *)
(* represent the first f[_][x,t] that is zero, but it does show the first   *)
(* zero in what we assume will be a run of zero f[_][x,t] terms.            *)
      If[!FreeQ[f[i-1][x,t],\[Phi]],
         firstZeroOccurance = i;
         If[debugSolitonSoln,
            Print["The first zero was found at order: ", firstZeroOccurance];
           ];
        ];
      solnList[i] = {},
   (* else *)
      epsilonLHSderivs[i] = -1*epsilon[i] /. {\[Phi][1][x,t] -> 0, 
                                  \[Phi][2][x,t] -> 0, \[Phi][3][x,t] -> 0};
      If[debugSolitonSoln,
         Print["At Pt 10 in the Solitonsoln module, the LHS at ", 
               \[Epsilon]^i," is: "];
         Print[epsilonLHSderivs[i]];
        ];
      epsilonPhiList[i] = monomialList[epsilonRHS[i], 
                          {\[Phi][1][x,t],\[Phi][2][x,t], \[Phi][3][x,t]}];
      If[debugSolitonSoln,
         Print["At Pt 11 in the Solitonsoln module, the following " <>
               "\[Phi] terms are present on the RHS at ",\[Epsilon]^i," :"];
         Print[epsilonPhiList[i]];
        ];

(* f[i][x,t] will be modeled after what terms were found on the RHS so      *)
(* that the LHS can equal the RHS if our coefficients are chosen correctly. *)
      f[i][x_,t_] = Sum[a[i][k]*Part[epsilonPhiList[i],k], 
                      {k, 1, Length[epsilonPhiList[i]]}]; 
      Print["Because of the terms found at order ",\[Epsilon]^i, 
            ", f[",i,"][x,t] has been selected as: "];
      Print[f[i][x,t]];
      If[debugSolitonSoln,
         Print["At Pt 12 in the Solitonsoln module, after substitution \n" <>
               " of f[", i, "][x,t], the LHS is: "];
         Print[epsilonLHSderivs[i]];
        ];
      epsilonLHS[i] = epsilonLHSderivs[i] /. phiDerivRules;
      If[debugSolitonSoln,
         Print["At Pt 13 in the Solitonsoln module, after application of " <> 
               "the derivative rules for Phi, the LHS is: "];
         Print[epsilonLHS[i]];
        ];
      epsilonFinal[i] = epsilon[i] /. phiDerivRules;

(* This solution list will hold all a[i][j] coefficients at Epsilon^i.      *)
(* These will be found by the following do loop.                            *)
      solnList[i] = {};
      Do[
         If[debugSolitonSoln,
            Print["---------------------------------"];
            Print["At ",\[Epsilon]^i,", solving for a", i, j, " below,"];
           ];
         epsilonRHSTerms[i] = Factor[Coefficient[epsilonRHS[i], 
                              Part[epsilonPhiList[i],j]]];
         epsilonLHSTerms[i] = Factor[Coefficient[epsilonLHS[i], 
                              Part[epsilonPhiList[i],j]]];
         If[debugSolitonSoln,
            Print["At Pt 14 in the Solitonsoln module, the following " <>
                  "Phi term is being considered: "];
            Print[Part[epsilonPhiList[i],j]];
            Print["At Pt 15 in the Solitonsoln module, the LHS " <>
                  "coefficient for this Phi term is: "];
            Print[epsilonLHSTerms[i]];
            Print["At Pt 16 in the Solitonsoln module, the RHS " <>
                  "coefficients for this Phi term is: "];
            Print[epsilonRHSTerms[i]];
           ];
         solnList[i] = Flatten[Append[solnList[i],
          Flatten[Solve[epsilonRHSTerms[i] == epsilonLHSTerms[i], a[i][j]]]]];
         If[debugSolitonSoln,
            Print["At Pt 17 in the Solitonsoln module, the coefficients "<>
                     "\n at ",\[Epsilon]^i," are: "];
            Print[solnList[i]];
           ];
         epsilonFinal[i] = epsilonFinal[i] /. solnList[i],
         {j,1, Length[epsilonPhiList[i]]}
        ]; (* end do loop *)

      Print["The final list of a[",i,"][j] coefficients is: "]
      Print[solnList[i]];

(* If we have found our coefficients correctly, everything should reduce    *)
(* down to zero at this level of Epsilon.                                   *)
      If[debugSolitonSoln,
         Print["At Pt 18 in the Solitonsoln module ",\[Epsilon]^i," is: "];
         Print[Factor[epsilonFinal[i]]];
        ];
      epsilonPicker = epsilonPicker /. solnList[i];
      f[i][x,t] = f[i][x,t] /. solnList[i];
      solnListAll = Append[solnListAll, solnList[i]];
     ]; (* end if epsilonRHS[i] === 0 *)
   i++;
  ]; (* end while *)

Print["---------------------------------------------------------"];

(* Meeting the following condition implies that we were not able to find a  *)
(* run of zero f[_][x,t] and thus unable to compute a soliton solution.     *)
If[firstZeroOccurance === arbStoppingPt,
   Print["The code was unable to produce any zeros for f[i][x,t]. " <>
         "The code stopped searching for a solution and aborted."];
   fSoln[x,t] = 0;
   Print["Time used by 'solitonsoln' module: ", TimeUsed[]-time, " seconds"];
   Print["***** End 'solitonsoln' module *****"],
   (* else *)
   Print["The operator formalism guarantees that all further " <>
         "f[i][x,t] will be zero because the previous f[i][x,t] were zero. " <>
         "Thus, the code has found a solution:"];
   Do[
      f[i][x,t] = 0;
      Derivative[m_, n_][f[i]][x, t] = 0,
      {i,homOrder*(firstZeroOccurance-1)+1, numPieces}
     ];

(* Epsilon serves purely as a bookkeeping variable. Setting it equal to 1   *)
(* will leave the result unchanged.                                         *)
   \[Epsilon] = 1;
   Print[F[x,t]];
   Print["where the a[i][j]'s are: "];
   Print[Flatten[solnListAll]];
   fSolnTheta[x_,t_] = F[x,t] /. Flatten[solnListAll] /. 
            {\[Phi][k_][x,t] -> E^(\[Theta][k])};
   Print["After replacement of the \[Phi]'s, the solution F[x,t] becomes: "];
   Print["F[x,t] = ", fSolnTheta[x,t]];
   Print["where \[Theta][i] = k[i]*x - \[Omega][i]*t + \[Delta][i]"];

   If[debugCheckSolns,
      Print["-----------------------------------------------"];
      Print["Start of verification that F[x,t] is correct."];
      Print["After substitution of F[x,t] into the homogenized PDE the "<>
            "result should equal zero. After substitution one obtains: "];
      Print[Factor[homPDE /. phiDerivRules /. 
         Flatten[Table[solnList[i],{i,2,firstZeroOccurance}]]]];
     ];
   Print["Time used by solitonsoln module: ", TimeUsed[]-time, " seconds"];
   Print["***** End 'solitonsoln' module *****"];
   Return[F[x,t]];
  ];
]; (* closes solitonsoln module *)  

(************************************************************************** *)
(* solitonsolnOper[_,_]: Uses operators to find and return the soliton      *)
(*    solution to a given PDE, or returns zero if one cannot be found.      *)
(*                                                                          *)
(* Modules called upon - monomialList, (* stripper *)                       *)
(************************************************************************** *)

solitonsolnOper[homPDE_, homOrder_]:=
   Module[{numExponentials, time, replacementRulesNN1, replacementRulesNN2, 
          replacementRulesNN3, (* nonlinOp1Sub, nonlinOp2Sub, nonlinOp3Sub, *) 
          powersOfFList, (* linearOp, *) replacementRulesLL, noPowersRules, 
          changeVarRules, nonlinOp1, nonlinOp2, nonlinOp3, newHomPDE, 
          newHomPDEOnlyF, opZeroRules, linearGenRHS, 
          quadraticGenRHS, cubicGenRHS, quarticGenRHS, eps1LHS,
          f1Sub, eps1LHSFDerivs, eps1LHSPhiDerivs, eps1LHSphis, 
    (* omegaRuleSol, *) omegaRule, fSubAll, solnListAll, firstZeroOccurance, 
          epsRHS, epsLHS, epsLHSFDerivs, epsRHSFDerivs, 
          fSolnTheta, fSolnFinal, epsRHSPhiDerivs, epsRHSPhis, 
          epsSolnList, fSub, epsRHSTermsList, epsRHSPhiLists, (* phiTerms, *) 
          coeffList, coeffTermPiece, coeffTerms, coeffTermsAll, phiTermsAll, 
          jDelete, epsLHSFSub, epsLHSPhis, epsRHSTerms, epsLHSTerms},
Clear[i];
numExponentials = Input["Enter the number of exponentials you wish to " <>
   "have in the soliton solution. ENTER 1, 2, or 3"];
While[numExponentials != 1 && numExponentials != 2 && numExponentials != 3,
   numExponentials = Input["It appears you are either trying to quit," <>
                           "mistyped or misunderstood. " <>
                           "Once again, enter either 1, 2, or 3"];
  ];
Print["***** Begin 'solitonsolnOper' module *****"];
time = TimeUsed[];

(* Depending on 'homOrder', these may or may not be modified later on.      *)
replacementRulesNN1 = {};
replacementRulesNN2 = {};
replacementRulesNN3 = {};
nonlinOp1Sub = nonlinOp2Sub = nonlinOp3Sub = 0;

(* All PDEs will contain a linear operator with the same replacement rules. *)
(* The linear operator will always be found at the highest power of F[x,t]. *)
powersOfFList = CoefficientList[homPDE, F[x,t]];
linearOp = Part[powersOfFList, Length[powersOfFList]]; 
Print["The linear operator is: "];
Print[linearOp];
replacementRulesLL = {LL[z_[x,t]] :> (linearOp /. 
          {Derivative[m_, n_][F][x, t] -> Derivative[m, n][z][x,t]})};
If[debugSolitonSolnOper,
   Print["At Pt 1 in the SolitonsolnOper module, the replacement " <>
         "rules for linear operators are:"];
   Print[replacementRulesLL];
  ];

(* Rules for nonlinear operators below.                                     *)
(* noPowersRules will eliminate powers (ranging from 2-4) on derivs by      *)
(* changing the variables to other arbitrary variables. changeVarRules will *)
(* then look for 2 to 4 derivatives multiplied by each other and change the *)
(* variables according to the format seen in Hereman and Nuseir's work.     *)
noPowersRules = {Power[Derivative[m_, n_][a_][x, t], 2] -> 
      Derivative[m, n][a][x, t]*Derivative[m, n][b][x, t], 
   Power[Derivative[m_, n_][a_][x, t], 3] -> 
      Derivative[m, n][a][x, t]*Derivative[m, n][b][x, t]*
      Derivative[m, n][c][x, t],
   Power[Derivative[m_, n_][a_][x, t], 4] -> 
      Derivative[m, n][a][x, t]*Derivative[m, n][b][x, t]*
      Derivative[m, n][c][x, t]*Derivative[m, n][d][x, t]};
changeVarRules = {Times[con___Real, Derivative[a_, b_][__][x, t], 
      Derivative[c_, d_][__][x, t], Derivative[e_, f_][__][x, t], 
      Derivative[g_, h_][__][x, t]] -> 
      Times[con, Derivative[a, b][F][x, t], Derivative[c, d][G][x, t], 
      Derivative[e, f][H][x, t], Derivative[g, h][J][x, t]],
  Times[con___Real, Derivative[a_, b_][__][x, t], 
      Derivative[c_, d_][__][x, t], Derivative[e_, f_][__][x, t]] -> 
      Times[con, Derivative[a, b][F][x, t], Derivative[c, d][G][x, t], 
      Derivative[e, f][H][x, t]],
  Times[con___Real, Derivative[a_, b_][__][x, t], 
      Derivative[c_, d_][__][x, t]] -> 
      Times[con, Derivative[a, b][F][x, t], Derivative[c, d][G][x, t]]};

(* Finding nonlinear operators below.                                       *)
If[Length[powersOfFList] > 1,
(* then1 *)
   nonlinOp1 = Part[powersOfFList, Length[powersOfFList]-1];
   Print["The nonlinear operator 1 (NN1) is: "];
   Print[nonlinOp1];
   nonlinOp1Sub = nonlinOp1 /. noPowersRules /. changeVarRules;
   Print["After substitutions, NN1 is: "];
   Print[nonlinOp1Sub];
   replacementRulesNN1 = {NN1[a_[x,t], b_[x,t]] :> (nonlinOp1Sub /. 
        {Derivative[m_, n_][F][x, t] -> Derivative[m, n][a][x,t],
        Derivative[m_, n_][G][x, t] -> Derivative[m, n][b][x,t]})};
   If[debugSolitonSolnOper,
      Print["At Pt 2 in the SolitonsolnOper module, the replacement " <>
            "rules for the NN1 operator are:"];
      Print[replacementRulesNN1];
     ];
   If[Length[powersOfFList] > 2,
   (* then2 *)
      nonlinOp2 = Part[powersOfFList,Length[powersOfFList]-2];
      Print["The nonlinear operator 2 (NN2) is: "];
      Print[nonlinOp2];
      nonlinOp2Sub = nonlinOp2 /. noPowersRules /. changeVarRules;
      Print["After substitutions, NN2 is: "];
      Print[nonlinOp2Sub];
      replacementRulesNN2 = {NN2[a_[x,t],b_[x,t],c_[x,t]] :> (nonlinOp2Sub /. 
        {Derivative[m_, n_][F][x, t] -> Derivative[m, n][a][x,t],
        Derivative[m_, n_][G][x, t] -> Derivative[m, n][b][x,t],
        Derivative[m_, n_][H][x, t] -> Derivative[m, n][c][x,t]})};
      If[debugSolitonSolnOper,
         Print["At Pt 3 in the SolitonsolnOper module, the replacement " <>
               "rules for the NN2 operator are:"];
         Print[replacementRulesNN2];
        ];
      If[Length[powersOfFList] > 3,
      (* then3 *)
         nonlinOp3 = Part[powersOfFList,Length[powersOfFList]-3];
         Print["The nonlinear operator 3 (NN3) is: "];
         Print[nonlinOp3];
         nonlinOp3Sub = nonlinOp3 /. noPowersRules /. changeVarRules;
         Print["After substitutions, NN3 is: "];
         Print[nonlinOp3Sub];
         replacementRulesNN3 = {NN3[a_[x,t],b_[x,t],c_[x,t],d_[x,t]] :>
           (nonlinOp3Sub /. 
           {Derivative[m_, n_][F][x, t] -> Derivative[m, n][a][x,t],
           Derivative[m_, n_][G][x, t] -> Derivative[m, n][b][x,t],
           Derivative[m_, n_][H][x, t] -> Derivative[m, n][c][x,t],
           Derivative[m_, n_][J][x, t] -> Derivative[m, n][d][x,t]})};
         If[debugSolitonSolnOper,
            Print["At Pt 4 in the SolitonsolnOper module, the replacement " <>
                  "rules for the NN3 operator are:"];
            Print[replacementRulesNN3];
           ];
        ]; (* end If[Length[powersOfFList] > 3] *)
     ]; (* end If[Length[powersOfFList] > 2] *)
  ]; (* end If[Length[powersOfFList] > 1] *)

(* Checks to make sure that if F,G,H,J variables are changed back to all F  *)
(* variables, our PDE is unchanged from what we were originally given.      *)
If[debugSolitonSolnOper,
   Print["----------------------------------------"]
   Print["Start of verification that the new PDE matches the original PDE."];
   newHomPDE = F[x,t]^(homOrder-1)*linearOp+F[x,t]^(homOrder-2)*nonlinOp1Sub+
               F[x,t]^(homOrder-3)*nonlinOp2Sub + nonlinOp3Sub;
   Print["The new PDE, with multiple variables in it, is: "];
   Print[newHomPDE];
   newHomPDEOnlyF = newHomPDE /. {G->F, H->F, J->F};
   Print["After replacing all other variables back to F, the PDE becomes: "];
   Print[newHomPDEOnlyF];
   Print["If subtracting this result from the original PDE produces zero," <>
         " it is likely (although not guaranteed) that the" <>
         " PDE was modified correctly. Subtraction yields: "];
   Print[Factor[newHomPDEOnlyF - homPDE]];
  ];
opZeroRules = {LL[0] -> 0, NN1[0, __] -> 0, NN1[__, 0] -> 0, 
               NN2[0, __, __] -> 0, NN2[__, 0, __] -> 0, NN2[__, __, 0] -> 0, 
               NN3[0, __, __, __] -> 0, NN3[__, 0, __, __] -> 0, 
               NN3[__, __, 0, __] -> 0, NN3[__, __, __, 0] -> 0};

(* Based on knowledge that phi[i][x,t] = E ^ (k[i]*x - omega[i]*t + delta)  *)
phiDerivRules = {Derivative[m_, n_][\[Phi][z_]][x_, t_] -> 
                    k[z]^m*(-\[Omega][z])^n*\[Phi][z][x, t]};

(* General formulas developed to determine what operators will be present   *)
(* on the RHS at any level of n. Formulas vary based on order of given PDE. *)
linearGenRHS[n_] := -Sum[f[n - j][x, t]*LL[f[j][x, t]], {j, 2, n - 1}];
quadraticGenRHS[n_] := -Sum[f[n - j][x, t]*LL[f[j][x, t]], 
      {j, 2, n - 1}] - Sum[NN1[f[n - j][x, t], f[j][x, t]], {j, 1, n - 1}];
cubicGenRHS[n_] := -(Sum[2*f[n-j][x, t]*LL[f[j][x, t]], {j, 2, n-1}]+ 
   Sum[NN1[f[n-j][x, t], f[j][x, t]], {j, 1, n-1}]+ 
   Sum[Sum[f[n-m][x, t]*f[m - j][x, t]*LL[f[j][x, t]], 
      {j, 2, m-1}], {m, 2, n-1}]+ 
   Sum[Sum[f[n-m][x, t]*NN1[f[m-j][x, t], f[j][x, t]]+ 
     NN2[f[n-m][x, t], f[m-j][x, t], f[j][x, t]], {j, 1, m-1}], {m, 2, n-1}]);
quarticGenRHS[n_] := -(Sum[3*f[n-j][x, t]*LL[f[j][x, t]], {j, 2, n-1}]+ 
    Sum[NN1[f[n-j][x, t], f[j][x, t]], {j, 1, n-1}]+ 
    Sum[Sum[3*f[n-m][x, t]*f[m-j][x, t]*LL[f[j][x, t]], 
       {j, 2, m-1}], {m, 2, n-1}]+ 
    Sum[Sum[2*f[n-m][x, t]*NN1[f[m-j][x, t], f[j][x, t]]+NN2[f[n-m][x, t], 
       f[m-j][x, t], f[j][x, t]], {j, 1, m-1}], {m, 2, n-1}]+ 
    Sum[Sum[Sum[f[n-l][x, t]*f[l-m][x, t]*f[m-j][x, t]*
       LL[f[j][x, t]], {j, 2, m-1}], {m, 2, l-1}], {l, 3, n-1}]+ 
    Sum[Sum[Sum[
       f[n-l][x, t]*f[l-m][x, t]*NN1[f[m-j][x, t], f[j][x, t]]+
       f[n-l][x, t]*NN2[f[l-m][x, t], f[m-j][x, t], f[j][x, t]]+ 
       NN3[f[n-l][x, t], f[l-m][x, t], f[m-j][x, t], 
       f[j][x, t]], {j, 1, m-1}], {m, 2, l-1}], {l, 3, n-1}]);

Print["----------------------------------------"];
Print["The code seeks a solution of the form F[x,t] = 1 + ", 
       Sum[f[i][x,t], List[i,1, DirectedInfinity[1]]]];
Print["To find the solution, the code will compute each f[i][x,t] " <>
      "separately, using the schemes below."];
If[Length[powersOfFList] === 1,
(* then1 *)
   Print["Since the PDE contains only a linear operator, the code will use " <>
         "the following scheme to derive the RHS, where 'n' is " <>
         "chosen to find a given f[n][x,t]."];
   Print[linearGenRHS[n]],
(* else1 *)
   If[Length[powersOfFList] === 2,
   (* then2 *)
      Print["Since the PDE contains 2 operators, the code uses the " <>
            "following scheme to derive the RHS, where 'n' is chosen " <>
            "to find a given f[n][x,t]."];
      Print[quadraticGenRHS[n]],
   (* else2 *)
      If[Length[powersOfFList] === 3,
      (* then3 *)
         Print["Since our PDE contains 3 operators, the code uses " <>
               "the following scheme to derive the RHS, where 'n' is " <>
               "chosen to find a given f[n][x,t]."];
         Print[cubicGenRHS[n]],
      (* else3 *)
         If[Length[powersOfFList] === 4,
         (* then4 *)
            Print["Since our PDE contains 4 operators, the code uses " <>
                  "the following scheme to derive the RHS, where 'n' is " <>
                  "chosen to find a given f[n][x,t]."];
            Print[quarticGenRHS[n]],
         (* else4 *)
            Return[0];
           ];
        ];
     ];
  ];
Print["The scheme to derive the LHS is as follows: "];
Print[LL[f[n][x,t]]];

(**************** BEGIN f[1][x,t] LOGIC ********************)
Print["------------------------------------"];
Print["The code attempts to find f[1][x,t]."];
eps1LHS = LL[f[1][x,t]];
Print["According to the scheme, the LHS is: ", eps1LHS];
Print["According to our scheme, the RHS is: ", 0];

(* The true value of f[1][x,t] will be held in ff[1][x,t]. This is done     *)
(* because we do not wish to substitute f[i][x,t]'s in future RHS           *)
(* equations. Thus, our RHS equations will look cleaner.                    *)
ff[1][x_,t_] = Sum[\[Phi][i][x, t], {i, 1, numExponentials}];
f1Sub = {Derivative[m_,n_][f[1]][x,t] -> Derivative[m,n][ff[1]][x,t],
          f[1][x,t] -> ff[1][x,t]};
Print["As determined by the choice of exponentials. f[1][x,t] is: "];
Print[ff[1][x,t]];
Print["where \[Phi][i][x,t] represents the exponential E^(k[i]*x-\[Omega]" <>
      "[i]*t + \[Delta]). Thus, the Phi derivative rules are as follows: "];
Print[phiDerivRules];

eps1LHSFDerivs = eps1LHS /. replacementRulesLL;
If[debugSolitonSolnOper,
   Print["At Pt 5 in the SolitonsolnOper module, after the linear " <>
         "operator is applied, the LHS is: "];
   Print[eps1LHSFDerivs];
  ];
eps1LHSPhiDerivs = eps1LHSFDerivs /. f1Sub;
If[debugSolitonSolnOper,
   Print["At Pt 6 in the SolitonsolnOper module, the LHS with f[1][x,t] " <>
         "substituted is: "];
   Print[eps1LHSPhiDerivs];
  ];
eps1LHSphis = eps1LHSPhiDerivs /. phiDerivRules;
If[debugSolitonSolnOper,
  Print["At Pt 7 in the SolitonsolnOper module, the LHS with " <>
        "phiDerivRules applied is: "];
  Print[eps1LHSphis];
  ];
eps1LHS = eps1LHSphis /. {\[Omega][_] -> \[Omega][i] , k[_] -> k[i]};
If[debugSolitonSolnOper,
   Print["At Pt 8 in the SolitonsolnOper module, the LHS with generic " <>
         "k's and \[Omega]'s is: "];
   Print[eps1LHS];
  ];
omegaRuleSol = Flatten[Solve[eps1LHS == 0, \[Omega][i]]];
omegaRule = {\[Omega][i_] -> omegaRuleSol[[1]][[2]]};
Print["In order for the LHS = RHS, the replacement rule for the " <>
      "\[Omega]'s must be: "];
Print[omegaRule];

(* The following do loop ensures that any future Omega[i]'s present will    *)
(* automatically be changed into their respective k[i] values. Omega's      *)
(* should not appear in terms after this point.                             *)
Do[
   \[Omega][j] = (\[Omega][i] /. omegaRule) /. {k[i] -> k[j]},
   {j, 1, numExponentials}
  ];
If[debugSolitonSolnOper,
   Print["At Pt 9 in the SolitonsolnOper module, the LHS should reduce to " <>
         "zero. The LHS is: "];
   Print[Factor[eps1LHS /. omegaRule]];
  ];

(**************** BEGIN f[2][x,t] LOGIC AND ONWARD ********************)
(* These are a list of lists which will contain every nonzero value of      *)
(* f[i][x,t] in fSubAll and all a[i][j] values in solnListAll.              *)
fSubAll = {f1Sub};
solnListAll = {};

(* The value in firstZeroOccurance represents an arbitrary stopping point.  *)
(* Unlike in 'solitonsoln' module, this value will NOT be modified within   *)
(* the while loop. Please feel free to change this according to your needs. *)
firstZeroOccurance = 6;

(* The following while loop will run until all operators terms on the RHS   *)
(* have a zero present in them, or until the stopping point is reached.     *)
i = 2;
While[i <= homOrder*(firstZeroOccurance-1),
   Print["-------------------------------------"];
   Print["The code attempts to find f[", i, "][x,t]."];
   If[Length[powersOfFList] === 1,
      epsRHS[i] = linearGenRHS[i],
   (* else1 *)
      If[Length[powersOfFList] === 2,
         epsRHS[i] = quadraticGenRHS[i],
      (* else2 *)
         If[Length[powersOfFList] === 3,
            epsRHS[i] = cubicGenRHS[i],
         (* else3 *)
            If[Length[powersOfFList] === 4,
               epsRHS[i] = quarticGenRHS[i],
            (* else4 *)
                Return[0];
              ];
           ];
        ];
     ]; (* end If[homOrder === 1] *)
   epsLHS[i] = LL[f[i][x,t]];
   If[debugSolitonSolnOper,
      Print["The LHS, from the above scheme, is: "];
      Print[epsLHS[i]];
      Print["The RHS, from the above scheme, is: "];
      Print[epsRHS[i]];
     ];
   epsLHSFDerivs[i] = epsLHS[i] /. replacementRulesLL;
   epsRHSFDerivs[i] = epsRHS[i] /. replacementRulesLL /. 
                        replacementRulesNN1 /. replacementRulesNN2 /. 
                        replacementRulesNN3 /. opZeroRules;
   If[debugSolitonSolnOper,
      Print["At Pt 10 in the SolitonsolnOper module, LHS at level ", i, 
            "\n with the linear operator applied is: "];
      Print[epsLHSFDerivs[i]];
      Print["At Pt 11 in the SolitonsolnOper module, RHS at level ", i, 
            "\n with operators applied is: "];
      Print[epsRHSFDerivs[i]];
     ];

(* The following condition will only be met if all RHS operators contain a  *)
(* zero from a previously zero f[i][x,t] term. If satisfied, the module     *)
(* will return the found F[x,t] solution in terms of E ^ theta[i].          *)
   If[epsRHSFDerivs[i] === 0,
   (* then *)
      Print["The code has reached the point where all other future f[i] " <>
         "terms must be zero because a zero is present in each operator."];
      Print["-------------------------------------"];
      Print["Let F[x,t] be the sum of all f[i][x,t] + 1. " <>
            "Thus, F[x,t] is: "];
      FF[x_,t_] = 1 + Sum[ff[j][x,t], {j,1,i-1}];
      Print[FF[x,t]];
      Print["where our a[i][j]'s are: "];
      Print[Flatten[solnListAll]];
      Print["Substituting the coefficients a[i][j] into F[x,t] yields: "];
      Print[FF[x,t] /. Flatten[solnListAll]];
      fSolnTheta[x_,t_] = FF[x,t] /. Flatten[solnListAll] /. 
               {\[Phi][k_][x,t] -> E^(\[Theta][k])};
      Print["Without the \[Phi]'s, the solution F[x,t] is: "];
      Print["F[x,t] = ", fSolnTheta[x,t]];
      Print["where \[Theta][i] = k[i]*x - \[Omega][i]*t + \[Delta][i]"];
      If[debugCheckSolns,
         Print["-----------------------------------------------"];
         Print["Start of verification that F[x,t] is correct."];
         Print["The original homogenized PDE was: "];
         Print[homPDE];
         fSolnFinal[x_,t_] = fSolnTheta[x,t] /. 
            {\[Theta][a_] -> k[a]*x - \[Omega][a]*t + \[Delta][a]};
         Print["The F[x,t], with \[Theta] replaced, is: "];
         Print[fSolnFinal[x,t]];
         Print["After substituting F[x,t], the PDE should reduce to zero. " <>
               "The PDE is now: "];
         Print[Factor[homPDE /. {F[x,t] -> fSolnFinal[x,t], 
            Derivative[m_,n_][F][x,t] -> Derivative[m,n][fSolnFinal][x,t]}]];
        ];
      Print["Time used by solitonsolnOper module: ", TimeUsed[]-time, " sec"];
      Print["***** End 'solitonsolnOper' module *****"];
      Return[fSolnTheta[x,t]],

   (* else *)
      epsRHSPhiDerivs[i] = epsRHSFDerivs[i] /. Flatten[fSubAll];
      If[debugSolitonSolnOper,
         Print["At Pt 12 in the SolitonsolnOper module, RHS with " <>
               "previous f's substituted is: "];
         Print[epsRHSPhiDerivs[i]];
        ];
(* 'Expand' is used so that we can find each individual Phi term multiplied *)
(* by our coefficients. 'Factor' will also work, but is often slower.       *)
      epsRHSPhis[i] = Expand[epsRHSPhiDerivs[i] /. phiDerivRules /. 
                         Flatten[solnListAll], \[Phi]]; 
      If[debugSolitonSolnOper,
         Print["At Pt 13 in the SolitonsolnOper module, RHS with " <>
               "no Phi derivatives applied is: "];
         Print[epsRHSPhis[i]];
        ];
     ]; (* end If[epsRHSFDerivs[i] === 0] *)

(* A zero produced here on the RHS does NOT imply the same thing as a zero  *) 
(* produced above. A zero here is not enough to cause termination of our    *)
(* while loop, however it does allow us to move to the next iteration.      *)
   If[epsRHSPhis[i] === 0,
   (* then *)
      f[i][x_,t_] = 0;
      ff[i][x_,t_] = 0;
      epsSolnList[i] = {};
      fSub[i] = {};
      Print["Because the substitution of earlier f[i][x,t]'s results in " <>
            "the RHS being zero, the code sets f[", i, "][x,t] = 0"],
   (* else *)
      If[FreeQ[epsRHSPhis[i], Plus],
         epsRHSTermsList[i] = {epsRHSPhis[i]},
      (* else *)
        epsRHSTermsList[i] = Apply[List, epsRHSPhis[i]]; 
        ];
      If[debugSolitonSolnOper,
         Print["At Pt 14 in the SolitonsolnOper module, RHS with terms " <>
               "in lists is: "];
         Print[epsRHSTermsList[i]];
        ];
      Do[
         epsRHSPhiLists[i][j] = monomialList[Part[epsRHSTermsList[i], j], 
            {\[Phi][1][x, t], \[Phi][2][x, t], \[Phi][3][x, t]}],
         {j, 1, Length[epsRHSTermsList[i]]}
        ];
      phiTerms[i] = Flatten[Union[Table[epsRHSPhiLists[i][j], 
         {j, 1, Length[epsRHSTermsList[i]]}]]];
      If[debugSolitonSolnOper,
         Print["At Pt 15 in the SolitonsolnOper module, the following " <>
               "Phi terms were found on the RHS: "];
         Print[phiTerms[i]];
        ];
      Table[coeffList[i][j] = {}, {j, 1, Length[phiTerms[i]]}];

(* The following inner do loop will run through each part of the list of    *)
(* RHS terms, checking to see if it matches the Phi term controlled by the  *)
(* outer do loop. If a match is found, the coefficient is put into a list   *)
(* containing other coefficients of the same Phi term.                      *)
      Do[
         Do[
         If[Part[epsRHSPhiLists[i][m], 1] === Part[phiTerms[i], j],
               coeffList[i][j] = Append[coeffList[i][j], 
                     Coefficient[Part[epsRHSTermsList[i], m], 
                     epsRHSPhiLists[i][m]]];
               ],
            {m, 1, Length[epsRHSTermsList[i]]}
           ],
         {j, 1, Length[phiTerms[i]]}
        ];

(* The lists of coefficients are then replaced back into sums & factored.   *)
(* The factoring is done in hopes of some terms reducing down to zero.      *)
      Do[
         coeffTermPiece[i][j] = Factor[Flatten[coeffList[i][j]] /. 
                                  {List -> Plus}],
        {j, 1, Length[phiTerms[i]]}
        ];
      coeffTerms[i] = Flatten[Table[coeffTermPiece[i][j], 
                        {j, 1, Length[phiTerms[i]]}]];

(* These two are created solely for the purpose of use in the next do loop. *)
      coeffTermsAll[i] = coeffTerms[i];
      phiTermsAll[i] = phiTerms[i];

(* jDelete keeps a counter of what terms in coeffTerms and phiTerms can be  *)
(* deleted because the coefficient is zero. The do loop performs the        *)
(* appropriate deletions.                                                   *)
      jDelete = 1;
      Do[
         If[debugSolitonSolnOper,
            Print["At Pt 16 in the SolitonsolnOper module, the " <>
                  "coefficient for ", Part[phiTermsAll[i],j], " is: "];
            Print[Part[coeffTermsAll[i], j]];
           ];
         If[Part[coeffTermsAll[i], j] === 0,
            coeffTerms[i] = Delete[coeffTerms[i], jDelete];
            phiTerms[i] = Delete[phiTerms[i], jDelete];
            jDelete--;
           ];
         jDelete++,
        {j, 1, Length[phiTermsAll[i]]}
        ];

      If[Length[phiTerms[i]] === 0,
         f[i][x_,t_] = 0;
         ff[i][x_,t_] = 0;
         epsSolnList[i] = {};
         fSub[i] = {};
         Print["Because all coefficients are zero, the code sets f[", 
               i, "][x,t] = 0"],

      (* else *)
(* If terms are still present, then we can go on with constructing our      *)
(* f[i][x,t] solution. As mentioned before, the true value of f[i][x,t]     *)
(* will be held in ff[i][x,t] because we do not wish to substitute our      *)
(* f[_][x,t]'s in future RHS equations.                                     *)
         ff[i][x_,t_] = Sum[a[i][j]*Part[phiTerms[i],j], 
                          {j, 1, Length[phiTerms[i]]}];
         fSub[i] = {Derivative[m_,n_][f[i]][x,t]->Derivative[m,n][ff[i]][x,t],
                      f[i][x,t] -> ff[i][x,t]};
         Print["Because the terms found on the RHS, f[", i, "][x,t] " <>
               "has been chosen as: "];
         Print[ff[i][x,t]];
         epsLHSFSub[i] = epsLHSFDerivs[i] /. fSub[i];
         If[debugSolitonSolnOper,
            Print["At Pt 17 in the SolitonsolnOper module, the LHS \n" <>
                  " with f[", i, "][x,t] substituted is: "];
            Print[epsLHSFSub[i]];
           ];
         epsLHSPhis[i] = epsLHSFSub[i] /. phiDerivRules;
         If[debugSolitonSolnOper,
            Print["At Pt 18 in the SolitonsolnOper module, the LHS " <>
                  "without Phi derivatives is: "];
            Print[epsLHSPhis[i]];
           ];

(* This solution list will hold all a[i][j]'s at level 'i'. These will be   *)
(* found by the following do loop.                                          *)
         epsSolnList[i] = {};
         Do[
            If[debugSolitonSolnOper,
               Print["---------------------------------"];
               Print["At level ", i ,", solving for a[", i, "][",
                     j, "] below,"];
              ];
            epsRHSTerms[i] = Part[coeffTerms[i], j];
            epsLHSTerms[i] = Coefficient[epsLHSPhis[i], 
                                Part[phiTerms[i],j]];
            If[debugSolitonSolnOper,
               Print["At Pt 19 in the SolitonsolnOper module, the Phi " <>
                     "term \n multiplied by a[", i, "][", j, "] is: "];
               Print[Part[phiTerms[i],j]];
               Print["At Pt 20 in the SolitonsolnOper module, the LHS " <>
                     "coefficient for this Phi term is: "];
               Print[epsLHSTerms[i]];
               Print["At Pt 21 in the SolitonsolnOper module, the RHS " <>
                     "coefficient for this Phi term is: "];
               Print[epsRHSTerms[i]];
              ];
            epsSolnList[i] = Flatten[Append[epsSolnList[i],
                     Factor[Flatten[Solve[epsRHSTerms[i] == epsLHSTerms[i], 
                     a[i][j]]]]]];
            If[debugSolitonSolnOper,
               Print["At Pt 22 in the SolitonsolnOper module, the " <>
                     "coefficients \n at level ", i, " are: "];
               Print[epsSolnList[i]];
              ],
           {j,1, Length[phiTerms[i]]}
           ]; (* end do loop *)
         Print["To obtain LHS = RHS, the a[", i, "][j] terms must be: "];
         Print[epsSolnList[i]];
        ]; (* end If[Length[phiTerms] === 0] *)
     ]; (* end If[epsRHSPhis[i] === 0] *)
   solnListAll = Append[solnListAll, epsSolnList[i]];
   fSubAll = Append[fSubAll, fSub[i]];
   i++;
  ]; (* end while *)

Print["-------------------------------------"];
Print["The code was unable to find a stopping point guaranteed by the "<> 
      "operator scheme. The code terminated and returned nothing."];
Print["Time used by solitonsolnOper module: ", TimeUsed[]-time, " seconds"];
Print["***** End 'solitonsolnOper' module *****"];
]; (* closes solitonsolnOper module *)

(************************************************************************** *)
(* printpage[ , ]: a subroutine for menu                                    *)
(*                                                                          *)
(* Modules called upon - none                                               *)
(************************************************************************** *)
 
printpage[n_,page_] := Module[{i,choice,lenpage},
lenpage = Length[page[n]];

Print["  *** MENU INTERFACE ***  (page: ",n,")"];
Print["---------------------------------------------------"];
For[i=1,i <= lenpage,i++,
   Print[Part[Part[page[n],i],1]]
   ];
Print[" (nn) Next Page"];
Print[" (tt) Your System"];
Print[" (qq) Exit the Program"];
Print["---------------------------------------------------"];
choice = Input["ENTER YOUR CHOICE: "];
Return[choice] 
]; (* end printpage module *)

(************************************************************************** *)
(* menu: creates the menu                                                   *)
(*                                                                          *)
(* Modules called upon - printpage                                          *)
(************************************************************************** *)

menu := Module[
{counterpage = 1,menulist,numpages,page,
choice1,control,choice2,choice3,lenmenulist,i},
menulist = 
{
{"  (1) KdV Equation (h_kdv.m)"},
{"  (2) Modified KdV Equation (h_mkdv.m)"},
{"  (3) Burgers Equation (h_burg.m)"},
{"  (4) Kupershmidt Equation (h_kuper.m)"},
(* 
{"  (xxxxx) Boussinesq System (h_bous.m)"},
*)
{"  (5) Fifth Order KdV Equation-KK Case (h_5kk.m)"},
{"  (6) Fifth Order KdV Equation-Ito Case (h_5ito.m)"},
{"  (7) Fifth Order KdV Equation (h_pa5kdv.m)"},
{"  (8) Fifth Order Sawada-Kotera KdV Equation (h_skkdv.m)"},
{"  (9) Fifth Order Lax KdV Equation (h_laxkdv.m)"},
{" (10) Fifth Order Kaup-Kupershmidt KdV Equation (h_kkkdv.m)"},
{" (11) Fitzhugh-Nagumo Equation with Convection (h_fitz.m)"},
{" (12) Fisher Equation with Convection (h_fish.m)"},
{" (13) An Equation due to Calogero (h_cal.m)"},
{" (14) An Equation due to Calogero, with variable change (h_cal2.m)"}
(*
{" (15) Ito System-parameterized (d_pito.m)"},
{" (16) Ito System (d_ito.m)"},
{" (17) NLS System (real and imaginary parts) (d_nls.m)"},
{" (18) NLS System (equation and conjugate) (d_nlsalt.m)"},
{" (19) DNLS System (d_dnls.m)"},
{" (20) MVDNLS System (d_mvdnls.m)"},
{" (21) Kaup System-parameterized (d_pkaup.m)"},
{" (22) Kaup System (d_kaup.m)"},
{" (23) Kaup-Broer System (d_broer.m)"},
{" (24) Drinfel'd-Sokolov System (d_soko.m)"},
{" (25) Dispersive Long Wave System (d_disper.m)"}, 
{" (26) Non-dispersive Long Wave System (d_nodisp.m)"}, 
{" (27) 3-Component KdV System (d_3ckdv.m)"},
{" (28) 2-Component Nonlinear Schrodinger Equation (d_2cnls.m)"},
{" (29) Boussinesq System (d_bous.m)"}
*)
}; (* closes menulist *)

lenmenulist = Length[menulist];	
numpages = Ceiling[lenmenulist/10];
For[i = 1,i <= numpages,i++,
   page[i] = If[lenmenulist >= (i*10),
      menulist[[Table[k,{k,(i-1)*10+1,i*10}]]],
      menulist[[Table[k,{k,(i-1)*10+1,lenmenulist}]]]]
   ];

(* Call to printpage module *)
choice1 = printpage[counterpage,page];

control := (
Switch[choice1,
tt,Print["Make sure that you have prepared the data file for the system"];
   Print["you want to test (similar to the data files supplied)."];
   choice2 = Input["If your data file is ready, press 1, else 2: "];
   If[choice2 === 1,
      choice3 = InputString["Enter the name of your data file: "];
      Get[choice3],
      Print["Aborting the computations!"];
      Print["Prepare your data file first, then start over."];
      Abort[]
     ],
nn,If[counterpage < numpages,counterpage++;
      choice1 = printpage[counterpage,page]; control,
      counterpage = 1; choice1 = printpage[1,page]; control],
qq,Print["Aborting the computations!"];Abort[],
1,<<h_kdv.m,2,<<h_mkdv.m,3,<<h_burg.m,4,<<h_kuper.m,
(* 5,<<h_bous.m, *)
5,<<h_5kk.m,6,<<h_5ito.m,7,<<h_pa5kdv.m,8,<<h_skkdv.m,9,<<h_laxkdv.m,
10,<<h_kkkdv.m,11,<<h_fitz.m,12,<<h_fish.m,13,<<h_cal.m,14,<<h_cal2.m,_,
(* 
15,<<d_pito.m,16,<<d_ito.m,17,<<d_nls.m,19,<<d_nlsalt.m,              
*)
   Print["Aborting the computations!"];
   Abort[]
  ]; (* closes Switch *)
); (* end control *)
control;
(*
If[Not[StringQ[myfile]],
   myfile = InputString["Enter the name of the output file: "]]	
*)
]; (* end menu module *)

(****************************************************************************)
(* pdeform[_]: takes an expression & prints it with subscripted derivatives *)
(*            and removes the functional dependence (i.e. f[x,t] -> f)      *)
(*            Written by Tracy Otto & Tony Miller (CSM-1995).               *)
(*            Modified by U. Goktas.                                        *)
(*                                                                          *)
(* Modules called upon - none                                               *)
(************************************************************************** *)

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 *)

(****************************************************************************)
(* length[_]: a modified length module                                      *)
(*                                                                          *)
(* Modules called upon - none                                               *)
(************************************************************************** *)

length[a_] := Module[{l},
If[Head[a] === Plus,l = Length[a],
   If[Head[a] === Times, l = Length[a], l = 1 ] 
  ]; 
Return[l] 
]; (* end length module *)

(****************************************************************************)
(* highestorder[_]: returns the highest order of a PDE                      *)
(*                                                                          *)
(* Modules called upon - length                                             *)
(************************************************************************** *)
 
highestorder[gexpr_] := Module[
{numTerms,i,currentTerm,lenCurrentTerm,derivs,numXDerivs,numTDerivs,term,k,
partOfTerm,newexpr3,totalDerivs,newexpr4,orderlist1 ={},orderlist2 = {},
orderlist3 = {},highestlist1,highestlist2,highestlist3,expire,derivList = {}, 
highOrder = 0},

derivList = Append[derivList, {"X Derivs", "T Derivs", "Max Derivs", "Term"}];
expr = Expand[gexpr];

If[Not[FreeQ[expr, Derivative]], 
   numTerms = length[expr];
   numXDerivs = {}; numTDerivs = {}; totalDerivs = {}; term = 1;
   For[i = 1, i <= numTerms, i++,
      If[Not[FreeQ[expr, Plus]],
         numXDerivs = {};
         numTDerivs = {};
         totalDerivs = {};
         term = 1
        ];
      If[numTerms != 1, currentTerm = expr[[i]], currentTerm = expr];
      If[Not[FreeQ[currentTerm, Derivative]],
         lenCurrentTerm = length[currentTerm];
         For[k = 1, k <= lenCurrentTerm, k++,
	    If[lenCurrentTerm != 1,
               partOfTerm = currentTerm[[k]],
               partOfTerm = currentTerm
              ];
            If[Not[FreeQ[partOfTerm,Derivative]],       
               If[FreeQ[partOfTerm, Power], 
                  newexpr3 = FullForm[partOfTerm];
                  derivs = newexpr3[[1]][[0]][[0]];
                  numXDerivs = Append[numXDerivs,derivs[[1]]];
                  numTDerivs = Append[numTDerivs,derivs[[2]]];
                  totalDerivs = Append[totalDerivs,derivs[[1]]+derivs[[2]]];
                  newexpr4 = partOfTerm;
                  term = term*partOfTerm,

                  newexpr3 = FullForm[partOfTerm];
                  derivs = newexpr3[[1]][[1]][[0]][[0]]; 
                  numXDerivs = Append[numXDerivs,derivs[[1]]]; 
                  numTDerivs = Append[numTDerivs,derivs[[2]]];
                  totalDerivs = Append[totalDerivs,derivs[[1]]+derivs[[2]]];
                  newexpr4 = newexpr3[[1]];
                  term = term*newexpr3[[1]];
                 ];      (* Matches If[FreeQ[partOfTerm, Power]] *)
               ,term = partOfTerm
              ]; (* Matches If[Not[FreeQ[newexp2,Deriv]]] *)
           ];    (* Matches For[k = 1, k <= lenCurrentTerm, k++] statement *)

         numXDerivs = Max[numXDerivs];
         numTDerivs = Max[numTDerivs];
         totalDerivs = Max[totalDerivs];
         If[ totalDerivs > highOrder, highOrder = totalDerivs];
         derivList = Append[derivList, {numXDerivs, numTDerivs, totalDerivs, 
            term}];
        ];  (* Closes If[Not[FreeQ[currentTerm,Derivative]]] *)
     ];     (* Closes For statement *)   
  ];        (* Closes If[Not[FreeQ[expr,Derivative]]] *)

If[debugHighestOrder,
   Print["At Pt 1 in module highestorder, the table of derivatives is: "];
   Print[TableForm[derivList]];
];
Return[highOrder]
]; (* closes highestorder module *)

(**********************************************************************)
(* stripper[givenexpr, 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 weighted parameters if there are any            *)
(* Output:  A list of the terms of the expression with all constant   *)
(*          coefficients removed.                                     *)
(* Adapted From:  data/newzeala/reu2004/reu2004-0811/software/        *)
(*                whlstrip.m                                          *)
(* Code is in File:  conservationlaws/dpnstrip.m                      *)
(* Major Changes Made: 22 May, 2006, 14:30 by DP at CSM               *)
(* Last Modified:  29 September,?009, 17:39 by DP at home            *)
(**********************************************************************)

(* Clear[stripper]; *)

stripper[givenexpr_, paramswithweight_:{}] :=
Module[{strippedTermsList, expandExpr, expnList, termsNestedInFunctions,
    rules1, rules2, rules3, coefficientList, singleNoCoefficientList,
    expnFnList, possibleIndVarCombos, possibleDepVarCombos, permutDepVarCombos,
    possibleVarCombos, i1, n1, printSTRIP},

  If[debugStripper, printSTRIP = Print, Clear[printSTRIP], Clear[printSTRIP]];

  printSTRIP["debug code: STRIP, Function: stripper, File: dpnstrip.m"];

  printSTRIP["At STRIP IN, the expression given to the stripper:"];
  printSTRIP[givenexpr];
  printSTRIP["At STRIP IN, parameters given to the stripper:"];
  printSTRIP[paramswithweight];

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

  If[Not[MemberQ[givenexpr, u[_][__], {0, Infinity}]] &&
      Not[MemberQ[givenexpr, Derivative[__][u[_]][__], {0, Infinity}]],
    If[Not[MemberQ[givenexpr, n1_ /;
           MemberQ[{x, y, z, t}, n1], {0, Infinity}]],
       Return[{1}]
      ] (* end If Not MemberQ *)
    ]; (* end If Not MemberQ *)

  (* Input expression must be in expanded form.                         *)
  expandExpr = Expand[givenexpr];
  printSTRIP["At STRIP 1, the expression given to the stripper in"<>
      " expanded form:"];
  printSTRIP[expandExpr];

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

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

  (* Find all possible combinations of independent variables and        *)
  (* dependent variables that could be nested inside functions.  This   *)
  (* includes u[_][__] and derivatives of u.                            *)
  possibleIndVarCombos = {{t}, {x}, {y}, {z}, {t, x}, {t, y}, {t, z}, {x, t},
      {x, y}, {x, z}, {y, t}, {y, x}, {y, z}, {z, t}, {z, x}, {z, y},
      {t, x, y}, {t, x, z}, {t, y, x}, {t, y, z}, {t, z, x}, {t, z, y},
      {x, t, y}, {x, t, z}, {x, y, t}, {x, y, z}, {x, z, t}, {x, z, y},
      {y, t, x}, {y, t, z}, {y, x, t}, {y, x, z}, {y, z, t}, {y, z, x},
      {z, t, x}, {z, t, y}, {z, x, t}, {z, x, y}, {z, y, t}, {z, y, x},
      {t, x, y, z}, {t, x, z, y}, {t, y, x, z}, {t, y, z, x}, {t, z, x, y},
      {t, z, y, x}, {x, t, y, z}, {x, t, z, y}, {x, y, t, z}, {x, y, z, t},
      {x, z, t, y}, {x, z, y, t}, {y, t, x, z}, {y, t, z, x}, {y, x, t, z},
      {y, x, z, t}, {y, z, t, x}, {y, z, x, t}, {z, t, x, y}, {z, t, y, x},
      {z, x, t, y}, {z, x, y, t}, {z, y, t, x}, {z, y, x, t}};

  (* Find all u[_] and Derivatives of u[_], even if they are nested in  *)
  (* a function. Find all other functions of the independent            *)
  (* variables.                                                         *)
  termsNestedInFunctions = Union[Flatten[Table[Map[Cases[expnFnList,
      Nest[_, #, i1] | _[Nest[_, #, i1], Nest[_, #, i1]] |
      _[Nest[_, #, i1], Nest[_, #, i1], Nest[_, #, i1]] |
      _[Nest[_, #, i1], Nest[_, #, i1], Nest[_, #,i1], Nest[_, #, i1]] /.
      List -> Sequence, {0, Infinity}] &, possibleIndVarCombos], {i1, 1, 4}]],
      Flatten[Map[Cases[expnFnList, Power[_, _[#]] | _[Power[_, _[#]]] /.
      List -> Sequence, {0, Infinity}] &, possibleIndVarCombos]]];
  printSTRIP["At STRIP 4, all terms containing u[_], derivatives of u[_],"<>
      " or functions of the independent variables:"];
  printSTRIP[termsNestedInFunctions];

  (* If rational expressions have dependent variables in the            *)
  (* denominator, zeros can occur with replacement                      *)
  Off[Power::infy];

  (* Create a list of coefficients - not including any weighted         *)
  (* parameters. Two rules are developed to remove all unwanted parts.  *)
  rules1 = Map[Rule[#, {}] &, Flatten[termsNestedInFunctions]];
  rules2 = Apply[List, Map[Rule[#, {}] &, paramswithweight]];
  rules3 = {x -> {}, y -> {}, z -> {}, t -> {}};
  coefficientList = expnFnList /. rules1;
  printSTRIP["At STRIP 5, dependent variables and the derivatives of"<>
      " dependent variables have been changed to {}:"];
  printSTRIP[coefficientList];

  coefficientList = coefficientList /. rules2;
  printSTRIP["At STRIP 6, any weighted parameters have been changed to {}:"];
  printSTRIP[coefficientList];

  coefficientList = coefficientList /. rules3;
  printSTRIP["At STRIP 7, explicit independent variables, x, y, z, or t,"<>
      " have been changed to {}:"];
  printSTRIP[coefficientList];

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

  (* Combine the coefficients for each term.                            *)
  singleNoCoefficientList = Apply[Times, coefficientList, {1, 1}];
  printSTRIP["At STRIP 9, a combined coefficient list is:"];
  printSTRIP[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];
  printSTRIP["At STRIP 10, a list of stripped terms is:"];
  printSTRIP[strippedTermsList];

  (* Clear all local variables not being returned.                      *)
  Clear[expandExpr, expnList, rules1, rules2, rules3, coefficientList,
      singleNoCoefficientList, possibleIndVarCombos, possibleDepVarCombos,
      permutDepVarCombos, possibleVarCombos,expnFnList, n1, printSTRIP];

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

(************************************************************************** *)
(* welcome: prints a welcome message                                        *)
(*                                                                          *)
(* Modules called upon - none                                               *)
(************************************************************************** *)

welcome := Block[{},
Print["********************************************************************"];
Print["WELCOME TO THE MATHEMATICA SOFTWARE PACKAGE HOMOGENIZE-AND-SOLVE.M  "];
Print["        by ANDREW COOK, WILLY HEREMAN, and UNAL GOKTAS.             "];
Print["   THE CODE COMPUTES SOLITARY WAVE SOLUTIONS AND SOLITON SOLUTIONS  "];
Print["                OF SOME NONLINEAR EVOLUTION EQUATIONS               "];
Print["          BASED ON A SIMPLIFIED VERSION OF HIROTA'S METHOD.         "];
Print["               Version 1 released on October 14, 2012               "];
Print["********************************************************************"];
]; (* end welcome *)

(****************************************************************************)

(* Code to run program STARTS HERE! *)

Clear[u];
Clear[\[Omega]];
Clear[\[Epsilon]];
Clear[F];
Clear[f];
Clear[g];
welcome;
menu;
Print["Read in the menu."];

(* AC: 2012/06/04 put pdeform back when adjusted to Mathematica v8.         *)
Print["Working with the data file for the ",name,"."]; 
For[i = 1, i <= noeqs, i++,
   eq[i][x,t] = Expand[eq[i][x,t]];
   Print["Equation ", i, " of the system with ", noeqs, " equation(s):"];
   Print["  ", 
         (* pdeform[eq[i][x,t]], *)
         eq[i][x,t], " = 0"]; (* Call to pdeform *)
   ];
eq[x,t] = Sum[eq[i][x,t],{i,1,noeqs}];

highestordereq = highestorder[eq[x,t]];
Print["Highest order of derivative in the system is ", highestordereq, "."];
Print["---------------------------------------------------"];

(* Call to 'homogenize' module.                                             *)
homogenizedList = homogenize[eq[x,t]];
homogenizedEqn = Part[homogenizedList,1];
homogenizedOrder = Part[homogenizedList,2];
kRule = Part[homogenizedList,3];
uSolnWithF[x_,t_] = Part[homogenizedList,4];
Print["---------------------------------------------------"];
Print["The selected homogenized PDE is: "];
Print[homogenizedEqn, " = 0"];
Print["---------------------------------------------------"];
checkSolns = Input["The code attempts to find the F[x,t] that satisfies " <>
                   "the homogenized PDE. Once found, would you " <>
                   "like the program to verify that F[x,t] and " <>
                   "u[x,t] are correct? Enter 'y' for yes and 'n' for no."];
If[checkSolns === y,
   debugCheckSolns = True;
  ];

moduleChoice = 
  Input["The code has 2 different ways of solving for soliton " <>
        "solutions. If you wish to use the faster approach " <>
        "with operators, please enter '1', otherwise enter '2' " <>
        "for the more reliable (but slower) approach."];

(* Call to 'solitonsoln' or 'solitonsolnOper' module.                       *)
If[moduleChoice === 1,
(* then1 *)
   fSoln[x_,t_] = solitonsolnOper[homogenizedEqn,homogenizedOrder],
(* else1 *)
   If[moduleChoice === 2,
   (* then2 *)
      fSoln[x_,t_] = solitonsoln[homogenizedEqn,homogenizedOrder],
   (* else2 *)
   Print["Invalid response, code is aborting!"];
   Abort[];
     ];
  ];

(* Checks F[x,t] and u[x,t] solutions found in the code.                    *)
If[debugCheckSolns,
   time = TimeUsed[];
   Print["---------------------------------------------------"];
   Print["Start of verification that the original PDE is satisfied."];
   Clear[F];
   Clear[u];
   Print["The original PDE was: "];
   Print[eq[x,t]];
   uSoln[x_,t_] = Factor[uSolnWithF[x,t] /. kRule];
   Print["Recall that u[1][x,t], with the selected K substituted, is: "];
   Print[uSoln[x,t]];
   Print["The code had computed F[x,t]: "];
   Print[fSoln[x,t]];
   eqSubU[x_, t_] = Factor[eq[x, t] /. {u[1][x, t] -> uSoln[x, t], 
       Derivative[m_, n_][u[1]][x, t] -> Derivative[m, n][uSoln][x, t]}];
   Print["Substituting u[1][x,t] into the original PDE gives: "];
   Print[eqSubU[x,t]];
   eqSubF[x_, t_] = Factor[Numerator[eqSubU[x, t]] /. {F[x, t] -> fSoln[x, t], 
       Derivative[m_, n_][F][x, t] -> Derivative[m, n][fSoln][x, t]} /.
       phiDerivRules];
   Print["Substituting F[x,t] into this result should result in 0 (zero). " <>
         "Substitution gives (this verification may take some time!): "];
   Print[eqSubF[x,t]];
   Print["Time used to check solutions: ", TimeUsed[]-time, " seconds."]; 
  ]; (* End If[debugCheckSolns] *)

(* End of code to run program *)

Print["*********************************************************************"];
Print["Code HOMOGENIZE-AND-SOLVE.M of October 14, 2012 "<>
      "was successfully loaded."];
Print["*********************************************************************"];

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