   
(*****************************************************************************)
(*                                                                           *)
(*          *** M A T H E M A T I C A   P R O G R A M ***                    *)
(*                                                                           *)
(*      SYMBOLIC COMPUTATION of CONSERVED DENSITIES for SYSTEMS of           *)
(*                  NONLINEAR EVOLUTION EQUATIONS                            *)
(*                                                                           *)
(* program name: condens.m                                                   *)
(*                                                                           *)
(* purpose: computation of the conserved densities and corresponding flux    *)
(*          with possible compatibility conditions, verification of the      *)
(*          conservation law.                                                *)
(*                                                                           *)
(* input to condens.m : system of nonlinear evolution equations of any order,*)
(*                      any degree, polynomial type, in variables x and t,   *)
(*                      only constant parameters, functions in x and t are   *)
(*                      NOT allowed as parameters                            *)
(*                                                                           *)
(*                      u[i]_t = f(u[1],...,u[N],u[1]_{nx},...,u[N]_{nx})    *)
(*                      with i=1,...,N; and n=1,2,...                        *)
(*                                                                           *)
(* output : density and flux of desired rank (if it exists),                 *)
(*          and compatibility conditions on parameters, if applicable        *)
(*                                                                           *)
(* tested on : IBM RISC 6000, IBM Compatible PC 486, SGI Indigo2 XL          *)
(*                                                                           *)
(* language : Mathematica 3.0 and 2.2 (also versions 2.0 and 2.1)            *)
(*                                                                           *)
(* authors: Unal Goktas and Willy Hereman                                    *)
(*          Department of Mathematical and Computer Sciences                 *)
(*          Colorado School of Mines                                         *)
(*          Golden, CO 80401-1887, USA                                       *)
(*                                                                           *)
(* Version 3.2: May 4, 1998                                                  *)
(*                                                                           *)
(* Copyright 1998                                                            *)
(*                                                                           *)
(*****************************************************************************)

Clear["Global`*"]; 

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

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

commentinter[] := Block[{},
Print["*********************************************************"];
Print["        WELCOME TO THE MATHEMATICA PROGRAM               "];
Print["         by UNAL GOKTAS and WILLY HEREMAN                "];
Print["    FOR THE COMPUTATION OF CONSERVED DENSITIES.          "];
Print["      Version 3.2 released on May 4, 1998                "];
Print["                 Copyright 1998                          "];
Print["*********************************************************"]
]; (* end commentinter *)

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

(*****************************************************************************)
(* cls: clears the screen                                                    *)
(*****************************************************************************)

cls := Print["\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n\n"];

(*****************************************************************************)
(* printpage[ , ]: a subroutine for menu                                     *)
(*****************************************************************************)
 
printpage[n_,page_] := Module[
{i,choice,lenpage},
cls; 
lenpage = Length[page[n]];
Print[" "];
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 *)

(*****************************************************************************)
(* menu: creates the menu                                                    *)
(*****************************************************************************)

menu := Module[
{counterpage = 1,menulist,numpages,page,
choice1,control,choice2,choice3,lenmenulist,i},
menulist = {
{"  1) KdV Equation (d_kdv.m)"},
{"  2) Modified KdV Equation (d_mkdv.m)"},
{"  3) Fifth Order KdV Equation-parameterized (d_5kdv.m)"},
{"  4) Fifth Order KdV Equation-Lax Case (d_5lax.m)"},
{"  5) Fifth Order KdV Equation-SK Case (d_5sk.m)"},
{"  6) Fifth Order KdV Equation-KK Case (d_5kk.m)"},
{"  7) Fifth Order KdV Equation-Ito Case (d_5ito.m)"},
{"  8) Seventh Order KdV Equation-parameterized (d_7kdv.m)"},
{"  9) Seventh Order KdV Equation-Lax Case (d_7lax.m)"},
{" 10) Seventh Order KdV Equation-SK-Ito Case (d_7ski.m)"},
{" 11) Seventh Order KdV Equation KK Case (d_7kk.m)"},
{" 12) Schamel Equation (d_scham.m)"},
{" 13) Rosenau-Hyman Equation (d_roshy.m)"},
{" 14) Hirota-Satsuma System-parameterized (d_phrsat.m)"},
{" 15) Hirota-Satsuma System (d_hirsat.m)"},
{" 16) Ito System-parameterized (d_pito.m)"},
{" 17) Ito System (d_ito.m)"},
{" 18) NLS System (real and imaginary parts) (d_nls.m)"},
{" 19) NLS System (equation and conjugate) (d_nlsalt.m)"},
{" 20) DNLS System (d_dnls.m)"},
{" 21) MVDNLS System (d_mvdnls.m)"},
{" 22) Kaup System-parameterized (d_pkaup.m)"},
{" 23) Kaup System (d_kaup.m)"},
{" 24) Kaup-Broer System (d_broer.m)"},
{" 25) Drinfel'd-Sokolov System (d_soko.m)"},
{" 26) Dispersive Long Wave System (d_disper.m)"}, 
{" 27) Non-dispersive Long Wave System (d_nodisp.m)"}, 
{" 28) 3-Component KdV System (d_3ckdv.m)"},
{" 29) 2-Component Nonlinear Schrodinger Equation (d_2cnls.m)"},
{" 30) 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}]]]]];

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 we supplied)."];
   choice2 = Input["If your 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,<<d_kdv.m,2,<<d_mkdv.m,3,<<d_5kdv.m,4,<<d_5lax.m,5,<<d_5sk.m,
6,<<d_5kk.m,7,<<d_5ito.m,8,<<d_7kdv.m,9,<<d_7lax.m,10,<<d_7ski.m,
11,<<d_7kk.m,12,<<d_scham.m,13,<<d_roshy.m,14,<<d_phrsat.m,15,<<d_hirsat.m,
16,<<d_pito.m,17,<<d_ito.m,18,<<d_nls.m,19,<<d_nlsalt.m,20,<<d_dnls.m,
21,<<d_mvdnls.m,22,<<d_pkaup.m,23,<<d_kaup.m,24,<<d_broer.m,25,<<d_soko.m,
26,<<d_disper.m,27,<<d_nodisp.m,28,<<d_3ckdv.m,29,<<d_2cnls.m,30,<<d_bous.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 *)

(*****************************************************************************)
(* EulerD[  ]: applies the Euler operator for variational derivatives.       *)
(*             Written by Yu He in 1992, simplified by W. Hereman            *)
(*****************************************************************************)

EulerD[f_, (y_)[x_, r___], w:{x_, r___}] := 
  Module[{Dfuncs, Dtimes, dummyfunc}, 
    Dfuncs = Union[Cases[{f}, Derivative[__][y][__], Infinity]]; 
    Dtimes = (Head[Head[#1]] & ) /@ Dfuncs /. Derivative -> List; 
    Expand[D[f, y[x, r]] + (ReleaseHold[Thread[dummyfunc[(D[f, 
                 #1] & ) /@ Dfuncs, 
             (Hold[Apply[Sequence, #1]] & ) /@ 
              (Thread[{w, #1}] & ) /@ Dtimes]]] /. dummyfunc -> D) . 
        ((-1)^#1 & ) /@ (Apply[Plus, #1] & ) /@ Dtimes] ]; (* end EulerD *)
 
EulerD[f_, v:{(y_)[x_, r___], ___}, w:{x_, r___}] := 
   (EulerD[f, #1, w] & ) /@ v /; 
   If[Apply[And, (MatchQ[#1, _[Apply[Sequence, w]]] & ) /@ v], 
   True, Message[EulerD::argx, w]];

EulerD[f_, (y_)[x_], x_] := EulerD[f, y[x], {x}];
EulerD[f_, v:{(y_)[x_], ___}, x_] := EulerD[f, v, {x}];

(*****************************************************************************)
(* pdeform[]: takes an expression and 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.                                         *)
(*****************************************************************************)

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

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

(*****************************************************************************)
(* split[]: if free coefficients c[i] remain in the density, the density     *)
(*          is further split in independent pieces corresponding to each     *)
(*          coefficient c[i]                                                  *)
(*****************************************************************************)

split[expr_] := Module[{listofc,lenlistofc,ci,rhoi,rest}, 
    rest=Expand[expr];
    If[FreeQ[rest,c[__]], 
    Print[" "];
    Print["The density has no free coefficients."], 
    listofc=Union[Cases[rest,c[__],12]];
    lenlistofc=Length[listofc];
    Print[" "];
    Print["There is/are ", lenlistofc, " free coefficient(s) in the density."];
    Print["These free coefficients are: ", listofc, "."];
    Print["Splitting the density in independent pieces."];
     Do[
       ci=Part[listofc,i];
       rhoi=Expand[Coefficient[rest,ci]]; 
       rest=Expand[rest-ci*rhoi];
       Print[" "];
       Print["Part of the density with coefficient ", ci, ":"];
       Print[" "];
       Print[rhoi],{i,1,lenlistofc}]; (* end do *)
     Print[" "];
     Print["Part of the density that has no free coefficient:"];
     Print[" "];
     Print[rest];
     Print[" "] ] (* end if *) ]; (* end module *)

(*****************************************************************************)
(* normal[]: normalizes the density according to the coefficient of the      *)
(*           first term in rho[x,t], returns normalized forms of rho[x,t]    *)
(*           and j[x,t]. Does not divide by the coefficient if the           *)
(*           coefficient is not a number. Returns a list with rho in the     *)
(*           first place, and j in the second place. To be used on the list  *)
(*           that comes out of stripper[rho[x,t]]                            *)
(*****************************************************************************)

normal[mylist_List] := Module[{leadingcoef,facexpr,denrho,rhojlist}, 
(* Cancel common denominators first, keep track of these with parameters *)
facexpr = Factor[rho[x,t]];
denrho = Denominator[facexpr];
If[Not[NumberQ[denrho]], denrho = Factor[denrho]];
If[Not[FreeQ[Head[denrho],Times]],
   If[NumberQ[Part[denrho,1]] && Part[denrho,1] =!= 0, 
     denrho = denrho/Part[denrho,1]]];
If[Not[NumberQ[denrho]],  
  Print["WARNING!"];
  Print[" "];
  Print["Both rho[x,t] and j[x,t] were multiplied by: ",denrho,"."]; 
  Print[" "];
  rho[x,t] = Expand[Factor[rho[x,t]*denrho]]; 
  j[x,t] = Expand[Factor[j[x,t]*denrho]]];
leadingcoef = Coefficient[rho[x,t],Part[mylist,1]];
If[NumberQ[leadingcoef] && leadingcoef =!= 0,
  rho[x,t] = Expand[rho[x,t]/leadingcoef];
  j[x,t] = Expand[j[x,t]/leadingcoef]]; 
rhojlist = {rho[x,t],j[x,t]};
Return[rhojlist]; ]; (* end normal *)

(*****************************************************************************)
(* stripper[]: cancels the numerical factors of each term of the argument    *)
(*             and puts these into a list after cancellation                 *)
(*****************************************************************************)

stripper[givenexpr_] := Module[
{lenexpr1,list = {},expr1,expr2,expr3,expr3s,lenexpr3,i,k,s},
expr1 = Expand[givenexpr];
If[Head[expr1] === Plus,lenexpr1 = Length[expr1],lenexpr1 = 1];
For[k = 1,k <= lenexpr1,k++,
   If[lenexpr1 == 1,expr2 = expr1,expr2 = Part[expr1,k]];
   expr2 = expr2 /. c[i_] :> 1;
   expr3 = FactorList[expr2];
   lenexpr3 = Length[expr3];
   For[s = 1,s <= lenexpr3,s++,
      expr3s = Part[expr3,s];
      If[FreeQ[expr3s,u],
        If[FreeQ[weightpars,Part[expr3s,1]],
          expr2 = expr2/Part[expr3s,1]^Part[expr3s,2];
          ],
        s = lenexpr3+1;
        ];
      ];
   list = Append[list,expr2];
   ];
Return[list]; ]; (* end stripper *)

(*****************************************************************************)
(* highestpair[]: a subroutine for highest order function                     *)
(*****************************************************************************)

highestpair[pairlist_] := Module[
{length,k,highestlist},
If[pairlist != {},
  length = Length[pairlist];
  highestlist = {Part[pairlist,1]};
  If[length >= 2,
    For[k = 2,k <= length,k++,
      If[Part[pairlist,k][[1]] > Part[highestlist,Length[highestlist]][[1]],
        highestlist = {Part[pairlist,k]},
        If[Part[pairlist,k][[1]] === Part[highestlist,Length[highestlist]][[1]],
	  highestlist = Union[Append[highestlist,pairlist[[k]]]]
          ];
        ];
      ];
   ];
 Return[highestlist],
 If[debug,
   Print["The expression (given to the function that computes the"]; 
   Print["highest order) does not have any derivatives!"]];
 Return[{{0,0}}];
 ]  ]; (* end highestpair *)
 
(*****************************************************************************)
(* highestorder[]: returns the highest order (or the term with the highest   *)
(*                 order) of the argument                                    *)
(*****************************************************************************)
 
highestorder[gexpr_,debug11_,debug12_,debug21_,debug22_,debug31_,debug32_] := 
Module[
{length1,i,newexpr1,length2,newexpr6,temporder1,temporder2,newexpr5,k,
newexpr2,newexpr3,temporder3,newexpr4,orderlist1 ={},orderlist2 = {},
orderlist3 = {},highestlist1,highestlist2,highestlist3,expr},
expr = Expand[gexpr];
If[Not[FreeQ[expr,Derivative]], 
  length1 = length[expr];
  temporder1 = {}; temporder2 = {}; temporder3 = {}; newexpr5 = 1;
  For[i = 1,i <= length1,i++,
     If[Not[FreeQ[expr,Plus]],temporder1 = {}; temporder2 = {};
            temporder3 = {};newexpr5 = 1];
     If[length1 != 1,newexpr1 = expr[[i]],newexpr1 = expr]; 
       If[Not[FreeQ[newexpr1,Derivative]],  
	 length2 = length[newexpr1];
	 For[k = 1,k <= length2,k++,
	    If[length2 != 1,newexpr2 = newexpr1[[k]],newexpr2 = newexpr1];
	    If[Not[FreeQ[newexpr2,Derivative]],
	      If[FreeQ[newexpr2,Power],
		newexpr3 = FullForm[newexpr2];
                newexpr6 = newexpr3[[1]][[0]][[0]];
		temporder1 = Append[temporder1,newexpr6[[1]]];
		temporder2 = Append[temporder2,newexpr6[[2]]];
		temporder3 = Append[temporder3,newexpr6[[1]]+newexpr6[[2]]];
		newexpr4 = newexpr2;
		newexpr5 = newexpr5*newexpr2, 
		newexpr3 = FullForm[newexpr2];
                newexpr6 = newexpr3[[1]][[1]][[0]][[0]];  
		temporder1 = Append[temporder1,newexpr6[[1]]]; 
		temporder2 = Append[temporder2,newexpr6[[2]]];
		temporder3 = Append[temporder3,newexpr6[[1]]+newexpr6[[2]]];
		newexpr4 = newexpr3[[1]];
		newexpr5 = newexpr5*newexpr3[[1]];
		];  
	      ];  
	    ];  
         temporder1 = Max[temporder1];
         temporder2 = Max[temporder2];
         temporder3 = Max[temporder3];
	 orderlist1 = Append[orderlist1,{temporder1,newexpr4}];
	 orderlist2 = Append[orderlist2,{temporder2,newexpr4}];
	 orderlist3 = Append[orderlist3,{temporder3,newexpr5}];
	 ];  
     ];        
  ];                  
Which[
debug11,highestlist1 = highestpair[orderlist1];
  If[debug12,Return[Table[highestlist1[[k]][[2]],{k,1,Length[highestlist1]}]],
     Return[highestlist1[[1]][[1]]]], 
debug21,highestlist2 = highestpair[orderlist2];
  If[debug22,Return[Table[highestlist2[[k]][[2]],{k,1,Length[highestlist2]}]],
      Return[highestlist2[[1]][[1]]]],
debug31,highestlist3 = highestpair[orderlist3];
  If[debug32,Return[Table[highestlist3[[k]][[2]],{k,1,Length[highestlist3]}]],
      Return[highestlist3[[1]][[1]]]] ]; ]; (* end highestorder *)

(*****************************************************************************)
(* setuniformrank[]: a subroutine for the scaling function                   *)
(*****************************************************************************)

setuniformrank[list_] := Module[
{i,lengthlist,syslist={}},
lengthlist = Length[list];
For[i = 1,i < lengthlist,i++,
   syslist = Append[syslist,Part[list,i] == Part[list,lengthlist]]; 
   ];
Return[syslist]; ]; (* end setuniformrank *)

(*****************************************************************************)
(* scaling[]: determines the scaling properties, if the equations are        *)
(*            incompatible it will print the appropriate messages            *)
(*****************************************************************************)

scaling[eqlist_List] := Module[
{pointsym,expr,list,syslist,msyslist = {},i,j,k,scalesol,
lenmsyslist,tempscalesol = {},tempmsyslist,trouble,troubleleft,
troubleright,posleft,posright,troublelist},
If[weightpars =!= {},
  pointsym = Union[Table[weightpars[[i]] -> E^weight[weightpars[[i]]],
                     {i,1,Length[weightpars]}],
       {u[i_Integer][x,t] :> E^weightu[i],
        D[u[i_Integer][x,t],{x,j_Integer}] :> E^(weightu[i]+j),
        D[u[i_Integer][x,t],{t,j_Integer}] :> E^(weightu[i]+j*weight[d/dt]),
        x -> E, t -> E^weight[d/dt]}],
  pointsym = {u[i_Integer][x,t] :> E^weightu[i],
        D[u[i_Integer][x,t],{x,j_Integer}] :> E^(weightu[i]+j),
        D[u[i_Integer][x,t],{t,j_Integer}] :> E^(weightu[i]+j*weight[d/dt]),
        x -> E, t -> E^weight[d/dt]};
  ];
For[i = 1, i <= noeqs, i++,
   expr[i] = Part[eqlist,i];
   list[i] = stripper[expr[i]];
   list[i] = list[i] /. pointsym;
   list[i] = PowerExpand[Log[list[i]]];
   syslist[i] = setuniformrank[list[i]];
   msyslist = Union[msyslist,syslist[i]];
   ];
scalesol = Flatten[Solve[msyslist]];
If[MemberQ[msyslist,False],
  Print[" "];
  Print["Fatal Error! Weights are incorrect."];
  Print["Check your choice(s) for the weights in the data file."];
  Print["Aborting the computations!"];
  CloseLog[]; Abort[],
  If[scalesol === {} && msyslist =!= {True},
    Print[" "];
    Print["In the given system there is at least one equation with terms"];
    Print["of unequal rank. Scaling properties can not be determined for"];
    Print["this system. The program will try to find the conflict, and,"];
    Print["if successful, provide suggestions to help resolve the conflict."];
    lenmsyslist = Length[msyslist];
    i = 1;
    While[tempscalesol === {} && i <= lenmsyslist,
         tempmsyslist = Complement[msyslist,{Part[msyslist,i]}];
         tempscalesol = Flatten[Solve[tempmsyslist]];
         i++;
         ];
    i--;
    If[tempscalesol =!= {},
      trouble = Part[msyslist,i];
      troubleleft = Part[trouble,1];
      troubleright = Part[trouble,2];
      For[j = 1, j <= noeqs, j++,
         posleft = Flatten[Position[list[j],troubleleft,{1}]];
         posright = Flatten[Position[list[j],troubleright,{1}]];
         If[posleft =!= {} && posright =!= {},
           troublelist = Union[
             Table[Part[Expand[eq[j][x,t]],posright[[k]]],
                     {k,1,Length[posright]}],
             Table[Part[Expand[eq[j][x,t]],posleft[[k]]],
                     {k,1,Length[posleft]}]
                ];
           Print["The terms"];
           Print[troublelist];
           Print["in equation ",j," are incompatible. Try to introduce an"];
           Print["auxiliary parameter with weight as coefficient of one of"];
           Print["these terms. Aborting the computations! "];
           ];
         ],
    Print["Try to introduce auxiliary parameters with weight as coefficients"];
    Print["into the system. Aborting the computations!"];
    ];
    CloseLog[];
    Abort[],
    Return[scalesol];
    ];
  ]; ]; (* end scaling *)

(*****************************************************************************)
(* subroutine4[]: a subroutine for the function evaluate                     *)
(*****************************************************************************)

subroutine4[givenexpr_] := Module[
{expr1,nexpr1,freerule = {},part1,lennexpr1,i},
expr1 = Expand[givenexpr];
nexpr1 = Numerator[Factor[expr1]];
lennexpr1 = length[nexpr1];
For[i = 1,i<= lennexpr1,i++,
   part1 = Part[nexpr1,i];
   If[MemberQ[unknownlist,part1],
     freerule = Union[freerule,{part1 -> 1}]];
   ];
If[freerule =!= {},
  Print[" "];
  Print["Caution! There is a common (free) factor in the density."];
  Print[" "];
  Print["Setting ",freerule,"."];
  ]; 
Return[freerule]; ]; (* end subroutine4 *)

(*****************************************************************************)
(* subroutine5[]: a subroutine for the function constructformrho             *)
(*****************************************************************************)

subroutine5[nodims_,givenlist_] := Module[
{m,boundary,scale,var,len,tempwei,temppar,list,pair,ctm = 0,k,s},
list[0] = {};
For[m = 1,m <= nodims,m++,
   pair = Part[givenlist,m];
   var = Part[pair,1];
   scale = Part[pair,2];
   If[scale =!= 0,
     If[m === 1,
       boundary = Floor[rhorank/scale];
       list[m] = Table[{i*scale,var^i},{i,0,boundary}];
       ctm = m,
       len = Length[list[m-1]];
       list[m] = {};
       For[k = 1,k <= len,k++,
          tempwei = Part[Part[list[m-1],k],1];
          temppar = Part[Part[list[m-1],k],2];
          boundary = Floor[(rhorank-tempwei)/scale];
          For[s = 0,s <= boundary,s++,
             list[m] = Union[list[m],{{tempwei+s*scale,temppar*var^s}}];
             ctm = m;
             ];
          ];
       ];
     ];
   ];
Return[list[ctm]]; ]; (* end subroutine5 *)

(*****************************************************************************)
(* constructformrho[]: returns the form of rho                               *)
(*****************************************************************************)
 
constructformrho[nodims_,varscalelist_] := Module[
{list,len,tempwei,temppar,inttest,formrholist = {},triviallist = {},
difflist1,difflist2,lendifflist1,i,j}, 
Print[" "];
Print["Starting the construction of the form of the density."];
list = subroutine5[nodims,varscalelist];
len = Length[list];
For[i = 1,i <= len,i++,
   tempwei = Part[Part[list,i],1];
   temppar = Part[Part[list,i],2];
   inttest = rhorank-tempwei;
   If[IntegerQ[inttest] && Not[FreeQ[temppar,u]],
     If[inttest === 0, 
       formrholist = Union[formrholist,{temppar}],
       difflist1 = stripper[D[temppar,{x,inttest-1}]];
       lendifflist1 = Length[difflist1];
       For[j = 1, j <= lendifflist1,j++,
          difflist2 = stripper[D[Part[difflist1,j],x]];
          triviallist = Union[triviallist,{Last[difflist2]}];
          formrholist = Union[formrholist,Complement[difflist2,triviallist]]; 
          ];
       ];
     ];
   ];

lenformrho = Length[formrholist];
formrho[x,t] = 0;
For[i = 1,i <= lenformrho,i++,
   formrho[x,t] = formrho[x,t]+c[i]*Part[formrholist,i];
   ];
Print[" "];
Print["For RANK = ",rhorank,", this is the form of density rho[x,t]: "];
Print[" "];
Print[pdeform[Expand[formrho[x,t]]]];
Print[" "];

If[formrho[x,t] === 0,
  Print["The only density is the trivial density rho[x,t] = 0."];
  Print[" "];
  Print["Aborting the computations! "]; 
  CloseLog[]; Abort[] ] ]; (* end constructformrho *)               

(*****************************************************************************)
(* evaluate[]: computes the density and flux, with verification              *)
(*****************************************************************************)

evaluate[expr1_,expr2_,solc_,rules_] := Module[
{test,factortest,newrhot,newformj,freerule,rhojlist}, 
rho[x,t] = expr1 /. solc;
rho[x,t] = Expand[rho[x,t]];
newrhot[x,t] = expr2 /. Flatten[{rules,solc}]; 
newrhot[x,t] = Expand[newrhot[x,t]];
newformj[x,t] = (-1)*Integrate[newrhot[x,t],x];
j[x,t] = Expand[newformj[x,t]];
Print[" "];
If[FreeQ[j[x,t],Integrate] && rho[x,t] =!= 0,
  If[rho[x,t] =!= 0,freerule = subroutine4[rho[x,t]],freerule = {}];
  rho[x,t] = rho[x,t] /. freerule;
  j[x,t] = j[x,t] /. freerule;
  test = Expand[D[rho[x,t],t]+D[j[x,t],x] /. rules];
  factortest = Expand[test];
  If[factortest =!= 0,factortest = Factor[factortest]];
  If[factortest =!= 0,Print[" "];
    Print["Automatic verification of the conservation law FAILED!"]];
  Print[" "];
  rhojlist = normal[stripper[rho[x,t]]];
  rho[x,t] = Part[rhojlist,1];
  rho[x,t] = Factor[rho[x,t]];
  rho[x,t] = Expand[rho[x,t]]; 
  j[x,t] = Part[rhojlist,2];
  j[x,t] = Factor[j[x,t]];
  j[x,t] = Expand[j[x,t]];
  Print["*******************************************************"];
  Print[" "];
  Print["The normalized density rho[x,t] is:"];
  Print[" "];
  Print[pdeform[rho[x,t]]];
  split[pdeform[rho[x,t]]];
  Print[" "];
  If[Length[j[x,t]] <= 20 || debug,
    Print["*******************************************************"];
    Print[" "];
    Print["The corresponding flux j[x,t] is:"];
    Print[" "];
    Print[pdeform[j[x,t]]],
    Print["*******************************************************"];
    Print[" "];
    Print["Flux j[x,t] has been computed. It has ",Length[j[x,t]]," terms."];
    Print["Since there are more than 20 terms, j[x,t] will not be shown."];
    Print["To see the flux at the end, type: pdeform[j[x,t]] or j[x,t]"]
    ];
  Print[" "];
  Print["*******************************************************"];
  Print["Result of explicit verification: (rho_t + J_x) = ",factortest];
  Print["*******************************************************"], 
  Print["*******************************************************"];
  Print["There is NO DENSITY of this form!"]; 
  Print["*******************************************************"];
  ] ]; (* end evaluate *) 

(*****************************************************************************)
(* constructeqlist[]: forms the system of equations for coefficients c[i]    *)
(*****************************************************************************)

constructeqlist[expr_] := Module[
{expr1,rules1,rules2,eqlist = {},tempexpr,expr2,templist,lentemplist,s,
tempcoef,parttemplist,expr3},
rules1 = D[u[i_][x,t],{x,j_}] :> 0;
rules2 = u[i_][x,t] :> 0;
expr1 = Expand[expr];
tempexpr = expr1 //. rules1;
tempexpr = tempexpr //. rules2;
eqlist = Union[eqlist,{tempexpr}]; 
expr2 = expr1-tempexpr;
If[expr2 =!= 0,
  templist = Union[stripper[expr2]];
  lentemplist = Length[templist];
  expr3 = expr2;
  For[s = 1,s <= lentemplist,s++,
     parttemplist = Part[templist,s];
     tempcoef = Coefficient[Expand[expr2],parttemplist];
     tempcoef = tempcoef //. rules1;
     tempcoef = tempcoef //. rules2;
     If[debug,expr3 = expr3 - (tempcoef*parttemplist)];
     eqlist = Union[eqlist,{tempcoef}];
     ];
  If[debug,
    expr3 = Factor[expr3];
    Print["Remainder of the expression during the construction of"];
    Print["the system for the coefficients c[i]: ",expr3]];
  If[debug && expr3 =!= 0,
    Print["Error in the construction of the system for the c[i]!"];
    Print["Aborting the computations!"];
    CloseLog[]; Abort[] ];
  ];
Return[eqlist]; ]; (* end constructeqlist *)

(*****************************************************************************)
(* mysimplify[]: simplifies the system of equations for the c[i]             *)
(*****************************************************************************)
 
mysimplify[list_] := Module[
{newlist,lennewlist,k,fterm,simplelist = {}},
Print["Starting the simplification of the system."];
Print[" "];
newlist = Complement[list,{0}];
lennewlist = Length[newlist];
For[k = 1 ,k <= lennewlist,k++,
   fterm = Factor[newlist[[k]]];
   If[Part[fterm,0] === Times && NumberQ[Part[fterm,1]],  
     fterm = fterm/Part[fterm,1] ];
   If[Intersection[Expand[simplelist],Expand[{fterm*(-1)}]] === {},
     simplelist = Union[simplelist,{fterm}]];
   ];
Return[simplelist]; ]; (* end mysimplify *)

(*****************************************************************************)
(* analyzer[]: determines the coefficients c[i] that must be zero, used in   *)
(*             the search for compatibility conditions                       *)
(*****************************************************************************)

analyzer[list_,len_] := Module[
{i,partlist,analyzelist = {}},
Print[" "];
Print["Starting the analysis of the system."];
Print[" "];
For[i = 1,i <= len,i++,
   partlist = Part[list,i];
   If[Length[partlist] === 1 && MemberQ[unknownlist,partlist], 
     analyzelist = Union[analyzelist,{partlist}],
     If[Length[partlist] === 2 && MemberQ[unknownlist,partlist[[2]]] &&
      (MemberQ[parameters,partlist[[1]]] || MemberQ[weightpars,partlist[[1]]]),
      analyzelist = Union[analyzelist,{partlist[[2]]}] ] ];
   ];
Return[analyzelist]; ]; (* end analyzer *)

(*****************************************************************************)
(* picklhs[]: makes a system, setting LHS == 0                               *)
(*****************************************************************************)

picklhs[list_,k_] := Part[list,k] == 0;

(*****************************************************************************)
(* coefmat[]: creates the coefficient matrix of the system for the c[i]      *)
(*****************************************************************************)

coefmat[system_List,unknowns_List] := Module[
{i,mat = {},lensys}, 
If[system === {} || unknowns === {}, Print["Fatal Error"],
lensys = Length[system];
For[i = 1,i <= lensys,i++, 
   mat = Append[mat,Table[Coefficient[Expand[Part[Part[system,i],1]],
   Part[unknowns,k]],{k,1,Length[unknowns]}]]]];
Return[mat]; ]; (* end coefmat *)

(*****************************************************************************)
(* main[]: solves the system for the c[i], computes the compatibility        *)
(*         conditions (if applicable)                                        *)
(*****************************************************************************)

main[utlist_] := Module[
{rules = {},trules,ressubst,formj,tformj,dtformj,comcond,k,comcondfactab,
myeqlist,maineqlist,lengthpar,syscondlist1={},syscond,inputval,inputrule,sol1,
lengthsol1,nrules1,par,newmaineqlist,solc,nrules2,nsolc,funclist,j,seqlist,
lenseqlist,analyzelist,inputpart,inputlist,i,iii,complexanalysis,comcondfac},
rhseq[i_] := D[u[i][x,t],t] /. Part[utlist,i];
rule[i_,n_] := D[D[u[i][x,t],{x,n}],{t,1}] -> D[rhseq[i],{x,n}];
For[i = 1,i <= noeqs,i++,
   trules[i] = Table[rule[i,n],{n,0,highestorderrho}];
   rules = Union[rules,trules[i]];
   ];
If[debug,
  Print["These are the replacement rules from the given system: ",rules]];
(* Computation of j *)
rhot[x,t] = D[formrho[x,t],t];
ressubst = rhot[x,t] /. rules;
ressubst = Expand[-ressubst]; 
If[debug,Print["The Euler operator will be applied to: ",
  pdeform[ressubst]]]; 
funclist = Table[u[i][x,t],{i,1,noeqs}];
Print[" "];
Print["Applying the Euler operator (the variational derivative)."];
formj = EulerD[ressubst,funclist,{x,t}];
Print[" "];
Print["Finished with variational differentiation."];
tformj = formj;
(* begin If point 1 *)
If[Length[Cases[tformj,x_0]] === noeqs,
  Print[" "];
  Print["End of the computations."];
  evaluate[formrho[x,t],rhot[x,t],{},rules], (* else point 1 *)
  dtformj = tformj; 
  dtformj = Expand[dtformj];
  Print[" "];
  Print["Starting the construction of the system for the c[i]."];
  Print[" "];
  myeqlist = {};
  For[j = 1,j <= noeqs,j ++,
     If[Part[dtformj,j] =!= 0,
       myeqlist = Union[myeqlist,constructeqlist[Part[dtformj,j]]]];
     ];
 seqlist = mysimplify[myeqlist];
 lenseqlist = Length[seqlist];
 maineqlist = Flatten[Table[picklhs[seqlist,i],{i,1,lenseqlist}]];
 system = MapAll[Factor,maineqlist];
If[
(Not[debug] && Length[maineqlist]>50) && (parameters==={} && weightpars==={}), 
 Print["Since there are no parameters, and the linear system"];
 Print["for the c[i] has ",Length[maineqlist],
       " equations, which is more than 50,"];
 Print["the system will not be shown."];
 Print["To see the system at the end of the computations, type: system"];
 Print[" "],
 Print["This is the system for the c[i]:"]; 
 Print[" "];
 Print[system];
 Print[" "];
 Print["List of unknown coefficients c[i]:"];
 Print[" "];
 Print[unknownlist];
 Print[" "] 
  ];
 Print["Starting the solution of this system with ",
 Length[maineqlist], " equations."];
 If[parameters =!= {}, (* then start for parameter test *)
   analyzelist = analyzer[seqlist,lenseqlist];
   lengthpar = Length[parameters];
   Print["The system for the coefficients c[i] has ",
          lengthpar," parameter(s)."];
   complexanalysis=lengthpar+Length[maineqlist]+Length[unknownlist];
   If[complexanalysis > 20, 
   coefmatrix = coefmat[system,unknownlist];
   Print["The system of PDEs, rank and form of rho, the system for"];
   Print["c[i], its coefficient matrix, and the lists of unknowns and"]; 
   Print["parameters, will all be saved in the file worklog.m for"];
   Print["further analysis. For analysis use the Mathematica functions"]; 
   Print["Reduce, Solve, etc. Load the file worklog.m with the command"]; 
   Print["<<worklog.m ."];
   Save["worklog.m",myfile,eqlist,rhorank,formrho,system,coefmatrix,
   unknownlist,parameters];
   Clear[eqlist,system,coefmatrix];
   ];
 For[i=1,i <= lengthpar,i++,
    syscondlist1 = Append[syscondlist1,parameters[[i]] != 0];
    ];
   syscond = Apply[And,syscondlist1];
   Print[" "];
   Print["Starting the search for compatibility conditions."];
   Print[" "];
   Print["The program will try all possible choices for non-vanishing"]; 
   Print["c[i]."];
   Print[" "];
   inputlist = Complement[unknownlist,analyzelist];
   If[inputlist === {},
     Print["All of the coefficients c[i] in the density have to vanish."];
     Print["The search for compatibility conditions ends."]; 
     ];
   While[inputlist =!= {},
    inputpart = Part[inputlist,1];
    inputval = inputpart == 1;
    inputlist = Complement[inputlist,{inputpart}];
    inputrule = ToRules[inputval];
    Print[" "];
    Print["-------------------------------------------------------"];
    Print[" "];
    Print["* Setting ",inputrule[[1]],":"];
    Print[" "];
    Print["Computation of the compatibility conditions."];
    comcond = Eliminate[maineqlist /. inputrule,unknownlist];
    Print[" "];
    comcondfac = MapAll[Factor,comcond];
    (* Remove possible duplicates *)
    If[Head[comcondfac] === Or, 
      comcondfactab = Table[Part[comcondfac,k], {k,1,Length[comcondfac]}];
      comcondfac = Union[comcondfactab]
      ];
    Print["This is the compatibility condition:"];
    Print[" "];
    Print[comcondfac];
    sol1 = Reduce[comcond && syscond,parameters];
    If[comcondfac === True,   
      Print[" "];
      Print["The compatibility condition is satisfied without constraints"];
      Print["on the parameters."] 
      ];
    If[sol1 =!= False,
      If[!FreeQ[sol1,Or],
        lengthsol1 = Length[sol1];
        For[iii = 1,iii <= lengthsol1,iii++,
           nrules1 = ToRules[Part[sol1,iii]];
           If[comcondfac=!=True, 
             For[k=1,k <= lengthpar,k++,
                par[k]=parameters[[k]] /. nrules1;   
                Print[" "];
                Print["For ",parameters[[k]]," = ",par[k] ];
                ] 
             ];
           newmaineqlist = maineqlist /.nrules1;
           solc = Flatten[Solve[newmaineqlist /. inputrule,unknownlist]];
           solc = Union[solc,inputrule];
           Print[" "];
           Print["The solution of the system is:"];
           Print[" "];
           Print[solc];
           nrules2 = rules /. nrules1;
           evaluate[formrho[x,t],rhot[x,t],solc,nrules2];
           inputlist = Intersection[inputlist,
             Union[solc[[Flatten[Position[solc,_->0]]]] /. x:(c[a_]->0)->c[a],
                   Complement[inputlist,First[#]& /@ solc]]];
           Clear[nrules1,nrules2,solc,newmaineqlist];
           ] (* for *),
        nrules1 = ToRules[sol1];
        If[comcondfac=!=True, 
           For[k = 1,k <= lengthpar,k++,
              par[k]=parameters[[k]] /. nrules1;
              Print[" "];
              Print["For ",parameters[[k]]," = ", par[k]];
              ]
           ];
        newmaineqlist = maineqlist /. nrules1;
        solc = Flatten[Solve[newmaineqlist /. inputrule,unknownlist]];
        solc = Union[solc,inputrule];
        Print[" "];
        Print["The solution of the system is "];
        Print[" "];
        Print[solc];
        nrules2 = rules /. nrules1;
        evaluate[formrho[x,t],rhot[x,t],solc,nrules2];
        inputlist = Intersection[inputlist,
             Union[solc[[Flatten[Position[solc,_->0]]]] /. x:(c[a_]->0)->c[a],
                   Complement[inputlist,First[#]& /@ solc]]];
        Clear[nrules1,nrules2,solc,newmaineqlist];
        ] (* if *),
      Print[" "];
      Print["The system becomes inconsistent, or the compatibility conditions"];
      Print["require that one or more of the parameters are zero."];
      Print["Not acceptable!"];
      ] ] (* closes While *),   
   nsolc = Flatten[Solve[maineqlist,unknownlist]];
   Print[" "];
   Print["Solution of the system: "];
   Print[" "];
   Print[nsolc];
   evaluate[formrho[x,t],rhot[x,t],nsolc,rules];
   Clear[solc1,solc2,solc3,nsolc];
   ];
  ]; (* end point 1 *) ]; (* end main *)
 
(* ------------------------------------------------------------------------- *)

(* set-up for collecting the data in a log file, transcript of computations *)
logfile="";
OpenLog[filename_String] := (logfile = OpenWrite[filename];
  If[logfile === $Failed, Return[]];
  AppendTo[$Echo,logfile]; AppendTo[$Output,logfile];);
CloseLog[] := ($Echo = Complement[$Echo,{logfile}];
  $Output = Complement[$Output,{logfile}];
  Close[logfile];);

(* End of the auxiliary functions and procedures *)
(* Start of the executable code lines *)

menu;
OpenLog[myfile];
commentinter[]; 
If[formrho[x,t] === 0, formrho[x,t] = {}];
If[Not[ListQ[formrho[x,t]]], formrho[x,t] = {formrho[x,t]}];

Print[" "];
Print["Working with the data file for the ",name,"."]; 
Print[" "];

For[i = 1,i <= noeqs,i++,
   eq[i][x,t] = Expand[eq[i][x,t]];
   Print[" "];
   Print["Equation ",i," of the system with ",noeqs," equation(s):"];
   Print[" "];
   Print["  ",pdeform[eq[i][x,t]]," = 0"];
   ];

eq[x,t] = Sum[eq[i][x,t],{i,1,noeqs}];
highestordereq = highestorder[eq[x,t],False,False,False,False,True,False];
Print[" "];
Print["Highest order of derivative in the system is ",highestordereq,"."];

If[formrho[x,t] === {}, (* start point rhogiven *)
  eqlist = Flatten[Table[eq[i][x,t],{i,1,noeqs}]];
  scalerules = scaling[eqlist];
  counternegweight = 0;
  listfreeweights = {};
  varscalelist = {};
  allchoices = {};
  Print[" "];
  Print["Program determines the weights of the variables (and parameters)."];
  For[i = 1,i <= noeqs,i++,
     weightu[i] = weightu[i] /. scalerules;
     If[i === 1, Print[" "]; Print["For the given system:"]; Print[" "]];
     Print["* weight of u[",i,"] is ",weightu[i],"."]; 
     varscalelist = Union[varscalelist,{{u[i][x,t],weightu[i]}}];
     If[Not[NumberQ[weightu[i]]],
       If[Length[weightu[i]] === 1,
         listfreeweights = Union[listfreeweights,{weightu[i]}]],
       If[weightu[i] < 0,counternegweight++];
       ];
     ];
  lenweightpars = Length[weightpars];
  For[i = 1,i <= lenweightpars,i++,
     weight[Part[weightpars,i]] = 
             weight[Part[weightpars,i]] /. scalerules; 
     Print["* weight of ",Part[weightpars,i]," is ",
           weight[Part[weightpars,i]],"."];
     varscalelist = Union[varscalelist,
                    {{Part[weightpars,i],weight[Part[weightpars,i]]}}];
     If[Not[NumberQ[weight[Part[weightpars,i]]]],
       If[Length[weight[Part[weightpars,i]]] === 1,
         listfreeweights = 
            Union[listfreeweights,{weight[Part[weightpars,i]]}]],
       If[weight[Part[weightpars,i]] < 0,counternegweight++];
       ];
     ]; 
  weight[d/dt] = weight[d/dt] /. scalerules;
  Print["* weight of d/dt is ",weight[d/dt],"."];
  tvarscalelist = Union[varscalelist,{{"d/dt",weight[d/dt]}}];
  If[Not[NumberQ[weight[d/dt]]],
    If[Length[weight[d/dt]] === 1,
      listfreeweights = Union[listfreeweights,{weight[d/dt]}]],
    If[weight[d/dt] < 0,
      Print[" "]; 
      Print["CAUTION! Weight of d/dt is negative. Computations"];
      Print["will continue, but the program is not designed for this"]'
      Print["case. Results may not be completely reliable."];
      ];
    ];
  counterfreeweight = Length[listfreeweights];
  If[counternegweight =!= 0,        (* start of If 1 *)
    Print[" "];
    Print["One or more of the weights are negative."];
    Print["Negative weights are not allowed."];
    Print["Aborting the computations!"];
    CloseLog[]; Abort[], 
    If[counterfreeweight =!= 0,
      If[counterfreeweight > 1,
        Print[" "];
        Print["Two or more of weights have freedom."];
        Print["Enter your values for the free weights or type Abort[];"];
        Print["Your values for the free weights can be entered by typing"];
        Print["`weightu[label] = val' or `weight[variable] = val' and"];
        Print["putting `;' between your entries."];
        Input[": "], (* else of counterfreeweight > 1 *)
        Print[" "];
        Print["Program will try to determine CHOICES for ",
              listfreeweights[[1]],":"];
        Print[" "];
        tempvarscalelist = Sort[Table[Part[tvarscalelist,k][[2]], 
                                 {k,1,Length[tvarscalelist]}]];
        tempvarscalelist = 
               Complement[tempvarscalelist,Select[tempvarscalelist,NumberQ]];
        lentempvarscalelist = Length[tempvarscalelist];
        Do[ (* begin Do loop *)
          Print["* CHOICE ",k];
          Print[" "];
          Print["Solving the equation: ",Part[tempvarscalelist,k]," = 1."];
          Print[" "];
          attempt = Flatten[Solve[Part[tempvarscalelist,k] == 1]];
          If[attempt =!= {},
            solfreeweight = Part[listfreeweights,1] /. attempt;
            While[solfreeweight <= 0, 
                 solfreeweight = solfreeweight + 1;
                 Print["Since the weight was zero or negative,"];
                 Print["it was incremented with 1."];
                 Print[" "]
                 ];
            Print[Part[listfreeweights,1]," = ",solfreeweight];
            weightrule = Part[listfreeweights,1] -> solfreeweight;
            tempvarscalelistval = tempvarscalelist /. weightrule;
            If[MemberQ[NonNegative[tempvarscalelistval],False],
              Print["Choice is rejected!"],
              varscalelistval = varscalelist /. weightrule;
              tvarscalelistval = tvarscalelist /. weightrule;    
              Print[" "];
              Print["List of all the variables (and parameters) with their"]; 
              Print["weights:"];
              Print[" "];
              Print[tvarscalelistval];
              Print[" "];
              testvarscalelist = Table[Part[varscalelistval,k][[2]],
                        {k,1,Length[varscalelist]}];
              If[Union[Positive[testvarscalelist]] === {False},
                Print["Since all the weights of u[i] are zero"];
                Print["this choice is rejected!"];
                Print[" "],
                allchoices = Union[allchoices,{solfreeweight}];
                ];
              ];
            ],{k,1,lentempvarscalelist} ];  (* end do *)
        Print[" "];
        Print["Simple POSITIVE choices are considered."];
        Print[" "];
        intchoices = Select[allchoices,IntegerQ];
        fracchoices = Complement[allchoices,intchoices];
        If[intchoices =!= {},
          (* Picking the minimum integer choice *)
          intchoices = Min[intchoices]; 
          choicerule = Part[listfreeweights,1] -> intchoices, (* else *)
          If[fracchoices =!= {},
           (* Picking the minimum fractional choice *)
           fracchoices = Min[fracchoices]; 
           choicerule = Part[listfreeweights,1] -> fracchoices, (* else *)
           Print["Enter your choice for the value of ",Part[listfreeweights,1],
                 "by typing its value. NO semi-colon at the end!"];
           choicerule = Part[listfreeweights,1]-> Input[": "];
           ]; 
          ]; 
        varscalelistchoice = varscalelist /. choicerule;
        tvarscalelistchoice = tvarscalelist /. choicerule;
        If[Length[allchoices] > 1,
          Print["List of all positive choices considered: "];
          Print[" "];
          Print[Union[Table[Part[listfreeweights,1] -> Part[allchoices,k], 
                             {k,1,Length[allchoices]}]]];
          Print[" "];
          Print["Some of the weights of u[i], but not all, could be zero."];
          Print["In the data file you could enter your choices in the format"];
          Print["`weightu[label] = value;' or `weight[variable] = value;'."];
          Print[" "]; 
          Print["Program continues with the choice: ", choicerule];
          Print[" "];
          Print["corresponding to the weights:"];
          Print[" "]; 
          Print[tvarscalelistchoice]; 
          Print[" "]; 
          ];
        varscalelist = varscalelistchoice;
      ] ] ]; (* end if 1 *)

  Print[" "];
  Print["The rank of rho should be an integer multiple of the lowest weight"];
  Print["of the DEPENDENT variable(s). Fractional weights are allowed."];
  Print[" "];
  If[
    formrho[x,t] === {},
    If[Not[NumberQ[rhorank]], rhorank = Input["Enter the rank of rho: "]];
    Print[" "]
    ];
  If[formrho[x,t] === {},
    Print["Computation of the density (and flux) of RANK = ",rhorank]; 
    Print["  "];
    ];  

  nodims = noeqs + Length[weightpars];
  varscalelist = Reverse[Sort[varscalelist,OrderedQ[{#1[[2]],#2[[2]]}]&]];
  constructformrho[nodims,varscalelist], (* else point for rhogiven *)
  formrho[x,t] = Part[formrho[x,t],1];
  formrho[x,t] = Expand[formrho[x,t]];
  If[Head[formrho[x,t]] === Times, 
    lenformrho = 1,lenformrho = Length[formrho[x,t]];
    ];
  Print[" "];
  Print["The program will only test the given form of the density rho."];
  Print["No determination of scaling properties."];
  Print["You have given this form for rho :"];
  Print[" "]; 
  Print[pdeform[formrho[x,t]]];
  Print[" "];
  ]; (* end point rhogiven *)

highestorderrho=highestorder[formrho[x,t],False,False,False,False,True,False];
Print["The highest-order of the terms in rho is ",
highestorderrho,"."];
eqlist = Flatten[Table[eq[i][x,t],{i,1,noeqs}]];
unknownlist = Table[c[i],{i,1,lenformrho}];
parameters = Union[parameters,weightpars];
weightpars = {};
 
Clear[i,k,lenformrho,counternegweight,counterfreeweight,scalerules,
     lenweightpars,varscalelist,tempvarscalelist,tvarscalelist,
     attempt,listfreeweights,allchoices,varscalelistval,temparscalelistval,
     tvarscalelistval,intchoices,fracchoices,choicerule,weightrule,
     nodims,highestordereq,weightu,lentempvarscalelist,solfreeweight,
     varscalelistchoice,tvarscalelistchoice,testvarscalelist,name,weight];

(* Solving for u[i]_t and calculation of the x-derivatives *)
mainsolforlist = Table[D[u[i][x,t],t],{i,1,noeqs}]; 
If[debug,
  Print["*************************************************"];
  Print["Trying to solve the equation(s) for ",mainsolforlist,","];
  Print["respectively."]; 
  Print[" "];
  Print["*************************************************"] ];
utlist = {};
For[i = 1,i <= noeqs,i++,
   ut[i] = Solve[eq[i][x,t] == 0,mainsolforlist[[i]]][[1]];
   utlist = Flatten[Append[utlist,ut[i]]];
   ];

Clear[ut,mainsolforlist,i];

main[utlist];

Print[" "];
Print["*********************** SUMMARY ***********************"];
Print[" "];              
Print["Total CPU time used in the current session is ",Round[TimeUsed[]],
      " seconds."];
Print[" "];
Print["To see the density type: rho[x,t], or use pdeform[rho[x,t]]."];
Print[" "];
Print["To see the flux type: j[x,t], or pdeform[j[x,t]]."];
Print[" "];
Print["To split the density in independent pieces type: split[rho[x,t]],"];
Print["or split[pdeform[rho[x,t]]]."];

CloseLog[];

(* ******************************* END ************************************* *)
