###############################################################################
#
#  The procedure varder calculates the variational derivative of
#  a functional L with respect to u.
#
###############################################################################

varder:=proc(L,u,x)
option `Authors: B. Deconinck`;
local N,k,lag,v,diffvec,m;
N:=hdifforder(L,u,x);
# N is the order of the highest derivative of u w.r.t. x occuring in L
lag:=L;
for k from N by -1 to 1 do
 lag:=subs(diff(u,[seq(x,j=1..k)])=v[k],lag)
od;
lag:=subs(u=v[0],lag);
diffvec:=[seq(diff(lag,v[k]),k=0..N)];
diffvec:=map2(subs,v[0]=u,diffvec);
diffvec:=map2(subs,[seq(v[k]=diff(u,[seq(x,m=1..k)]),k=1..N)],diffvec);
# The last result serves as the output of the procedure
(expand(diffvec[1]+add((-1)^k*diff(diffvec[k+1],[seq(x,m=1..k)]),k=1..N)))
end:

###############################################################################
#
# The procedure hdifforder finds the order of the highest
# derivative of u w.r.t. x occuring in L.
#
###############################################################################

hdifforder:=proc(L,u,x)
option `Authors: B. Deconinck`;
local k,lijst,flag;
k:=0; flag:=1; lijst:=[];
while flag=1 do
 if has(L,diff(u,lijst)) then
  lijst:=[op(lijst),x];
  k:=k+1
 else flag:=0
 fi
od;
# the last result serves as the output of the procedure
if k=0 then 0
else k-1
fi;
end:

###############################################################################
#
# The procedure varder calculates the higher-order
# variational derivative of a functional L with respect to u.
# This procedure cannot be used to compute the variational
# derivative. We require n to be at least 1.
#
###############################################################################

vardern:=proc(L,u,n,x)
option `Authors: B. Deconinck`;
local N,k,lag,v,diffvec,m;
N:=hdifforder(L,u,x);
if N<n then return(0)
fi;
# N is the order of the highest derivative of u w.r.t. x occuring in L
lag:=L;
for k from N by -1 to n do
 lag:=subs(diff(u,[seq(x,j=1..k)])=v[k],lag)
od;
diffvec:=[seq(diff(lag,v[k]),k=n..N)];
diffvec:=map2(subs,[seq(v[k]=diff(u,[seq(x,m=1..k)]),k=n..N)],diffvec);
# The last result serves as the output of the procedure
expand(diffvec[1]+add((-1)^(k-n)*binomial(k,n)*diff(diffvec[k-n+1],[seq(x,m=n..k-1)]),k=n+1..N))
end:

#####################################################################################################
#
# The homotopy procedure, to compute integrals of functions that
# contain unknown functions. Here u is a vector of unknown functions
# that depend on an independent variable x.
#
#####################################################################################################

homotopyintegral:=proc(LL,u,x,a0)
    option `Authors: B. Deconinck`;
    local aa,limset,ans,i,r,m,lambda,lambda0,flux,scaledflux,summedflux,N,n,L,h,a,s;
    L:=expand(eval(LL));
    n:=nops(u); # number of unknown functions
    # N is the vector of differential orders of u in L
    N:=[seq(hdifforder(L,u[r],x),r=1..n)];
    # flux is the integrand before scaling of the arguments
    flux:=expand([seq(u[r]*vardern(L,u[r],1,x)+add(diff(u[r]*vardern(L,u[r],i+1,x),[seq(x,m=1..i)]),i=1..N[r]-1),r=1..n)]);
    summedflux:=expand(add(flux[r],r=1..n));
    scaledflux:=expand(subs({seq(u[r]=lambda*u[r],r=1..n)},summedflux));
    h:=r->eval(expand(int(expand(scaledflux/lambda),lambda)),lambda=r);#print(h(a));
    ans:=expand(expand(h(1)-limit(h(a),a=a0)));
end:

################################################################################################
#
# This is a modified homotopy integral, a list of terms is returned
# Here u is a vector of unknown functions
# that depend on an independent variable x.
#
################################################################################################

vechomotopyintegral:=proc(LL,u,x)
option `Authors: B. Deconinck, M. Nivala`;
local sh,result,newans,newflux,aa,limset,ans,i,r,m,lambda,lambda0,flux,scaledflux,summedflux,N,n,L,h,a,s;
L:=expand(eval(LL));
n:=nops(u); # number of unknown functions
# N is the vector of differential orders of u in L
N:=[seq(hdifforder(L,u[r],x),r=1..n)];
# flux is the integrand before scaling of the arguments
# to obtain the maximum number of terms, we take this as a list and a sum
flux:=expand([seq(op([u[r]*vardern(L,u[r],1,x),seq(diff(u[r]*vardern(L,u[r],i+1,x),[seq(x,m=1..i)]),i=1..N[r]-1)]),r=1..n)]);
scaledflux:=expand(subs({seq(u[r]=lambda*u[r],r=1..n)},flux));
sh:=normal(convert(scaledflux,`+`));
ans:=expand({seq(int(expand(scaledflux[i]/lambda),lambda),i=1..nops(scaledflux)),int(expand(sh/lambda),lambda)});
#remove expressions that were not integrated
for s in ans do
 if has(s,int) then
  ans:=ans minus {s}
 fi;
od;
#get as many terms as possible
newans:=[NULL];
for i from 1 to nops(ans) do
 if whattype(ans[i])=`+` then
  newans:=[op(newans),op(ans[i])];
 else
  newans:=[op(newans),ans[i]];
 fi;
od;
#normalize the solution and remove constants
result:={NULL};
for i from 1 to nops(newans) do
 if has(newans[i],x) then
  result:={op(result),SPLIT(newans[i],x)[1]};
 fi;
od;
[op(subs(lambda=1,result))];
end:

###############################################################################
#
# whatfun picks out undefined functions in an expression
#
###############################################################################

whatfun:=proc(E)
option `Authors: M. Nivala`;
local n,s,sid;
s:={NULL};
sid:=[op(indets(E,`function`))];
for n from 1 to nops(sid) do
 if not(has(sid[n],`diff`)) and not(type(op(0,sid[n]),`procedure`))
  then s:={op(s),sid[n]};
 fi;
od;
[op(s)];
end:

##############################################################################
#
# orderset orders a set from large to small
#
##############################################################################

orderset:=proc(s)
option `Authors: B. Deconinck, M. Nivala`;
    local ss,leftover,k;
    ss:=NULL;
    leftover:=s;
    for k to nops(s) do
        ss:=ss,max(op(leftover));
        leftover:=leftover minus {max(op(leftover))}
    od;
    [ss]
end:

#################################################################################
#
# funlexorderconverts a list of functional expressions to strings, orders them
# according to lexorder, and returns the funtional expressions
#
#################################################################################

funlexorder:=proc(LIST)
option `Authors: M. Nivala`;
local Newlist,i,subrules,s;
subrules:=[seq(convert(LIST[i],`string`)=LIST[i],i=1..nops(LIST))];
Newlist:=[seq(op(1,subrules[i]),i=1..nops(LIST))];
s:=sort(Newlist,`lexorder`);
subs(subrules,s);
end:

#######################################################################################
#
# orders a set according to order of derivatives, from highest to lowest,
# grouping those of the same order
#
#######################################################################################

dorderlist:=proc(L,u,x)
option `Authors: M. Nivala`;
local finlist,j,diffset,dlist,i,MAX,difflist,ord;
difflist:=[seq(hdifforder(L[i],u,x),i=1..nops(L))];
diffset:=orderset({op(difflist)});
finlist:=[NULL];
for i in diffset do
 dlist:=[NULL];
 for j from 1 to nops(L) do
  if i=difflist[j] then
   dlist:=[op(dlist),L[j]];
  fi;
 od;
finlist:=[op(finlist),dlist];
od;
finlist;
end:

###################################################################################
#
# torderset orders a set according to number of products in a term
#
###################################################################################

torderset:=proc(E)
option `Authors: M. Nivala`;
local w,fin,n,s;
s:={seq(nops(op(1,E[n]))+(n/(nops(E)+1)),n=1..nops(E))};
s:=orderset(s);
w:={seq(nops(op(1,E[n]))+(n/(nops(E)+1))=E[n],n=1..nops(E))};
fin:=subs(w,s);
end:

###################################################################################
#
# SPLIT splits the product T into an x-dependent part and a constant part
#
###################################################################################

SPLIT:=proc(T,x)
option `Authors: M. Nivala`;
local xpart,coeffpart,s,opT,xset,coeffset;
xset:={NULL};
coeffset:={NULL};
if whattype(T)=`*` then
 opT:=[op(T)];
else
 opT:=[T];
fi;
for s in opT do
 if has(s,x) then
  xset:={op(xset),s};
 else
  coeffset:={op(coeffset),s};
 fi;
od;
xpart:=convert(xset,`*`);
coeffpart:=convert(coeffset,`*`);
return([xpart,coeffpart]);
end:

#################################################################################
#
# SUBS substitutes ssubrules into EE
#
#################################################################################

SUBS:=proc(ssubrules,EE,x)
option `Authors: M. Nivala`;
local E,Elist,subrules;
subrules:=torderset(ssubrules);
E:=expand(EE);
if not(whattype(E)=`+`) then
 Elist:=[E];
else
 Elist:=convert(E,`list`);
fi;
Elist:=[seq(SPLIT(Elist[i],x),i=1..nops(Elist))];
Elist:=subs(subrules,Elist);
E:=add(convert(Elist[i],`*`),i=1..nops(Elist));
end:

##################################################################################
#
# GROUP goups terms with the same x-dependence in E
#
##################################################################################

GROUP:=proc(E,x)
option `Authors: M. Nivala`;
local EL,subrules,EE,Esubs,newsubrules,i,b;
if not(whattype(E)=`+`) then
  EL:=[E];
else
 EL:=convert(E,`list`);
fi;
subrules:=[op({seq(SPLIT(EL[i],x)[1]=b[i],i=1..nops(EL))})];
subrules:=torderset(subrules);
Esubs:=SUBS(subrules,E,x);
Esubs:=collect(Esubs,{seq(b[i],i=1..nops(subrules))});
newsubrules:=[seq(op(2,subrules[i])=op(1,subrules[i]),i=1..nops(subrules))];
EE:=subs(newsubrules,Esubs);
end:

################################################################################################
#
# RRED performs gaussian elimination with a minimal number of row switching operations
# it row-reduces A while preserving the order of the rows
# returns the first maxsstep-1 linearly independent rows
#
################################################################################################

RRED:=proc(A,sstep,B,maxsstep)
option `Authors: M. Nivala`;
local step,i,j,newA,flag1,rowdims,newB;
rowdims:=linalg[rowdim](A);
if sstep=maxsstep then
 return(B);
fi;
step:=sstep;
newA:=A;
flag1:=0;
for i from 1 to rowdims do
  if newA[i,step]<>0 and flag1=0 then
   flag1:=i;
  fi;
od;
if flag1<>0 then
 newA:=linalg[pivot](newA,flag1,step);
 if step=1 then
  newB:=linalg[matrix]([linalg[row](newA,flag1)]);
 else
  newB:=linalg[stackmatrix](B,linalg[row](newA,flag1));
 fi;
 if rowdims=1 then
  return(newB);
 else
  newA:=linalg[delrows](newA,flag1..flag1);
  RRED(newA,step+1,newB,maxsstep);
 fi;
else
 RRED(newA,step+1,newB,maxsstep);
fi;
end:

#################################################################################################################################
#
# INT is used to compute integrals of expressions that depend on unknown functions that depend on an independent variable x.
# All terms that are total derivatives are integrated.
# The order of the highest derivative in the remaining integral is minimized.
#
#################################################################################################################################

INT:=proc(EE,x)
option `Authors: M. Nivala`;
local dcheck,limset,eqnset2,RE,param, eqnset,v,slnset,novan,Flist,intnogood,nogood,subDFL,intxset,sln,Esub,Fsub,DFLsub,termlist,FL,xset,E,u,goo,A,Ar,Ag,bob,syst,i,j,m,DFL,s,n,fresult,result,a,E1,DF1,R,subrules,b,EL,F1,F,DF,k;

E:=expand(eval(EE));

#determine the unknown functions
u:=whatfun(E);
if u=[NULL] then return(int(E,x))
fi;

#check if E is exact
dcheck:=[seq(varder(E,u[i],x),i=1..nops(u))];
if dcheck=[seq(0,i=1..nops(u))] then
 return(expand(int(E,x)));
fi;

#convert E to a list
if not(whattype(E)=`+`) then
 EL:=[E];
else
 EL:=convert(E,`list`);
fi;

#remove terms which have no functional or x dependence to be integrated directly
xset:={NULL};
for k in EL do
if not(has(k,u)) or not(has(k,x)) then
 xset:={op(xset),k};
fi;
od;
xset:=expand(convert(xset,`+`));
intxset:=int(xset,x);
E:=expand(eval(E-xset));
if not(whattype(E)=`+`) then
 EL:=[E];
else
 EL:=convert(E,`list`);
fi;

#FL is a list of terms from the homotopy integral of E
FL:=vechomotopyintegral(E,u,x); #print(FFFFFFF,FL);
if FL=[NULL] then
 return(int(E,x)+intxset);
fi;

#DFL is the derivative of FL with respect to x
DFL:=[seq(expand(diff(FL[n],x)),n=1..nops(FL))];

#form substitution rules
#subDFL is for substitution rules
subDFL:={NULL};
for k in DFL do
 if not(whattype(k)=`+`) then
   subDFL:={op(subDFL),k};
 else
  subDFL:={op(subDFL),op(k)};
 fi;
od;
termlist:=[op({op(EL),op(FL),op(subDFL)})];
#remove constant coefficients
termlist:=[op({seq(SPLIT(termlist[i],x)[1],i=1..nops(termlist))})];
#orders termlist according to the order of the highest derivative of each element, from highest to lowest, groouping those of the
#same order
termlist:=dorderlist(termlist,u,x);#print(dord,termlist);
#for consistency, implement an ordering of funtions and their derivatives according to lexorder for those terms which are of the
#same differential order
termlist:=[seq(op(funlexorder(termlist[i])),i=1..nops(termlist))];#print(orderedsub,termlist);

#perform the substitutions
subrules:=(seq(termlist[n]=b[n],n=1..nops(termlist)));
#order subrules according to number of products in each term
subrules:=torderset({subrules});#print(SUBRULES,subrules);
Esub:=SUBS(subrules,E,x);
DFLsub:=[seq(SUBS(subrules,DFL[i],x),i=1..nops(DFL))];

#R is what remains after integrating E
#R is optimized, removing higest order derivaites first
R:=expand(Esub-sum(a[m]*DFLsub[m],m=1..nops(DFL)));#print(E-sum(a[m]*DFL[m],m=1..nops(DFL)));
R:=collect(R,{seq(b[m],m=1..nops(termlist))});#print(RRR,R);
#form the linear system to be solved
syst:={NULL};
for m from 1 to nops(termlist) do
 if has(coeff(R,b[m]),[seq(a[n],n=1..nops(DFL))]) then
  syst:=[op(syst),coeff(R,b[m])=0];
 fi;
od;

#A is the overdetermined system for the coefficients, a
#Ar is the first nops(a) linearly independent rows of A
A:=linalg[genmatrix](syst,[seq(a[j],j=1..nops(DFLsub))],flag);#print(AAAA,A);
Ar:=RRED(A,1,A,nops(DFLsub)+1);#print(ArAr,Ar);

slnset:=linalg[backsub](Ar,`false`,`v`);#print(slnset);
sln:=seq(a[m]=slnset[m],m=1..nops(DFLsub));#print(solution,sln);
sln:=eval(sln,v=0);#print(solution,sln);
assign(sln);

#all terms which were not integrated using homotopy, attempt to integrate each independently without homotopy
nogood:=expand(E-add(a[i]*DFL[i],i=1..nops(DFL)));#print(nogood);
if not(whattype(nogood)=`+`) then
  nogood:=[nogood];
else
 nogood:=convert(nogood,`list`);
fi;
intnogood:=expand(add(int(nogood[i],x),i=1..nops(nogood)));#print(intnogood);
#print(expand(E-add(a[i]*DFL[i],i=1..nops(DFL))+xset+novan+RE),"Was Integrated By Maple");

fresult:=add(a[j]*FL[j],j=1..nops(FL))+combine(intnogood)+intxset;#print(fdfdfdfd,fresult);
unassign(a);

#check if correct
if simplify(expand(EE-diff(fresult,x)))<>0 then
 print("No integration done.");
 return(int(E,x)+intxset)
else
 return(fresult);
fi;
end:

########################################################################################################
#
# ouputs the expression G modulo total derivatives
#
########################################################################################################

DiffMod:=proc(G,x)
local IG,n,intpart,Glist;
IG:=INT(G,x);
if not(whattype(IG)=`+`) then
 Glist:=[IG];
else
 Glist:=convert(IG,`list`);
fi;
intpart:=0;
for n from 1 to nops(Glist) do
 if has(Glist[n],`int`) then intpart:=op(1,Glist[n]);
 fi;
od;
intpart;
end:

#########################################
#
# checks if INT is correct
#
#########################################

INTcheck:=proc(E,F,x)
normal(expand(E-diff(F,x)));
end:

#####################################################################################################
#
# the procedure weight calculates the weight of every
# unknown dependent and independent variable, choosing
# the weight of the derivate w.r.t x to be 1. Here L is the
# differential expression, w is a list of all the functions and
# variables appearing in L, with the independent variables
# usually given last. Then "x" denotes which one of these should
# have weight 1.
#
######################################################################################################

weight:=proc(LL,w,x)
local ans,v,s,i,A,Ar,sln,slnset,L,j,GM,m,M,ters,syscoef,firstrules,firstbob,uu,tt,N,M1,M2,k,n,t,rules,bob,sys,newbob,result;
L:=expand(LL);
N:=nops(L); # the number of equations and dependent variables
t:=op({op(w[1])} minus {x}); # the other independent variable
M1:=PDEtools[difforder](L,x);
M2:=PDEtools[difforder](L,t);
# rules is a list of replacement for w and its derivatives in terms
# of powers of a scaling variable lambda
firstrules:=[seq(conjugate(w[k])=w[k],k=1..nops(w))];
firstbob:=subs(firstrules,L);
rules:=[seq(seq(diff(w[k],t$n)=lambda^(uu[k]+n*tt),n=1..M2),k=1..nops(w)),seq(seq(diff(w[k],x$n)=lambda^(uu[k]+n*1),n=1..M1),k=1..nops(w)),seq(w[k]=lambda^uu[k],k=1..nops(w))];
bob:=subs(rules,firstbob);
# newbob is normalizing the coefficient of the first term to be 1
for k from 1 to N do
 newbob[k]:=expand(bob[k]/(op(1,bob[k])))
od ;

# Using the derivative as a way of getting the exponent
sys:=simplify([seq(seq(lambda*diff(op(n,newbob[k]),lambda)/op(n,newbob[k]),n=1..nops(newbob[k])),k=1..N)]);
ters:=[tt,seq(uu[k],k=1..nops(w))];

A:=linalg[genmatrix](sys,ters,flag);
Ar:=linalg[gausselim](A);#print(Ar);
slnset:=traperror(linalg[backsub](Ar,'false','v'));
if slnset="inconsistent system" then
 print("A scaling symetry does not exist, please add parameters.");
 return([NULL]);
fi;
sln:=[seq(slnset[jj],jj=2..nops(w)+1)];
#print(sadf,seq(plexsim(sln,w,N,i),i=0..N));
s:=simplex[minimize](add(sln[i],i=1..nops(sln)),{seq(sln[i]>=1,i=1..N),seq(sln[j]>=0,j=(N+1)..nops(w))});#print(sln,s,subs(s,sln),eval(subs(s,sln),v=0));
ans:=eval(subs(s,sln),v=0);
end:

#################################################################################################################################
#
# weightcheck checks for fractional weights, returning integer weights
#
#################################################################################################################################

weightcheck:=proc(w)
local wnew,s,maxd,neww,flag,fracset;
wnew:=w;
fracset:={NULL};
flag:=0;
for s in w do
 if type(s,fraction) then
  flag:=1;
 fi;
od;
if flag=1 then
 for s in w do
  fracset:={op(fracset),denom(s)};
 od;
 maxd:=lcm(op(fracset));
 neww:=[[seq(maxd*w[n],n=1..nops(w))],maxd];
else
 [wnew,1];
fi;
end:

##############################################################################################
#
# determines the weight of an expression given the weight (of the functions u) w.
#
##############################################################################################

weighteqn:=proc(E,u,w,x,N)
local newE,M1,M2,rules,lambda,k,n;
M1:=PDEtools[difforder](E,x);
rules:=[seq(seq(diff(u[k],x$n)=lambda^(w[k]+n*N),n=1..M1),k=1..nops(u)),seq(u[k]=lambda^w[k],k=1..nops(u))];
newE:=subs(rules,E);
degree(newE,lambda);
end:

#################################################################
#
# forms all monomials of weight W in u and its derivatives
#
#################################################################

endlesspossibilities:=proc(u,W,w,x,N)
local fset,tset,s,M,i;
M:=[seq(floor(W/w[i]),i=1..nops(u))];
tset:={op(expand(mul((add(u[i]^n,n=0..M[i])),i=1..nops(u))))} minus {1};
fset:={NULL};
for s in tset do
 if weighteqn(s,u,w,x,N)<=W then
  fset:={op(fset),s};
 fi;
od;
fset;
end:

#########################################################
#
#
#
#########################################################

rmvcon:=proc(slnset,coset,param)
local finset,i;
finset:={NULL};
for i in slnset do
 if not(has(i,coset)) and has(i,param) then
  finset:={op(finset),i}
 fi;
od;
finset;
end:

###########################################################
#
#
#
###########################################################

rmvineq:=proc(slnset)
local finset,i;
finset:={NULL};
for i in slnset do
 if whattype(i) = `=` then
  finset:={op(finset),i};
 fi;
od;
finset;
end:

############################################################
#
#
#
############################################################

paramfactor:=proc(result,param)
local n,s,fs,newresult;
newresult:={op(result)};
for s in result do
 fs:=factors(s)[2];
 for n from 1 to nops(fs) do
  if has(param,fs[n][1]) then newresult:=newresult minus {s}
  fi;
 od;
od;
[op(newresult)];
end:

##################################################################################
#
#################################################################################
degreduce:=proc(E,u,x,M)
local m,h,s,i,N,check;
N:=nops(u);
h:=[seq(hdifforder(E,u[i],x),i=1..N)];
m:=max(op(h));
if m=0 then
 return(0)
fi;
if numboccur(h,m)=1 then
 for i from 1 to N do
  if h[i]=m then
   check:=i;
  fi;
 od;
 if check=M then
  s:=degree(E,diff(u[check],x$m));
  if s>1 then
   return(0);
  else
   return(1);
  fi;
 else
  return(0);
 fi;
else
 return(0);
fi;
end:

##################################################################################
#
##################################################################################
degreducerun:=proc(S,u,x,M)
local i,SS,N,j,U,finS,finSS;
N:=nops(u);
U:=[seq(op(0,u[i]),i=1..nops(u))];
SS:=[seq({NULL},i=1..N),{NULL}];
for j from 1 to N do
 for i from 1 to nops(S) do
  if type(S[i],freeof({op(U)} minus {U[j]})) then
   SS[j]:={op(SS[j]),S[i]};
  fi;
 od;
od;
SS[N+1]:=S minus {seq(op(SS[i]),i=1..N)};
finS:={NULL};finSS:={NULL};
for j from 1 to N do
 for i from 1 to nops(SS[j]) do
  if degreduce(SS[j][i],u,x,j)=0 then
   finS:={op(finS),SS[j][i]};
  fi;
 od;
od;
for i from 1 to nops(SS[N+1]) do
 if degreduce(SS[N+1][i],u,x,M)=0 then
  finSS:={op(finSS),SS[N+1][i]};
 fi;
od;
[finS,finSS];
end:

##########################################################################################################################
#
# In conservation, kdv is a list of eqns, u is a list of functions, and W is the total weight
#
##########################################################################################################################

conservation:=proc(kdv,u,x,W)
local tempsln,bb,check,count,butsln,newsln,c,t,flagNx,Nd,Nx,j,aaa,ar,flagp,param,n,flagflux,Nu,Np,flagw,wt,possw,possu,paramw,i,a,aw,newa,newnewa,rho,rhot,ut,rhotsubs,vrho,rhotset,termlist,subrules,syst,s,slnsets,fluxslnsets,sln,fluxer,newfluxsln,fluxsln;

ar:=[args];
t:=op({op(u[1])} minus {x});

flagp:=0;

# check if parameters are given
for n from 1 to nops(ar) do
  if has(ar[n],parameters) then
   param:=op(2,ar[n]);
   flagp:=1;
  fi;
od;
if flagp=0 then
 param:=[NULL];
fi;

# check if flux is wanted
if has(ar,flux) then
 flagflux:=1
else
 flagflux:=0;
fi;

Nu:=nops(u);
Np:=nops(param);

flagw:=0;
flagNx:=0;
# check if weights are given (weights of dependent variables and parameters must be given), if not calculate them
for n from 1 to nops(ar) do
  if has(ar[n],weights) then
   flagw:=1;
   wt:=op(2,ar[n]);
  elif has(ar[n],weightx) then
    flagNx:=1;
    Nx:=op(2,ar[n]);
  fi;
od;

# if the weight of x is not given, then set it equal to 1
if flagNx=0 then
 Nx:=1;
fi;

# calculate weights and convert them to integer values.
if flagw=0 then
 wt:=weight(kdv,[op(u),op(param)],x);
 if wt=[NULL] then
  return(NULL);
 fi;
 wt:=weightcheck(wt);
 Nx:=wt[2];
 wt:=wt[1];
else
 wt:=weightcheck(wt);
 Nx:=wt[2]*Nx;
 wt:=wt[1];
fi;
if has(ar,returnweights) then
 print(weights=wt,Nx);
fi;

# check which parameters have non-zero weight
possw:=[NULL];
possu:=[NULL];
paramw:=[NULL];
for i from 1 to Nu do
 if wt[i]<>0 then
  possu:=[op(possu),u[i]];
  possw:=[op(possw),wt[i]];
 fi;
od;
for i from Nu+1 to Np+Nu do
 if wt[i]<>0 then
  possu:=[op(possu),param[i-Nu]];
  paramw:=[op(paramw),param[i-Nu]];
  possw:=[op(possw),wt[i]];
 fi;
od;

# form possible density
newa:={NULL};
a:=endlesspossibilities(possu,W,possw,x,Nx);
for s in a do
 if has(s,x) then
  newa:={op(newa),s}
 fi;
od;
a:=[op(newa)];#print(a,nops(a));
aw:=[seq(weighteqn(a[i],possu,possw,x,Nx),i=1..nops(a))];#print(aw);
a:=subs({seq(param[i]=1,i=1..nops(param))},a);#print(a,nops(a));
newa:={NULL};
for i from 1 to nops(a) do
 Nd:=(W-aw[i])/Nx;
 if Nd<>0 and evalb(whattype(Nd)=integer) then
  newa:={op(newa),diff(a[i],x$(Nd))};
 elif Nd=0 then
  newa:={op(newa),a[i]};
 fi;
od;
newa:=expand(newa);#print(newa);
newnewa:={NULL};
for i from 1 to nops(newa) do
 if not(whattype(newa[i])=`+`) then
  newnewa:={op(newnewa),newa[i]};
 else
  newnewa:={op(newnewa),op(newa[i])};
 fi;
od;
newnewa:=newnewa minus {0};

#print(nops(newnewa),startmod);

# reduce the number of terms that must go into DiffMod
newa:=degreducerun(newnewa,u,x,1);#print(newa);

#print(nops({op(newa[1]),op(newa[2])}),nops(newa[2]),enddegreduce);

newa[2]:=expand({seq(DiffMod(newa[2][i],x),i=1..nops(newa[2]))} minus {1});
newa:={op(newa[1]),op(newa[2])};
a:={NULL};
for i from 1 to nops(newa) do
 if not(whattype(newa[i])=`+`) then
  a:={op(a),newa[i]};
 else
  a:={op(a),op(newa[i])};
 fi;
od;
a:=newa;
a:={seq(SPLIT(a[i],x)[1],i=1..nops(a))} minus {1};

# for consistency, choose an ordering of possible terms
# orders a according to the order of the highest derivative of each element, from highest to lowest, groouping those of the
# same order
a:=dorderlist(a,u,x);#print(dord,a);
# for consistency, implement an ordering of functions and their derivatives according to lexorder for those terms which are of the
# same differential order
a:=[seq(op(funlexorder(a[i])),i=1..nops(a))];

#print(nops(a),enddiffmodstuff);

rho:=add(c[k]*a[k],k=1..nops(a));#print(rho);
rho:=subs({seq(param[i]=1,i=1..nops(param))},rho);
rhot:=diff(rho,t);
ut:={seq(diff(u[k],t)=-(kdv[k]-coeff(kdv[k],diff(u[k],t))*diff(u[k],t))/coeff(kdv[k],diff(u[k],t)),k=1..nops(u))};
rhotsubs:=expand(subs(ut,rhot));
vrho:=expand([seq(varder(rhotsubs,u[i],x),i=1..nops(u))]);

# form substitution rules
rhotset:={NULL};
for i from 1 to nops(vrho) do
 if not(whattype(vrho[i])=`+`) then
  rhotset:={op(rhotset),vrho[i]};
 else
  rhotset:={op(rhotset),op(vrho[i])};
 fi;
od;
termlist:=[op(rhotset minus {0})];
# remove constant coefficients
termlist:=[op({op({seq(SPLIT(termlist[i],x)[1],i=1..nops(termlist))})} minus {1})];
subrules:=(seq(termlist[n]=b[n],n=1..nops(termlist)));
subrules:=torderset({subrules});#print(subrules);

# set up the linear system
syst:=[seq(SUBS(subrules,vrho[i],x),i=1..nops(vrho))];
syst:=[seq(coeffs(syst[i],[seq(b[n],n=1..nops(termlist))]),i=1..nops(syst))];

# solve the system
if flagp=1 and {op(syst)}<>{0} then
 s:=[PDEtools[casesplit](syst,[seq(c[k],k=1..nops(a)),op(param)])];
 slnsets:=[seq({op(op(1,s[i])),op(op(2,s[i]))},i=1..nops(s))];
 slnsets:=[seq(solve(slnsets[i]),i=1..nops(slnsets))];
 sln:={seq([{coeffs(expand(subs(slnsets[i],rho)),[seq(c[k],k=1..nops(a))])} minus {0},rmvcon(slnsets[i],{seq(c[k],k=1..nops(a))},param)],i=1..nops(s))};
 newsln:=sln;
 for i from 1 to nops(sln) do
  if sln[i][1]={NULL} or sln[i][1]={0} then
   newsln:=newsln minus {sln[i]};
  fi;
 od;
 butsln:=[seq(op([seq([newsln[i][1][j],newsln[i][2]],j=1..nops(newsln[i][1]))]),i=1..nops(newsln))];
 newsln:=butsln;
 if flagflux=1 then
  for i from 1 to nops(newsln) do
   newsln[i]:=[newsln[i][1],homotopyintegral(expand(-eval(subs(subs(newsln[i][2],ut),diff(newsln[i][1],t)))),u,x,0),newsln[i][2]];
  od;
 fi;
else
 s:=op(solve(syst,[seq(c[k],k=1..nops(a))]));
 newsln:=[op({coeffs(expand(subs(s,rho)),[seq(c[k],k=1..nops(a))])} minus{0})];
 newsln:=[seq([newsln[i]],i=1..nops(newsln))];
 if flagflux=1 then
  newsln:=[seq([newsln[i][1],homotopyintegral(expand(-eval(subs(ut,diff(newsln[i][1],t)))),u,x,0)],i=1..nops(newsln))];
 fi;
fi;

# normalize solution
if newsln<>[NULL] then
 for i from 1 to nops(newsln) do
  check:=0;
  count:=1;
  while check=0 do
   if has(newsln[i][1],a[count]) then
    subrules:={(seq(a[n]=b[n],n=1..nops(a)))};
    subrules:=torderset(subrules);
    tempsln:=SUBS(subrules,newsln[i][1],x);
    if coeff(tempsln,b[count])<>0 then
     newsln[i][1]:=expand(newsln[i][1]/coeff(tempsln,b[count]));
     check:=1;
    else
     count:=count+1;
    fi;
   else
    count:=count+1;
   fi;
  od;
 od;
 return(op(newsln));
else
 return([NULL]);
fi;

end:

###################################################################
#
# ConservationLaw runs conservation for a list of weights, Wlist
#
###################################################################

ConservationLaw:=proc(kdv,u,x,Wlist)
local i,opts;
opts:=seq([args][i],i=5..nops([args]));
seq(conservation(kdv,u,x,Wlist[i],opts),i=1..nops(Wlist));
end:

#########################################################################################

modder1:=proc(E,u,x)
local subrules,numu,M,N,n,i,j,nx;
M:=nops(u);
N:=PDEtools[difforder](E,x);
subrules:=[seq([seq(diff(u[j],x$(N-n))=1,n=0..N-1)],j=1..M)];print(subrules);
numu:=[seq(0,i=1..M)];
for i from 1 to M do
 numu[i]:=[degree(subs(subrules[i],E),u[i]),seq(degree(subs([seq(subrules[i][j],j=1..(n-N))],E),diff(u[i],x$n)),n=1..N)];
 for j from 1 to nops(numu[i]) do
  if numu[i][j]=FAIL then
   numu[i][j]:=0;
  fi;
 od;
od;
nx:=0;
for i from 1 to M do
 nx:=nx+add((n-1)*numu[i][n],n=1..N+1);
od;
[[seq(convert(numu[n],`+`),n=1..M),nx],N];
end:

#######################################################################################

modder2:=proc(S,u,x)
local i,j,SS,s,count,flag;
s:=[seq(modder1(S[i],u,x),i=1..nops(S))];
SS:={op(S)};
for i from 1 to nops(s) do
 count:=1;
 flag:=0;
 while flag=0 and count<=nops(s) do
  if s[i][1]=s[count][1] and s[i][2]>s[count][2] then
   SS:=SS minus {S[i]};
   flag:=1;
  else
   count:=count+1;
  fi;
 od;
od;
SS;
end: