(* char.m - Updated: January 17, 1996 *)

characteristic[expr1_,expr2_] := Module[
{res,rhotjx,i,j,coef,charoper},

res = rhotjx = D[expr1,t]+D[expr2,x];

Print[" "];
Print["Characteristic operator(s) of this conservation law st."];
Print[" "];
Print["charoper[1] * Equation-1 +.....+ charoper[N] * Equation-N == 0"];
Print[" "];
Print["(where 1 <= N <= ",noeqs,") is/are :"];

For[i = 1,i <= noeqs,i++,
   charoper[i] = 0;
   For[j = 1,j <= highestorderrho,j++,
      coef = Coefficient[rhotjx,D[u[i][x,t],t,{x,j}]];
      charoper[i] = charoper[i]+coef*Dx^j;
      res = res - coef*D[eq[i][x,t],{x,j}];
      ];
   ];
For[i = 1,i <= noeqs,i++,
   coef = Coefficient[res,D[u[i][x,t],t]];
   charoper[i] = charoper[i] + coef;
   res = res - coef*eq[i][x,t];
   Print[" "];
   Print[" charoper[",i,"] = ",charoper[i]];
   ];
(* Print["res ",Simplify[res]]; *)
];

characteristic[rho[x,t],j[x,t]];