(* ************************************************************************ *)
(* : Title        : PDESolutionTester.m                                     *)
(* : Authors      : Jason Blevins, Jeff Heath, Matthew Porter-Peden         *)
(*                : Department of Mathematical and Computer Sciences        *)
(*                : Colorado School of Mines                                *)
(*                : Golden, Colorado 80401-1887, U.S.A                      *)
(* : Contact      : whereman@mines.edu, jrblevin@math.ncsu.edu              *)
(* : Last updated : 2002-06-24 15:45 by Blevins @ CTLM                      *)
(*                :                                                         *)
(* : Summary      : PDESolutionTester tests solutions to see if they        *)
(*                : satisfy the given system of PDEs. If there are          *)
(*                : unresolved parameters in the system or solution,        *)
(*                : PDESolutionTester attempts to solve them.               *)
(*                :                                                         *)
(* : Warnings     : Does not work well with Log solutions because of the    *)
(*                : lack of a nice series expansion.  If the package hangs  *)
(*                : while solving the coefficient system, consider reducing *)
(*                : the SolveOrder option to 3.                             *)
(* ************************************************************************ *)

BeginPackage["Calculus`PDESolutionTester`"]

Unprotect[PDESolutionTester]

PDESolutionTester::usage = 
"Usage:
\[NewLine]PDESolutionTester[system, depVars, indepVars, systemParams,
solutions, solutionParams, nonZeros, paramConstraints, simplificationRules,
options]
\[NewLine]
\[NewLine]Short Form (scans for variables):
\[NewLine]PDESolutionTester[system, solutions, nonZeros, options]
\[NewLine]
\[NewLine]PDESolutionTester tests the solutions to see if they satisfy the
given system of PDEs. If there are unresolved parameters in the system or
solution, PDESolutionTester attempts to solve them.
\[NewLine]PDESpecialSolutions takes the option Verbose with the default 
value False. If you wish to see output of the results, set Verbose true.
\[NewLine]PDESolutionTester takes the options NumericTest and 
SymbolicTest with the default values False. If NumericTest is set to True, 
PDESpecialSolutions tests the possible solutions based on 13 different sets of 
random numbers ranging from zero to one for all parameters and 
remaining constants. Solutions are accepted if they pass one or more of 
these tests.
\[NewLine]If SymbolicTest is set to True the solutions are tested truly
symbolically and the result of the test is shown in factored form. 
\[NewLine]In addition, PDESpecialSolutions has the option PrintSolutions with 
the default value False. If PrintSolutions is set to True, the correct
solutions will be printed in a more readable way.
\[NewLine]PDESolutionTester also takes the options SeriesOrder and SolveOrder
which control the order of the taylor series that is substituted and the
number of coefficients from the resulting system that are solved.  Default
values are SeriesOrder = 4 and SolveOrder = 4."

(* Options definitions. *)
Options[PDESolutionTester] = {
    Verbose -> False,
    DebugLevel -> 0,
    NumericTest -> False,
    SymbolicTest -> False,
    PrintSolutions -> False,
    SeriesOrder -> 4,
    SolveOrder -> 4
};

Begin["`Private`"]

(* ************************************************************************ *)
(* : Title   : PDESolutionTester                                            *)
(* : Author  : Jason Blevins                                                *)
(* : Input   : system - a system of PDEs as a List                          *)
(*           : depVars - a List containing the dependent variables          *)
(*           : indepVars - a List containing the independent variables      *)
(*           : systemParams - a List of parameters in the system            *)
(*           : solutions - a List of solutions                              *)
(*           : solutionParams - a List containing parameters in solutions   *)
(*           : nonZeros - a List of parameters that can be nonZero          *)
(*           : paramConstraints - a List of rules for parameter constraints *)
(*           : simplificationRules - a List of simplification rules         *)
(*           : opts - a List of options                                     *)
(* : Output  : A List of constraints to form correct solutions              *)
(* : Summary : PDESolutionTester has two function prototypes, the first,    *)
(*           : with 4 parameters runs the system and solution scanner to    *)
(*           : find the parameters and variables in the system, the second  *)
(*           : has more parameters, but is more flexible and reliable.      *)
(* ************************************************************************ *)

PDESolutionTester[system_List, solutions_List, nonZeros_List,
    opts___?OptionQ]:=
    Module[
        { debugLevel  = DebugLevel /. {opts} /. Options[PDESolutionTester],
          verboseFlag = Verbose /. {opts} /. Options[PDESolutionTester],
          depVars,
          indepVars,
          systemParams,
          solutionParams,
          result
        },

        (* ******************************************************** *)
        (* Scan the system for independent, dependent & parameters  *)
        (* ******************************************************** *)
        If[debugLevel >= 3,
            Print["**DEBUG3** (ST0, PDESolutionTester): Sending system to",
                  " the system scanner..."];
        ];

        If[verboseFlag,
            Print["Scanning system and solution for variables..."];
        ];

        {depVars, indepVars, systemParams} =
            PDESolutionTesterSystemScanner[system, opts];

        If[debugLevel >= 3,
            Print["**DEBUG3** (ST0.1, PDESolutionTester): System scanner",
                  " results:  depVars = ", depVars, ", indepVars = ",
                  indepVars, ", systemParams = ", systemParams];
        ];

        (* ******************************************************** *)
        (* Scan solution for independent, dependent & parameters    *)
        (* ******************************************************** *)
        If[debugLevel >= 3,
            Print["**DEBUG3** (ST0.2, PDESolutionTester): Sending solution",
                  " to the solution scanner..."];
        ];

        solutionParams = PDESolutionTesterSolutionScanner[solutions,
            depVars, indepVars, opts];

        If[debugLevel >= 3,
            Print["**DEBUG3** (ST0.3, PDESolutionTester): Solution Scanner",
                  " results: solutionParams = ", solutionParams];
        ];

        (* ******************************************************** *)
        (* Call the main SolutionTester method results from scanner *)
        (* ******************************************************** *)
        result =
            PDESolutionTester[system, depVars, indepVars,
                systemParams, solutions, solutionParams, nonZeros,
                {}, {}, opts];
        
        Return[result];

    ]; (* End Module *)

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

PDESolutionTester[system_List, depVars_List, indepVars_List,
    systemParams_List, solutions_List, solutionParams_List,
    nonZeros_List, paramConstraints_List, simplificationRules_List,
    opts___?OptionQ]:=
    Module[
        { debugLevel  = DebugLevel /. {opts} /. Options[PDESolutionTester],
          verboseFlag = Verbose /. {opts} /. Options[PDESolutionTester],
          printSolutionsFlag = PrintSolutions /. {opts} /.
              Options[PDESolutionTester],
          numericFlag = NumericTest /. {opts} /. Options[PDESolutionTester],
          symbolicFlag = SymbolicTest /. {opts} /. Options[PDESolutionTester],
          seriesOrder = SeriesOrder /. {opts} /. Options[PDESolutionTester],
          solveOrder = SolveOrder /. {opts} /. Options[PDESolutionTester],
          newSystem,
          seriesSystem,
          potentialSolutions,
          validConstraints,
          numConstraints
        },

        (* ******************************************************** *)
        (* Sanity Checks                                            *)
        (* ******************************************************** *)
        If[seriesOrder < solveOrder,
            CellPrint[
                Cell["Warning: SeriesOrder is less than SolveOrder.",
                     "Message"]
            ];
        ];

        (* ******************************************************** *)
        (* If verbose, prints system in readable form.              *)
        (* ******************************************************** *)
        If[verboseFlag,
            PDESolutionTesterPrintSystem[system,
                "The given system of PDEs is: ",
                solutions,
                "The given solution set is: ",
                depVars, opts
            ];
        ];

        (* ******************************************************** *)
        (* Convert system into homogeneous expressions              *)
        (* ******************************************************** *)
        newSystem = system /. (a__==b__):>(a-b);

        If[debugLevel >=3,
            Print["**DEBUG3** (ST2, PDESolutionTester): newSystem = ",
            newSystem];
        ];

        (* ******************************************************** *)
        (* Step 1:                                                  *)
        (* Substitute solutions, series, & constraints into system  *)
        (* ******************************************************** *)
        If[verboseFlag,
            Print["Step 1: Applying Taylor Approximations to the solutions",
                  " and substituting them into the system..."];
        ];

        seriesSystem = PDESolutionTesterApplySolutions[newSystem, indepVars,
            solutions, paramConstraints, opts];

        (* ******************************************************** *)
        (* Step 2:                                                  *)
        (* Solve for any obstructions                               *)
        (* ******************************************************** *)
        If[verboseFlag,
            Print["Step 2: Solving obstructions..."];
        ];

        potentialSolutions = PDESolutionTesterFindConstraints[seriesSystem,
            indepVars, systemParams, solutionParams, nonZeros, opts];

        If[verboseFlag,
            Print["Potential constraints are: ", potentialSolutions];
        ];

        (* ******************************************************** *)
        (* Step 3:                                                  *)
        (* Numeric and symbolic tests on potential solutions        *)
        (* ******************************************************** *)
        If[verboseFlag,
            Print["Step 3: Test potential solutions..."];
        ];

        validConstraints = PDESolutionTesterTestSolutions[newSystem,
            solutions, paramConstraints, potentialSolutions, indepVars,
            Union[systemParams, solutionParams], simplificationRules, opts];

        (* ******************************************************** *)
        (* Analyze and print results                                *)
        (* ******************************************************** *)
        If[verboseFlag,
            If[Length[validConstraints] === 0,
                Print["Solution is incorrect and could not be corrected."];
            ];

            If[validConstraints === {{}},
                Print["Solution is correct."];

                (* Else *),
                If[(numConstraints = Length[validConstraints]) > 0,
                    If[(!numericFlag && !symbolicFlag),
                        Print["There were ", numConstraints,
                              " UNTESTED potential solutions found: ",
                              validConstraints];

                        (* Else *),
                        Print["There were ", numConstraints,
                              " valid constraints found: ",
                              validConstraints];
                    ];

                    If[printSolutionsFlag,
                        Do[
                            Print["******************************"];
                            Print["Applying constraint ", i, " yields:"];

                            PDESolutionTesterPrintSystem[
                                (system /. paramConstraints) /.
                                validConstraints[[i]],
                                "System:",
                                (solutions /. paramConstraints) /.
                                validConstraints[[i]],
                                "Solution:",
                                depVars, opts
                            ];
                            (* Iterator *),
                            {i, 1, numConstraints}
                        ];
                    ];
                ];
            ];
        ];

        Return[validConstraints];

    ]; (* End Module *)

(* End of function PDESolutionTester. *)

(* ************************************************************************ *)
(* : Title   : PDESolutionTesterApplySolutions                              *)
(* : Author  : Jason Blevins and Jeff Heath                                 *)
(* : Input   : system, indepVars, solutions, paramConstraints, opts         *)
(* : Output  : A polynomial system                                          *)
(* : Summary : ApplySolutions applies any constraints to the system and     *)
(*           : solution, then it replaces elementary functions in the       *)
(*           : solution with their taylor series approximatins and then     *)
(*           : substitutes the new solution into the system.                *)
(* ************************************************************************ *)

PDESolutionTesterApplySolutions[system_, indepVars_, solutions_,
    paramConstraints_, opts___?OptionQ] :=
    Module[
        { debugLevel  = DebugLevel /. {opts} /. Options[PDESolutionTester],
          verboseFlag = Verbose /. {opts} /. Options[PDESolutionTester],
          seriesOrder = SeriesOrder /. {opts} /. Options[PDESolutionTester],
          highestDeriv,
          trigRules,
          jacobiRules,
          seriesRules,
          newSystem,
          newSolutions,
          toPure
        },

        (* ******************************************************** *)
        (* Calculate highest derivative of the system and create    *)
        (* taylor series rules of seriesOrder orders higher         *)
        (* (this sums mixed derivatives)                            *)
        (* ******************************************************** *)
        highestDeriv = Max[Apply[Plus, #] & /@
                (Cases[Variables[system], Derivative[__][_][__]] /.
                Derivative[m___][_][__] -> {m})];
        seriesOrder = highestDeriv + seriesOrder;

        If[debugLevel >= 2,
            Print["**DEBUG2** (AS1, ApplySolutions): Highest Derivative in",
                  " the system, highestDerivative = ", highestDeriv,
                  ", taking Taylor series of order, seriesOrder = ",
                  seriesOrder];
        ];

        trigRules = #[z_] -> Normal[Series[#[z], {z, 0, seriesOrder}]]& /@
            {Exp, Sin, Cos, Tan, Cot, Sec, Csc, Sinh, Cosh, Tanh, Sech,
            Csch, Coth};
        jacobiRules = #[z_, m_] -> Normal[Series[#[z,m],
            {z, 0, seriesOrder}]]& /@ {JacobiSN, JacobiCN, JacobiDN};

        seriesRules = Union[trigRules, jacobiRules];

        If[debugLevel >= 3,
            Print["**DEBUG3** (AS2, ApplySolutions): Taylor Series Rules,",
                  " seriesRules = ", PDESolutionTesterCleanUp[seriesRules]];
        ];

        (* ******************************************************** *)
        (* Apply constraints to system and solutions                *)
        (* ******************************************************** *)
        newSystem = system /. paramConstraints;

        If[debugLevel >=3,
            Print["**DEBUG3** (AS3, ApplySolutions): After substitution of",
                  " constraints, newSystem = ", newSystem];
        ];

        newSolutions = solutions /. paramConstraints;

        If[debugLevel >=3,
            Print["**DEBUG3** (AS4, ApplySolutions): After substitution of",
                  " constraints, newSolutions = ", newSolutions];
        ];

        (* ******************************************************** *)
        (* Apply series rules to solutions                          *)
        (* ******************************************************** *)
        newSolutions = MapAll[Expand, (newSolutions /. seriesRules)];

        If[debugLevel >= 3,
            Print["**DEBUG3** (AS5, ApplySolutions): After applying",
                  " seriesRules and expanding",
                  ", newSolutions = ", newSolutions];
        ];

        (* ******************************************************** *)
        (* Convert solutions to pure functions & substitute in sys  *)
        (* toPure function by Douglas Baldwin                       *)
        (* ******************************************************** *)
        toPure = (# /. (a__[var__]->temp__) :> (a :> Function[{var}, temp]))&;

        newSystem = newSystem /. toPure[newSolutions];

        If[debugLevel >=2,
            Print["**DEBUG2** (AS6, ApplySolutions): After substitution",
                  " of solutions into system, newSystem = ", newSystem];
        ];

        Return[newSystem];

    ]; (* End Module *)

(* End of function PDESolutionTesterApplySolutions *)

(* ************************************************************************ *)
(* : Title   : PDESolutionTesterFindConstraints                             *)
(* : Author  : Jason Blevins and Jeff Heath                                 *)
(* : Input   : system - A List of PDE expressions,                          *)
(*           : indepVars - A List containing the independent variables      *)
(*           : systemParams - A List containing the system's parameters     *)
(*           : solutionParams - A List containing the solution's parameters *)
(*           : nonZeros - A List containing parameters that are nonzero     *)
(*           : opts - A List of options passed to the package               *)
(* : Output  : A List of Lists of Rules, potential constraints on params    *)
(* : Summary : FindConstraints takes a system of polynomials in indepVars   *)
(*           : sets indepVars to zero finds the first 5 coefficients and    *)
(*           : solves any obstructions it finds                             *)
(* ************************************************************************ *)

PDESolutionTesterFindConstraints[system_, indepVars_, systemParams_,
    solutionParams_, nonZeros_, opts___?OptionQ] :=
    Module[
        { debugLevel  = DebugLevel /. {opts} /. Options[PDESolutionTester],
          verboseFlag = Verbose /. {opts} /. Options[PDESolutionTester],
          solveOrder = SolveOrder /. {opts} /. Options[PDESolutionTester],
          zeroRules,
          newSystem,
          coeffSys,
          time,
          result
        },

        (* ******************************************************** *)
        (* Create rules to individually consider polynomials of     *)
        (* each independent variable                                *)
        (* ******************************************************** *)
        zeroRules =
            Table[# -> 0 & /@
                Complement[indepVars, {Part[indepVars, i]}]

                (* Iterator *),
                {i, 1, Length[indepVars]}
            ];

        If[debugLevel >= 3,
            Print["**DEBUG3** (FC1, FindConstraints): zeroRules = ",
            zeroRules];
        ];

        (* ******************************************************** *)
        (* Apply zeroRules to the system, returns a list of systems *)
        (* of polynomials in each variable of form:                 *)
        (* {{List of eqs in x1}, {List of eqs in x2}, etc... }      *)
        (* ******************************************************** *)
        newSystem = system /. zeroRules;
        newSystem = MapAll[Numerator, MapAll[Factor,
            MapAll[Expand, newSystem]]];

        If[debugLevel >= 2,
            Print["**DEBUG2** (FC2, FindConstraints): After applying",
            " zeroRules, newSystem = ", newSystem];
        ];

        (* ******************************************************** *)
        (* Creating a system from the coefficients on the           *)
        (* independent variables in each equation of form:          *)
        (* { { {coeff list eq1 x1}, {coeff list eq2 x1}, ... },     *)
        (*   { {coeff list eq1 x2}, ... }, ... }                    *)
        (* ******************************************************** *)
        coeffSys =
            Table[
                CoefficientList[#, indepVars[[i]]]& /@
                newSystem[[i]],
                {i, 1, Length[indepVars]}
            ];

        If[debugLevel >= 3,
            Print["**DEBUG3** (FC2.5, FindConstraints): coeffSys = ",
            coeffSys];
        ];

        (* ******************************************************** *)
        (* Because CoefficientList[0] returns {} and                *)
        (* CoefficientList[constant] returns {constant}, we must    *)
        (* Pad each coefficient list to make it solveOrder long     *)
        (* ******************************************************** *)
        coeffSys = Map[PadRight[#, solveOrder]&, coeffSys, {2}];

        If[debugLevel >= 3,
            Print["**DEBUG3** (FC2.75, FindConstraints): coeffSys = ",
            coeffSys];
        ];

        (* ******************************************************** *)
        (* Only take the first solveOrder terms from each eqn       *)
        (* ******************************************************** *)
        coeffSys = Map[Take[#, solveOrder] &, coeffSys, {2}];

        If[debugLevel >= 2,
            Print["**DEBUG2** (FC3, FindConstraints): coeffSys = ", coeffSys];
        ];

        (* ******************************************************** *)
        (* Factor the list, remove any constant coefficient then    *)
        (* Union to remove duplicates                               *)
        (* ******************************************************** *)
        coeffSys = MapAll[Factor, Union[Flatten[Map[Drop[#, 1]&,
                       Map[FactorTermsList, coeffSys, {3}], {3}]]]];

        If[debugLevel >= 1,
            Print["**DEBUG1** (FC3.5, FindConstraints): coeffSys = ",
                  coeffSys];
        ];

        (* ******************************************************** *)
        (* Solve the obstructions using AnalyzeAndSolve             *)
        (* ******************************************************** *)
        time = TimeUsed[];
        result = Algebra`AnalyzeAndSolve`AnalyzeAndSolve[coeffSys,
                solutionParams, systemParams,{}, nonZeros];
        result = MapAll[Factor, result];
        time = TimeUsed[] - time;

        If[debugLevel >= 1,
            Print["**DEBUG1** (FC4, FindConstraints): AnalyzeAndSolve",
                  " took ", time, " seconds to find ", Length[result],
                  " potential constraints, result = ", result];
        ];

        Return[result];

    ]; (* End Module *)

(* End of function PDESolutionTesterFindConstraints. *)

(* ************************************************************************ *)
(* : Title   : PDESolutionTesterTestSolutions                               *)
(* : Author  : Jason Blevins                                                *)
(* : Input   : system, solutions, paramConstraints, potentialSolutions,     *)
(*           : indepVars, params, userSimpRules, opts                       *)
(* : Output  : returns elements from potentialSolutions that passed tests   *)
(* : Summary : TestSolutions performs a series of numeric and symbolic      *)
(*           : tests to determine which of the potentialSolutions are valid *)
(* ************************************************************************ *)

PDESolutionTesterTestSolutions[system_, solutions_, paramConstraints_,
    potentialSolutions_, indepVars_, params_, userSimpRules_,
    opts___?OptionQ] :=
    Module[
        { debugLevel  = DebugLevel /. {opts} /. Options[PDESolutionTester],
          verboseFlag = Verbose /. {opts} /. Options[PDESolutionTester],
          numericFlag = NumericTest /. {opts} /. Options[PDESolutionTester],
          symbolicFlag = SymbolicTest /. {opts} /. Options[PDESolutionTester],
          simplificationRules,
          solns = {},
          warn = {},
          solutionRules,
          newSystem,
          toPure,
          newSolutions,
          randomRules,
          numericResults,
          symbolicResults,
          time
        },

        (* ******************************************************** *)
        (* Apply constraints to system and solutions                *)
        (* ******************************************************** *)
        newSystem = system /. paramConstraints;

        newSolutions = solutions /. paramConstraints;

        (* ******************************************************** *)
        (* Convert solutions to pure functions to allow them to be  *)
        (* differentiated when applied to the system                *)
        (*                                                          *)
        (* toPure was written by Douglas Baldwin                    *)
        (* ******************************************************** *)
        toPure = (# /. (a__[var__]->temp__) :> (a :> Function[{var}, temp]))&;

        newSystem = newSystem /. toPure[newSolutions];


        (* ******************************************************** *)
        (* Combine built in simplification rules with user          *)
        (* specified rules to reducing the solution even            *)
        (* further so as to find more general solutions.            *)
        (* Then apply those simplification rules to the solutions   *)
        (*                                                          *)
        (* Written by Douglas Baldwin for PDESpecialSolutions       *)
        (* Modified by Jason Blevins to simplify in JacobiCN only   *)
        (* ******************************************************** *)
        simplificationRules = {
              Sqrt[xxx__^2] :> 
                  (warn = Append[warn, "Sqrt[a^2]->a"]; xxx), 
              Sqrt[Power[yyy__,2]] :> 
                  (warn = Append[warn, "Sqrt[a^2]->a"]; yyy), 
              Sqrt[-zzz__^2] :> 
                  (warn = Append[warn, "Sqrt[-a^2]->I*a"]; I*zzz),
              Sqrt[-ttt_*zzz__^2] :> 
                  (warn = Append[warn, "Sqrt[-a*b^2]->I*b*Sqrt[a]"];
                  I*zzz*Sqrt[ttt]), 
              (1-Sech[xx__]^2) -> Tanh[xx]^2,
              (1-Tanh[xx__]^2) -> Sech[xx]^2,
              JacobiDN[xx__, mm__]^2 :> 1-mm+mm*JacobiCN[xx,mm]^2,
              JacobiSN[xx__, mm__]^2 :> 1-JacobiCN[xx,mm]^2,
              JacobiDN[xx__, mm__] :> 
                  (warn = Append[warn,
                  "JacobiDN[x,mod]->Sqrt[1-mod+mod*JacobiCN[x,mod]^2]"];
                  Sqrt[1-mm+mm*JacobiCN[xx,mm]^2]),
              JacobiSN[xx__, mm__] :> 
                  (warn = Append[warn,
                  "JacobiSN[x,mod]->Sqrt[1-JacobiCN[x,mod]^2]"];
                  Sqrt[1-JacobiCN[xx,mm]^2])
        };

        simplificationRules = Union[simplificationRules, userSimpRules];

        solutionRules =
            FixedPoint[
                MapAll[Expand,
                    MapAll[Factor,
                        # //. simplificationRules
                    ] //. simplificationRules
                ]&,
                potentialSolutions,
                3
            ];

        If[warn != {},
            CellPrint[
                Cell["Warning: The following simplification rules were used: " <>
                     ToString[Union[warn]], "Message"]
            ];
        ];

        If[debugLevel >= 3,
            Print["**DEBUG3** (TS1, TestSolutions): After applying" <>
                  " simplification rules to the potential solutions," <>
                  " solutionRules = ", solutionRules];
        ];

        (* ******************************************************** *)
        (* If no tests are specified, return potential solutions    *)
        (* ******************************************************** *)
        If[(!numericFlag && !symbolicFlag),

            If[debugLevel >= 2,
                Print["**DEBUG2** (TS2.5, TestSolutions): Skipping numeric" <>
                      " and symbolic tests, function will now return" <>
                      " solutionRules..."];
            ];

            CellPrint[
                Cell["These solutions are not being tested numerically or" <>
                     " symbolically. To test the solutions set either the" <>
                     " NumericTest option to True, or set the SymbolicTest" <>
                     " option to True, or both.",
                     "Message"]
            ];

            Return[Union[MapAll[Factor, solutionRules]]];
        ];

        (* ******************************************************** *)
        (* Perform numeric testing                                  *)
        (* ******************************************************** *)
        If[numericFlag,
            If[verboseFlag,
                Print["Numerically testing solutions..."];
            ];

            (* ************************************************ *)
            (* Apply solutionRules to the system, tag the rule  *)
            (* on the end of each copy of the system in which   *)
            (* it was substituted                               *)
            (* ************************************************ *)
            numericResults = {newSystem /. #, #}& /@ solutionRules;

            If[debugLevel >= 3,
                Print["**DEBUG3** (TS3, TestSolutions):",
                      " numericResults = ", numericResults];
            ];

            (* ************************************************ *)
            (* Create 13 lists of random numbers to test with   *)
            (* ************************************************ *)
            randomRules = Table[
                (#->Random[])& /@
                Union[indepVars, params],
                {13}
            ];

            If[debugLevel >= 3,
                Print["**DEBUG3** (TS4, TestSolutions): randomRules = ",
                    randomRules];
            ];

            (* ************************************************ *)
            (* Apply these random numbers to the system,        *)
            (* while keeping the tagged solutionRule at the end *)
            (* ************************************************ *)
            numericResults = {#[[1]] /. randomRules, #[[2]]} & /@
                numericResults;

            If[debugLevel >= 2,
                Print["**DEBUG2** (TS5, TestSolutions): After applying",
                      " randomRules, numericResults = ",
                      numericResults];
            ];

            (* ************************************************ *)
            (* If any one of the results is less than 10^-10,   *)
            (* we will keep the solution                        *)
            (* ************************************************ *)
            If[Min[MapAll[Abs, #[[1]]]] < 10^-10,
                solns = Append[solns, #[[2]]];
            ] & /@ numericResults;

            If[debugLevel >= 1,
                Print["**DEBUG1** (TS6, TestSolutions): Correct solutions",
                      " after numeric testing are, solns = ", solns];
            ];
        ]; (* End numericFlag If *)

        (* ******************************************************** *)
        (* Perform symbolic testing                                 *)
        (*                                                          *)
        (* Based on code by Douglas Baldwin in PDESpecialSolutions  *)
        (* ******************************************************** *)
        If[symbolicFlag,
            If[verboseFlag,
                Print["Symbolically testing solutions..."];
            ];

            (* ************************************************ *)
            (* Create a list of the solutions applied to the    *)
            (* system, with the solution tagged to the end      *)
            (* ************************************************ *)
            symbolicResults = {newSystem /. #, #}& /@ solutionRules;

            If[debugLevel >= 3,
                Print["**DEBUG3** (TS6.5, TestSolutions): After building" <>
                      " symbolic test cases, symbolicResults = ",
                      symbolicResults];
            ];

            (* ************************************************ *)
            (* Symbolically test the solutions & simplify then  *)
            (* place those that simplify to zero in solns       *)
            (* ************************************************ *)
            time = TimeUsed[];

            symbolicResults =
                ExpToTrig[
                    FixedPoint[
                        MapAll[Factor,
                            MapAll[Expand,TrigToExp[#]] /. simplificationRules
                        ] /. simplificationRules &, #]&
                    /@ symbolicResults
                ];

            time = TimeUsed[] - time;

            If[debugLevel >= 2,
                Print["**DEBUG2** (TS7, TestSolutions): After Symbolic",
                    " testing, symbolicResults = ",
                    symbolicResults];
            ];

            solns =
                Join[solns,
                    Cases[symbolicResults,
                        {{(0)..}, _List},
                        Infinity
                    ] //. {{(0)..}, a_List} :> a
                ];

            If[debugLevel >= 1,
                Print["**DEBUG1** (TS8, TestSolutions): After Symbolic",
                    " testing (", time, " s), solns = ", solns];
            ];

            (* ************************************************ *)
            (* Removes all the good cases, so we may output the *)
            (* bad ones, if any, to the user.                   *)
            (* ************************************************ *)
            symbolicResults =
                DeleteCases[
                    DeleteCases[symbolicResults,
                        {{(0)..}, _List},
                        Infinity
                    ],
                    {}
                ];

            If[debugLevel >= 2,
                Print["**DEBUG2** (TS9, TestSolutions): Solutions which" <>
                      " did not simplify to zero, symbolicResults = ",
                      symbolicResults];
            ];

            If[symbolicResults != {},
                CellPrint[
                    Cell["These solutions did not simplify to zero." <>
                        " Please test these equations by hand or" <>
                        " interactively with Mathematica.",
                        "Message"
                    ]
                ];
                Print[symbolicResults];
            ];
        ]; (* End symbolicFlag If *)

    Return[Union[MapAll[Factor, solns]]];

    ]; (* End Module *)

(* End of function PDESolutionTesterTestSolutions. *)

(* ************************************************************************ *)
(* : Title   : PDESolutionTesterSystemScanner                               *)
(* : Author  : Matthew Porter-Peden                                         *)
(* : Input   : system (a list of equations)                                 *)
(* : Output  : Returns {depVars,indepVars, systemParams}                    *)
(* : Summary : finds the dependent variables and system params from system  *)
(* ************************************************************************ *)

PDESolutionTesterSystemScanner[system_, opts___?OptionQ]:=
    Module[
        { debugLevel  = DebugLevel /. {opts} /. Options[PDESolutionTester],
          verboseFlag = Verbose /. {opts} /. Options[PDESolutionTester],
          vari1,
          depVar,
          indepVar,
          depVars,
          indepVars,
          systemParams,
          newsystem,
          res,
          result,
          depPos
        },

        If[debugLevel >= 3, Print["At Pt. Scan1 system is: ", system]];

        (* Converts System into expressions that this function can use. *)
        newsystem = system /. (a__ == b__) :> (a - b);

        If[debugLevel >= 3, Print["At Pt. Scan2 newsystem is: ", newsystem]];

        (*Turns Powers into multiplication so the Variables function *)
        (*could find terms like the x in e^x*)
        newsystem = newsystem /. Power -> Times;

        If[debugLevel >= 3,
            Print["At Pt. Scan3 modified newsystem is: ", newsystem]
        ];

        (*Picks out all the symbols in the new system*)
        vari1 = Variables[newsystem];

        If[debugLevel >= 3,
            Print["At Pt. Scan4 vari0=Variables[newsystem]: ", vari1]
        ];

        (*Picks out the System's Parameters by removing everthing that is *)
        (*a function of something else*)
        systemParams = DeleteCases[vari1, u_[a___]];

        If[debugLevel >= 3,
            Print["At Pt. Scan5 systemParams: ", systemParams]
        ];

        (*Picks out the everything that is a function of something else*)
        depVar = Cases[vari1, u_[a___]];

        If[debugLevel >= 3,
            Print["At Pt. Scan6 depVars: ", depVar]
        ];

        (*Converts derivatives, i.e. u' to u*)
        depVar = depVar /. Derivative[a___][u_][b___] -> u[b];

        If[debugLevel >= 3, Print["At Pt. Scan7 depVars: ", depVar]];

        (*Removes repeted variables*)
        depVar = Union[depVar];

        If[debugLevel >= 3, Print["At Pt. Scan8 depVars: ", depVar]];

        (*This Function removes the heads of functions*)
        (*Leaving the independant variables behind*)
        removeHead[u_[a___]] := a;

        (*Removes The Head on Each part of the depVar list*)
        (*It adds the contents together to get rid of the *)
        (*sequence and make it a single expression*)
        Do[
            res[i] = Plus[removeHead[depVar[[i]]]];

            If[debugLevel >= 3,
                Print["At Pt. Scan9.",i," res[", i, "] is: ",
                      res[i]];
            ];

            (* Iterator *),
            {i, 1, Length[depVar]}
        ];

        (*Checks to see if any of the contents are numbers*)
        Do[
            result[i] = NumericQ[res[i]]; 

            If[debugLevel >= 3,
                Print["At Pt. Scan10.",i," result[", i,"] is Number? ",
                      result[i]];
            ];

            (* Iterator *),
            {i, 1, Length[depVar]}
        ];

        (*If all the contents of an item are numbers, it is a perameter such *)
        (*as k[1] or a[1,2]*)
        (*It is then added to the systemParams list*)
        Do[
            If[result[i],
                systemParams = Append[systemParams, depVar[[i]]];
            ];

            (* Iterator *),
            {i, 1, Length[depVar]}
        ];

        (*Removes the system Parameters from the depVars list*)
        depVar = Complement[depVar, systemParams];

        If[debugLevel >= 3, Print["At Pt. Scan11 depVars: ", depVar]];

        (*The next three steps use the removeHead function to get at the *)
        (*independent variables*)
        indepVar = Union[MapAll[removeHead, depVar]];
        indepVar = List[removeHead[indepVar]];

        If[debugLevel >= 3,Print["At Pt. Scan12 indepVars: ",indepVar]];

        (*Removes any numbers in the independent vars. list*)
        indepVars=Variables[indepVar];

        If[debugLevel >= 3,Print["At Pt. Scan13 indepVar: ",indepVar]];

        (*Substitutes random numbers in for the independent vars.*)
        (*So it can find false Dependent vars. such as Cos[x] or *)
        (*JacobiCN[x,mod]*)
        depPos = Position[depVar /. ((# -> Random[]) & /@ Join[indepVar]),
            a_?NumericQ, 1];

        If[debugLevel >= 3,
            Print["At Pt. Scan14 position of false depVar: ", depPos];
            If[Length[depPos] > 0, Print["Deleting false depVar"]];
        ];

        (*And then it deletes any false Dependent vars.*)
        depVars = Delete[depVar, depPos];

        (*Makes certain that no independent variables got picked up as *)
        (*Parameters*)
        systemParams = Union[Complement[systemParams, indepVars]];


        (*Prints all the results of this part*)
        If[verboseFlag,
            Print["Dependent Variables = ", depVars];
            Print["Independent Variables = ", indepVars];
            Print["System Parameters = ", systemParams];
        ];

        Return[{depVars, indepVars, systemParams}];

    ]; (* End Module *)

(* End of function PDESolutionTesterSystemScanner. *)

(* ************************************************************************ *)
(* : Title   : PDESolutionTesterSolutionScanner                             *)
(* : Author  : Matthew Porter-Peden                                         *)
(* : Input   : solutions, depVars, indepVars, opts                          *)
(* : Output  : List of parameters in solutions                              *)
(* : Summary : Run this after SystemScanner on solutions to determine any   *)
(*           : parameters in the solutions.                                 *)
(* ************************************************************************ *)

PDESolutionTesterSolutionScanner[solutions_, depVars_, indepVars_, 
        opts___?OptionQ]:=
    Module[
        { debugLevel  = DebugLevel /. {opts} /. Options[PDESolutionTester],
          verboseFlag = Verbose /. {opts} /. Options[PDESolutionTester],
          vari2,
          intermediateRes,
          sol,
          solres,
          solresult,
          solParams,
          brackets,
          insideBrackets,
          functionToListRules,
          insideParams,
          solutionParams
        }, (* Protected Local Variables *)

        If[debugLevel >= 3,
            Print["At Pt. SolScan1 solution is: ", solutions];
        ];

        (*Converts the Solution Rules to expressions*)
        sol = depVars /. solutions;

        If[debugLevel >= 3, Print["At Pt. SolScan2 sol is: ", sol]];

        (*Turns Powers into multiplication so the Variables function *)
        (*could find terms like the x in e^x*)
        sol = sol /. Power -> Times;

        If[debugLevel >= 3, Print["At Pt. Scan3 modified sol is: ", sol]];

        (*Picks out all the symbols*)
        vari2 = Variables[sol];

        If[debugLevel >= 3, Print["At Pt. SolScan3 vari2 is: ", vari2]];

        (*Extracts the Parameters that don't have any brackets*)
        solParams = DeleteCases[vari2, u_[a___]];

        If[debugLevel >= 3,
            Print["At Pt. SolScan4 solParams is: ", solParams];
        ];

        (*Picks out everything with brackets*)
        brackets = Cases[vari2, u_[a___]];

        If[debugLevel >= 3,
            Print["At Pt. SolScan5 brackets are: ", brackets];
        ];
            
        (*Removes the Head on each part of the depVar list*)
        (*It adds the contents together to get rid of the sequence and make *)
        (*it a single expression*)
        Do[
            solres[i] = Plus[removeHead[brackets[[i]]]]; 

            If[debugLevel >= 3,
                Print["At Pt. SolScan6.",i," solres[", i, "] is: ",
                      solres[i]];
            ];

            (* Iterator *),
            {i, 1, Length[brackets]}
        ];

        (*Checks to see if any of the contents are numbers*)
        Do[
            solresult[i] = NumericQ[solres[i]];

            If[debugLevel >= 3,
                Print["At Pt. SolScan7.", i, " solresult[", i,
                      "] is Number? ", solresult[i]];
            ];

            (* Iterator *),
            {i, 1, Length[brackets]}
        ];

        (*If all the contents of an item are numbers, it is a perameter such *)
        (*as c[1] or a[1,2]*)
        (*They are then added to the systemParams list*)
        Do[
            If[solresult[i],
                solParams = Append[solParams, brackets[[i]]];
            ]

            (* Iterator *),
            {i, 1, Length[brackets]}
        ];

        brackets = Complement[brackets, solParams];

        If[debugLevel >= 3,
            Print["At Pt. SolScan8 solParams is: ",solParams];
        ];

        (*A list of rules to convert functions into lists.  i.e. *)
        (*Tanh[x+y]->{x,y}*)
        functionToListRules = {ArcCos -> List, ArcCosh -> List,
            ArcCot -> List, ArcCoth -> List, ArcCsc -> List, ArcCsch -> List,
            ArcSec -> List, ArcSech -> List, ArcSin -> List, ArcSinh -> List,
            ArcTan -> List, ArcTanh -> List, Cos -> List, Cosh -> List,
            Coth -> List, Csc -> List, Csch -> List, Divide -> List,
            InverseJacobiCD -> List, InverseJacobiCN -> List,
            InverseJacobiCS -> List, InverseJacobiDC -> List,
            InverseJacobiDN -> List, InverseJacobiDS -> List,
            InverseJacobiNC -> List, InverseJacobiND -> List,
            InverseJacobiNS -> List, InverseJacobiSC -> List,
            InverseJacobiSD -> List, InverseJacobiSN -> List,
            JacobiCD -> List, JacobiCN -> List, JacobiCS -> List,
            JacobiDC -> List, JacobiDN -> List, JacobiDS -> List,
            JacobiNC -> List, JacobiND -> List, JacobiNS -> List,
            JacobiSC -> List, JacobiSD -> List, JacobiSN -> List,
            Log -> List, Minus -> List, NonCommutativeMultiply -> List,
            Plus -> List, Power -> List, Sec -> List, Sech -> List,
            Sin -> List, Sinh -> List, Sqrt -> List, Subtract -> List,
            Tan -> List, Tanh -> List, Times -> List
        };

        (*Applies functionToListRules to the brackets list to get the *)
        (*contents of each part of the brackets list*)
        insideBrackets = Union[Flatten[brackets /. functionToListRules]];

        If[debugLevel >= 3,
            Print["At Pt. SolScan9 insideBrackets is: ",insideBrackets];
        ];

        (*Removes any number from the inside brackets list*)
        insideBrackets = Variables[insideBrackets];

        If[debugLevel >= 3,
            Print["At Pt. SolScan10 insideBrackets is: ",insideBrackets];
        ];

        (*Removes the independent vars. from  insideBrackets and calls it *)
        (*insideParams*)
        insideParams = Complement[insideBrackets, indepVars];

        If[debugLevel >= 3,
            Print["At Pt. SolScan11 insideParams is: ", insideParams];
        ];

        (*Adds the inside Parameters to the solution Parameters list*)
        solParams = Flatten[Join[solParams, insideParams]];
        solParams = Union[solParams];

        (*Removes independent variables from solution peramerters list*)
        solutionParams = Complement[solParams, indepVars];

        If[debugLevel >= 3,
            Print["At Pt. SolScan12 solParams is: ", solutionParams];
        ];

        If[verboseFlag, Print["Solution Parameters = ", solutionParams]];

        Return[solutionParams];

    ]; (* End Module *)

(* End of function PDESolutionTesterSolutionScanner. *)

(* ************************************************************************ *)
(* : Title   : PDESolutionTesterPrintSystem                                 *)
(* : Author  : Douglas Baldwin, Jason Blevins                               *)
(* : Input   : system, sysMsg, solutions, solMsg, depVars, opts             *)
(* : Output  : Only printed output                                          *)
(* : Summary : Prints system with sysMsg as header, then prints solution    *)
(*           : with solMsg as the header                                    *)
(* ************************************************************************ *)
PDESolutionTesterPrintSystem[system_, sysMsg_, solutions_, solMsg_,
    depVars_, opts___?OptionQ]:=
    Module[
        { debugLevel  = DebugLevel /. {opts} /. Options[PDESolutionTester],
          verboseFlag = Verbose /. {opts} /. Options[PDESolutionTester]
        },

        (* ******************************************************** *)
        (* Prints the system in a more readable form                *)
        (* Code by Douglas Baldwin                                  *)
        (* ******************************************************** *)
        Print[sysMsg];
        Print[
            TableForm[system /. 
                {Derivative[n__][F_][x__]:>
                SequenceForm[F,
                    Subscript[SequenceForm @@ 
                        Flatten[
                            Table[#[[1]], {#[[2]]}
                            ]& /@ Transpose[{{x}, {n}}]
                        ]
                    ]
                ],
                Sequence @@ ((#->Head[#])& /@ depVars)
                }
            ]
        ];

        (* ******************************************************** *)
        (* Prints the solutions in a more readable form             *)
        (* ******************************************************** *)
        Print[solMsg];
        Print[
            TableForm[solutions /.
                {Sequence @@ ((#->Head[#])& /@ depVars)}
            ]
        ];

    ]; (* End Module *)

(* End of function PDESolutionTesterPrintSystem. *)

(* ************************************************************************ *)
(* : Title   : PDESolutionTesterCleanUp                                     *)
(* : Author  : Douglas Baldwin                                              *)
(* : Input   : Expression                                                   *)
(* : Output  : Expression                                                   *)
(* : Summary : Remove "PDESolutionTesterPrivate`" from output.              *)
(* ************************************************************************ *)

PDESolutionTesterCleanUp = 
  ToExpression[
    StringReplace[ToString[InputForm[#]],
      "Calculus`PDESolutionTester`Private`"->""
    ]
  ]&

(* End of function PDESolutionTesterCleanUp. *)

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

End[] (* `Private` context *)

Protect[PDESolutionTester]

SetAttributes[PDESolutionTester, ReadProtected]

EndPackage[]

Print["Package PDESolutionTester.m was successfully loaded."];

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


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

(* Nonlinear algebraic solver by Unal Goktas (WRI) and Willy Hereman *)
(* written by Unal Goktas in October/November 2000 *)

BeginPackage["Algebra`AnalyzeAndSolve`"]

Unprotect[AnalyzeAndSolve]

Begin["`Private`"]

If[$VersionNumber < 4,
   Attributes[Internal`DeactivateMessages] = {HoldAll};
   Internal`DeactivateMessages[body_, msgs___MessageName] :=
      Module[{wasOn = Select[Hold[msgs], Head[#] =!= $Off &], result},
         CheckAbort[
            Scan[Off, wasOn];
            result = body;
            Scan[On, wasOn];
            result,
            Scan[On, wasOn];
            Abort[]
         ]
      ]
]

AnalyzeAndSolve[system: {HoldPattern[_ == _]..}, primaryunknowns_,
   secondaryunknowns_, parameters_, nonzeros_] :=
   AnalyzeAndSolve[(#[[1]]-#[[2]]&) /@ system, primaryunknowns,
      secondaryunknowns, parameters, nonzeros]

AnalyzeAndSolve[system_, primaryunknowns_, secondaryunknowns_,{},nonzeros_]:=
   Block[{constraints},
      constraints = (And @@ ((# != 0 &) /@ nonzeros));
      Internal`DeactivateMessages[
         Solve[(And @@ ((# == 0 &) /@ system)) && constraints,
            Join[primaryunknowns, secondaryunknowns], Sort -> False],
         Solve::svars
      ]
   ]

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

(* "system" is polynomial in "primaryunknowns", "secondaryunknowns", and
"parameters". *)
AnalyzeAndSolve[system_, primaryunknowns_, secondaryunknowns_, parameters_,
   nonzeros_] :=
   Block[{a, globalsol = {}, constraints},
      a = Together[({#}& /@ system) /. {{}}];
      constraints = (And @@ ((# != 0 &) /@ nonzeros));
      MapThread[
         RecursivelyKeepTrack[#1, #2, primaryunknowns, secondaryunknowns,
            parameters, nonzeros, constraints]&, {a, {{}}}];
      Union[globalsol]
   ]

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

(* RecursivelyKeepTrack[{}, __] := {} *)

RecursivelyKeepTrack[{}, sol_, __] := (AppendTo[globalsol, sol]; {})

RecursivelyKeepTrack[system_, sol_, __] /;
   (!FreeQ[system, DirectedInfinity] || !FreeQ[system, Indeterminate]) := {}

RecursivelyKeepTrack[system_, sol_, primaryunknowns_, secondaryunknowns_,
   parameters_, nonzeros_, constraints_] :=
   Block[{a, b, c},
      a = FactorListAndCleanUp[#, primaryunknowns, secondaryunknowns,
             parameters, nonzeros]& /@ system;
      a = Union[a];
      a = Sort[a, (ComplexityLevel[#1, primaryunknowns, secondaryunknowns,
                      parameters] <=
                   ComplexityLevel[#2, primaryunknowns, secondaryunknowns,
                      parameters]&)];
      b = a[[1]];
      b = SolveStepByStep[#, primaryunknowns, secondaryunknowns, parameters,
             constraints /. sol, sol]& /@ b;
      b = (Sequence @@ # &) /@ b;
      If[Length[b] == 0, Return[{}]];
      b = Transpose[b];
      c = (constraints && (And @@ #))& /@ b[[2]];
      b = b[[1]];
      a = Together[
             Internal`DeactivateMessages[Rest[a] /. b,
                Power::infy, General::indet
             ]
          ];
      a = Numerator[a];
      a = DeleteCases[a, {___, 0, ___}, {2}];
      MapThread[
         RecursivelyKeepTrack[#1, #2, primaryunknowns, secondaryunknowns,
            parameters, nonzeros, #3]&, {a, b, c}]
   ]

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

FactorListAndCleanUp[
   system_, primaryunknowns_, secondaryunknowns_, parameters_, nonzeros_] :=
   Block[{a},
      a = Flatten[Map[First[#]&, FactorList /@ system, {2}]];
      a = DeleteCases[a, _?(NumericQ[#] || MemberQ[nonzeros, #]&)];
      Union[a]
   ]

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

ComplexityLevel[expr_List, primaryunknowns_,secondaryunknowns_,parameters_]:=
   Max[ComplexityLevel[#, primaryunknowns, secondaryunknowns, parameters]& /@
      expr]

ComplexityLevel[expr_, primaryunknowns_, secondaryunknowns_, parameters_] :=
   Block[{a, b},
      a = Union[Cases[{expr}, q_ /; MemberQ[primaryunknowns, q], -1]];
      b = Length[a];
      (
       SubComplexityLevel[expr, a, b, 1]
      ) /; b != 0
   ]

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

ComplexityLevel[expr_, primaryunknowns_, secondaryunknowns_, parameters_] :=
   Block[{a, b},
      a = Union[Cases[{expr}, q_ /; MemberQ[secondaryunknowns, q], -1]];
      b = Length[a];
      (
       SubComplexityLevel[expr, a, b, 2] 
      ) /; b != 0
   ]

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

ComplexityLevel[expr_, primaryunknowns_, secondaryunknowns_, parameters_] :=
   Block[{a, b},
      a = Union[Cases[{expr}, q_ /; MemberQ[parameters, q], -1]];
      b = Length[a];
      (
       SubComplexityLevel[expr, a, b, 3] 
      ) /; b != 0
   ]

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

SubComplexityLevel[expr_, a_, b_, level_] :=
   Block[{c = Select[a, PolynomialQ[expr, #]&]},
      Which[
         Length[c] == 0,
         100*b*level+LeafCount[expr],
         True,
         Min[Exponent[expr, c]]
      ]
   ]

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

(* solves a single equation *)

SolveStepByStep[eqn_, primaryunknowns_, secondaryunknowns_, parameters_,
   constraints_, sol_] :=
   Block[{a},
      a = Union[Cases[{eqn}, q_ /; MemberQ[primaryunknowns, q], -1]];
      (
       SubSolveStepByStep[eqn, a, 
          Flatten[{secondaryunknowns, parameters}],
          constraints, sol
       ]
      ) /; Length[a] != 0
   ]

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

SolveStepByStep[eqn_, primaryunknowns_, secondaryunknowns_, parameters_,
   constraints_, sol_] :=
   Block[{a},
      a = Union[Cases[{eqn}, q_ /; MemberQ[secondaryunknowns, q], -1]];
      (
       SubSolveStepByStep[eqn, a, parameters, constraints, sol]
      ) /; Length[a] != 0
   ]

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

SolveStepByStep[eqn_, primaryunknowns_, secondaryunknowns_, parameters_,
   constraints_, sol_] :=
   Block[{a},
      a = Union[Cases[{eqn}, q_ /; MemberQ[parameters, q], -1]];
      (
       SubSolveStepByStep[eqn, a, {}, constraints, sol]
      ) /; Length[a] != 0
   ]

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

SolveStepByStep[___] := {}

SubSolveStepByStep[eqn_, primaryunknowns_, pars_, constraints_, sol_] :=
   Block[{a, b, c},
      a = Select[primaryunknowns, PolynomialQ[eqn, #]&];
      If[Length[a] != 0,
         a = Sort[a, Exponent[eqn, #1] <= Exponent[eqn, #2]&];
         a = Flatten[{a, Complement[primaryunknowns, a]}],
         a = primaryunknowns
      ];
      c = Reduce[eqn == 0 && constraints, Flatten[{a, pars}], Sort -> False];
      a = {ToRules[c]};
      If[Length[a] == 0, Return[{}]];
      c = Cases[{#}, p_ != q_, -1]& /@ If[Head[c] === Or, List @@ c, {c}];
      b = Internal`DeactivateMessages[sol /. a, Power::infy, General::indet];
      Union[Transpose[{Flatten /@ Thread[{b, a}, List, 2], c}]]
   ]

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

End[] (* `Private` context *)

SetAttributes[AnalyzeAndSolve, ReadProtected]

EndPackage[]

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