# Module syntax:
# [[x,y,...], [Dx,Dy,...], [the matrix of [Dx(b1),..,Dx(bn)], the matrix of [Dy(b1)..Dy(bn)], ...]]
# denote the matrix of [Dx(b1),..,Dx(bn)] as Mx, then (b1,..,bn)*Mx=(Dx(b1),..,Dx(bn))
# i.e., i'th colunm of Mx are coeff of Dx(bi) w.r.t. the basis {b1,..,bn}
 
F1 := proc(a,b1,b2,c,x,y)
	options remember;
	local L,aa,cc,bb,v,i;
	if type(a,`+`) then
		L := [op(a)];
		for i in L do if type(i,posint) then
			aa := a-1;
			return (
			b1*procname(aa,b1+1,b2,c,x,y)+
			b2*procname(aa,b1,b2+1,c,x,y)+
			(aa-b1-b2)*procname(aa,b1,b2,c,x,y))/aa
		fi od:
	fi;
	if type(c,`+`) then
		L := [op(c)];
		for i in L do if type(i,posint) then
			cc := c-1;
			return (
				(-y*a*x+y*cc*x-b2*x-y*b1)*procname(a,b1,b2,cc,x,y)
				-b2*x*(y-1)*procname(a,b1,b2+1,cc,x,y)
				-y*b1*(x-1)*procname(a,b1+1,b2,cc,x,y)
			) * cc/x/(a-cc)/y/(b1-cc+b2)
		fi od
	fi;
	if type(b1,`+`) and type(b2,`+`) and add(`if`(type(i,posint),1,0),i=b1)*add(`if`(type(i,posint),1,0),i=b2)>0 then
		return	(-y*procname(a,b1-1,b2,c,x,y)+procname(a,b1,b2-1,c,x,y)*x)/(-y+x)
	fi;
	if type(b1,`+`) then
		L := [op(b1)];
		for i in L do if type(i,posint) and i>1 then
			bb := b1-2;
			return (
				(y*bb+y-c*y+b2*y)*procname(a,bb,b2,c,x,y)+
				(-2*y+y*x+c*y-b2*x-b2*y+y*bb*x-y*a*x+b2*x*y-2*y*bb)*procname(a,bb+1,b2,c,x,y)+
				(-b2*x*y+b2*x)*procname(a,bb+1,b2+1,c,x,y)
			)/y/(x-1)/(bb+1)
		fi od
	fi:
	if type(b2,`+`) then
		L := [op(b2)];
		for i in L do if type(i,posint) and i>1 then
			bb := b2-2;
			return (
				(b1*x-c*x+bb*x+x)*procname(a,b1,bb,c,x,y)+
				(-y*b1-b1*x+y*b1*x-2*bb*x+bb*x*y-y*a*x-2*x+c*x+y*x)*procname(a,b1,bb+1,c,x,y)+
				(-y*b1*x+y*b1)*procname(a,b1+1,bb+1,c,x,y)
			)/x/(y-1)/(bb+1)
		fi od
	fi:
	'F1'(args)
end:

D1F1 := proc(a,b1,b2,c,x,y)
	a*b1/c * F1(a+1, b1+1, b2, c+1, x, y)
end:

D2F1 := proc(a,b1,b2,c,x,y)
	a*b2/c * F1(a+1, b1, b2+1, c+1, x, y)
end:






# Input:  parameters in Appell function F1
# Output: C(x,y)[Dx,Dy]-module of F1(a,b1,b2,c,x,y)

ModuleF1:=proc(a,b1,b2,c,x::name,y::name)
   #options trace;
   local Bas, i, B, bas, DxBas, DyBas,M,Mx,My;
   Bas[1] := F1(A, B1,   B2,   C, x, y);  
   Bas[2] := F1(A, B1+1, B2,   C, x, y);
   Bas[3] := F1(A, B1,   B2+1, C, x, y);
   bas:={seq(Bas[i],i=1..3)};  # basis of module of F1
   for i to 3 do
       DxBas[i]:=collect(D1F1(op(Bas[i])),bas,normal); # derivative of basis w.r.t. x
       DyBas[i]:=collect(D2F1(op(Bas[i])),bas,normal); # derivative of basis w.r.t. y
   od;
   M:=subs(A=a,B1=b1,B2=b2,C=c,subs(seq(Bas[i]=BAS||i,i=1..3),[[seq(Bas[i],i=1..3)],[x,y],[Dx,Dy],[[seq(DxBas[i],i=1..3)],[seq(DyBas[i],i=1..3)]]]));
   Mx:=Matrix([seq([seq(coeff(M[4][1][j],BAS||i),j=1..3)],i=1..3)]); # derivative matrix w.r.t. x
   My:=Matrix([seq([seq(coeff(M[4][2][j],BAS||i),j=1..3)],i=1..3)]); # derivative matrix w.r.t. y
   return simplify([M[2],M[3],[Mx,My]]);
end:



# Input: a module M with two variables, two rational functions u and v in two variables
# Output: the new module from M by changing two varibles to u and v

changeVmoduleF1:=proc(M,u,v)
 #options trace;
      local i, j, ind,newM,indM,M1,M2,m,mm;
      for j to 2 do m[j]:=subs(M[1][1]=l[1],M[1][2]=l[2],M[3][j]); mm[j]:=subs(l[1]=u,l[2]=v,m[j]); od;      
      ind:=indets(u) union indets(v); 
      for i in ind do for j in ind minus {i} do if has(i,j) then ind:=ind minus {i} fi; od; od; 
      
      for i to nops(ind) do
          newM[i]:=simplify(mm[1]*diff(u,ind[i])+mm[2]*diff(v,ind[i]));  # chain rule 
      od; 
      return [[op(ind)],[seq(D||i,i=ind)],[seq(newM[i],i=1..nops(ind))]];  
end:





# module2f1 deals with the case when t is a variable
# Input:  parameters in 2f1 function
# Output: C(t)[Dt]-module of 2F1(a,b,c,t)

module2f1:=proc(a,b,c,t::name)
   local L,DB1,DB2,f,i,eq,s,M,v1,v2;
   L:=t*(1-t)*Dt^2+(c-(a+b+1)*t)*Dt-a*b;
   DB1[0]:=hypergeom([a,b],[c],t);         # basis element of module of 2F1
   DB2[0]:=hypergeom([a+1,b+1],[c+1],t);
   f:=hypergeom([a+2,b+2],[c+2],t);
   for i to 2 do 
       DB1[i]:=diff(DB1[i-1],t)
   od;
   DB2[1]:=diff(DB2[0],t);
   eq:=add(coeff(L,Dt,i)*DB1[i],i=0..2); 
   s:=solve(eq,f); 
   DB2[1]:=collect(subs(f=s,DB2[1]),{DB1[0],DB2[0]},normal);
   M:=subs(DB1[0]=bas1,DB2[0]=bas2,[[DB1[0],DB2[0]],[t],[D||t],[DB1[1],DB2[1]]]);
   v1:=Vector([coeff(M[4][1],bas1),coeff(M[4][1],bas2)]);
   v2:=Vector([coeff(M[4][2],bas1),coeff(M[4][2],bas2)]);
   return [M[2],M[3],[Matrix([v1,v2])]];
end:



# input: a module M with one varible and any rational functions u in that varible
# output: new module from M by send the varible to the function u  

changeVmodule2F1:=proc(M,u)
 #options trace;
      local m1, m2,i,ind,newM,indM,M1,M2;
      indM:=M[1]; 
      m1:=M[3][1]; 
      ind:=[op(indets(u))];
      M1:=subs(indM[1]=u,m1);
      for i to nops(ind) do
          newM[i]:=simplify(M1*diff(u,ind[i]));
      od; 
      return [ind,[seq(D||i,i=ind)],[seq(newM[i],i=1..nops(ind))]];  
end:



# Apply the i'th derivative of m, m is a vector, (c1, ..., cn)'(so m=c1*bas1+...+cn*basn), M is the module
# Input: module M, number i and a vector m in M.
# Output: i^th derivative of m.

ApplyDi := proc(M,i,m)  local j,A,B,Mi,var;
        Mi:=M[3][i]; var:=M[1][i];
        A:=map(diff,m,var); B:= LinearAlgebra[Multiply](Mi,m); simplify(A+B);
end:


# Input: a module M, an operator L and an element m in M
# Output: result of applying L on m
# m must be a column vector denoting an element of the module M, like (c1, ..., cn)' (so m=c1*B1+...+cn*Bn)

ApplyOp := proc(M, L, m) 
#options trace;
local Dv, i,j, S, d; 
        Dv := indets(L,name) intersect {op(M[2])};
        if Dv={} then
                LinearAlgebra[ScalarMultiply](m, L)  # L is a constant
        else
                Dv := Dv[1];     
                for i do if M[2][i] = Dv then break fi od;
                S := procname(M, coeff(L, Dv, 0), m);  
                d := m;
                for j to degree(L,Dv) do
                        d := ApplyDi(M, i, d);
                        S := S + LinearAlgebra[Multiply](coeff(L, Dv, j),d)
                od;
                S
        fi
end:


# Input: a module M and a matrix L
# Output: L must be a matrix
FindLinearRelation := proc(M,L)
 #options trace;
        local nrow, ncol, b,r,s;
        nrow:=LinearAlgebra[RowDimension](L); ncol:=LinearAlgebra[ColumnDimension](L);# ncol=nrow+1
        b:=[seq(i*0,i=1..nrow)];
        r:=linalg[linsolve](L,Vector(b));
        r:=[seq(r[i],i=1..ncol)];
        s:=indets(r) minus indets(M);  
        if nops(s)<>1 then  # nops(s)=1 then the rank of columns of L = nrow = ncol-1, so m is cyclic
           return NULL; else return eval(r, s[1]=1) 
        fi;        
end:



# Input: M = module, m = element.
# Output: the minimal operator of m w.r.t. the i'th variable
MinOpM := proc(M, m, i) local n, L, j;
        n := LinearAlgebra[RowDimension](M[3][1]);
        L[0] := m;
        for j to n do L[j] := ApplyDi(M,i, L[j-1]) od:
        L:=Matrix([seq(L[j],j=0..n)]);
        L := FindLinearRelation(M,L);
        add(L[j+1]*M[2][i]^j, j=0..n)
end:


# Input: module M, d and f are matrices, like (c1,..,cn)'
# Output: result of d applying on f
ApplyDualElement := proc(M, d, f) 
   #options trace;
   local m;
   m:=LinearAlgebra[Transpose](d);
   LinearAlgebra[Multiply](m,f);
end:



# Input: module M
# Output: its dual module
dualmodule:=proc(M)
   #options trace;
   local n,i,m;
   n:=nops(M[3]);
   for i to n do
       m[i]:=(-1)*LinearAlgebra[Transpose](M[3][i])
   od; return [M[1],M[2],[seq(m[i],i=1..n)]]
end:



# Input: m is a function, list L gives all variables in m
# Output: the module of m with form [[B],L,DL,[seq(diff(m,L[i])*B],i=1..nops(L))]
OneDiModule:=proc(m,L)
    local l,i,r,Dm;
    l:=[seq(D||i,i=L)];
    for i to nops(l) do
        r[i]:= simplify(diff(m,L[i])/m);
    od;
    Dm:=[seq([r[i]],i=1..nops(l))];
    return [L,l,Dm];
end:



# Input: two modules M1 and M2, i'th derivative, j'th element in M1 basis, k'th element in M2 basis
# Output: tensor product of bj and Bk
ApplyDiTensor:=proc(M1,M2,i,j,k)
     local N,M,colmj1,colmk2,s1,s2,A1,A2,m;
     N:=LinearAlgebra[RowDimension](M1[3][1]); M:=LinearAlgebra[RowDimension](M2[3][1]);
     if i>nops(M1[1]) or j>N or k>M then return "error" fi;
     colmj1:=LinearAlgebra[Column](M1[3][i],j);
     colmk2:=LinearAlgebra[Column](M2[3][i],k);
     s1:=[seq(colmj1[m],m=1..N)];
     s2:=[seq(colmk2[m],m=1..M)];
     A1:=[0$M*N];A2:=[0$M*N];
     for m from 0 to N-1 do
         A1:=subsop(m*M+k=s1[m+1],A1)
     od;
     for m from 1 to M do
         A2:=subsop((j-1)*M+m=s2[m],A2)
     od;
     Vector(A1)+Vector(A2);
end:




# Input: two modules M1 and M2
# Output: tensor product module of the input
TProModule:=proc(M1,M2)
    local n,m,basis,i,j,k,Dibasis,D;
    if M1[1]<>M2[1] then return "error" fi;
    n:=LinearAlgebra[RowDimension](M1[3][1]); m:=LinearAlgebra[RowDimension](M2[3][1]); 
    for i to nops(M1[1]) do  
        for j to n  do
            for k to m do
                Dibasis[j,k]:=ApplyDiTensor(M1,M2,i,j,k)
            od;
        od;   D[i]:=Matrix([seq(seq(Dibasis[j,k],k=1..m),j=1..n)]);
    od;
    return [M1[1],M1[2],[seq(D[i],i=1..nops(M1[1]))]]
end:     



summodule:=proc(M1,M2)
   local RN1,RN2, CN1,CN2,A1,A2,Mat,i,j,k,r;
   if M1[1]<>M2[1] then return "error" fi;
   RN1:=LinearAlgebra[RowDimension](M1[3][1]);
   RN2:=LinearAlgebra[RowDimension](M2[3][1]);
   CN1:=LinearAlgebra[ColumnDimension](M1[3][1]);
   CN2:=LinearAlgebra[ColumnDimension](M2[3][1]);
   for i to nops(M1[1]) do 
       A1:=M1[3][i];
       A2:=M2[3][i];
       for j to RN1 do
           r[j]:=[seq(LinearAlgebra[Row](A1,j)[k],k=1..CN1),0$CN2];
       od;
       for j to RN2 do
           r[RN1+j]:=[0$CN1,seq(LinearAlgebra[Row](A2,j)[k],k=1..CN2)];
       od;
       Mat[i]:=Matrix([seq(r[j],j=1..RN1+RN2)]);
   od;
   return [M1[1],M1[2],[seq(Mat[i],i=1..nops(M1[1]))]];
end:



      

## Input: two matrices 
## Output: a matrix C s.t. A=B*C

findC:=proc(A,B) # A and B are homomorphism matrices of #basis by dim. try to find dim by dim matrix C s.t. A=B*C.
  #options trace;
  local dim, RN, r, C, D, s, detA, detB, p;
  dim:=LinearAlgebra[ColumnDimension](A); RN:=LinearAlgebra[RowDimension](A);
  if dim=RN #number of basis equals dimension, square matrix of homx
         then detA:=LinearAlgebra[Determinant](A);
              detB:=LinearAlgebra[Determinant](B);
              if detA*detB<>0  
                 then C:=LinearAlgebra[Multiply](LinearAlgebra[MatrixInverse](B),A);
                 elif (detA=0 and detB<>0) or (detA<>0 and detB=0) 
                      then return "not isomorphic"
              fi;
  fi;
  if dim<>RN or (dim=RN and detA=0 and detB=0) then
      for p to dim do
          r[p]:=Matrix([seq([g[q+(p-1)*dim]],q=1..dim)])
      od;
      C:=Matrix([seq(r[p],p=1..dim)]);
      D:=LinearAlgebra[Multiply](B,C)-A;
      s:=solve({seq(seq(numer(D[p][q]),p=1..RN),q=1..dim)},{seq(g[k],k=1..dim*dim)});
      if s=NULL then return "not isomorphic"
         else C:=subs(s,C);
      fi;
  fi;
  C;
end:

                  
                      


# Input: a module M, an element m in M and an indet number i. 
# Output: the minimal operator of m w.r.t. ith variable if m is cyclic w.r.t ith variable.
# Test if a given element is cyclic w.r.t the given variable. If yes, then compute its minimal operator.
IsCyclicI:=proc(M,m,i) 
     #options trace;
       local n,j,r,L,Op;
       L[0]:= m;
       n:=LinearAlgebra[RowDimension](M[3][1]);
       for j to n do 
           L[j]:=ApplyDi(M,i,L[j-1]);         
       od;
       L:=Matrix([seq(L[j],j=0..n)]);
       L := FindLinearRelation(M,L);
       if L<>NULL then 
             Op:= add(L[j+1]*M[2][i]^j, j=0..n);     # minimal operator of m w.r.t. i^th variable          
             return collect(Op/coeff(Op,M[2][i]^degree(Op,M[2][i])),{M[2][i]},factor); #coeff of highest degree of Dx^i = 1
             else return NULL; 
       fi;
end:


# Input: a module M and a variable 
# Output: a cyclic vector in M w.r.t. that variable and its minimal operator w.r.t. it.

CycVec:=proc(M,var)
       #options remember;
       #options trace;
       local i,j,m,n,N,v,roll,indet,e,ind;
       n:=LinearAlgebra[RowDimension](M[3][1]); # rank of M as D-module
       for i to nops(M[1]) do
           if M[1][i]=var then ind:=i; break; fi;  # get the index of var -- ind
       od;
       i:=0;
       while i<=n-1 do
           j:=n-1-i;
           m:=Matrix(Vector([0$i,1,0$j]));    # basis elements 
           if IsCyclicI(M,m,ind)<>NULL 
              then return([m,IsCyclicI(M,m,ind)]); break;          
           fi;
           i:=i+1;
       od;
       roll:=rand(-2..2);
       for i to n do
           v[i]:=roll();
       od;
       m:=Matrix(Vector([seq(v[i],i=1..n)]));
       if IsCyclicI(M,m,ind)<>NULL 
              then return [m,IsCyclicI(M,m,ind)]; break;          
       fi;
       roll:=rand(-9..9);
       for i to n do
           v[i]:=roll();
       od;
       m:=Matrix(Vector([seq(v[i],i=1..n)]));
       if IsCyclicI(M,m,ind)<>NULL 
              then return [m,IsCyclicI(M,m,ind)]; break;          
       fi;
       indet:=M[1][1];
       for i to n do
           v[i]:=randpoly([x],degree=1,coeffs=rand(-2..2));
       od;
       m:=Matrix(Vector([seq(v[i],i=1..n)]));
       if IsCyclicI(M,m,ind)<>NULL 
              then return [m,IsCyclicI(M,m,ind)]; break;          
       fi;
end:




# Input: a list of d n' by n matrices l
# output: one n'*n by d matrix A. 
rewriteMat:=proc(l)  
 #options trace;
   local d, n, k, a, A, nprime;
   d:=nops(l);  #dimension of hx and hy
   nprime:=LinearAlgebra[RowDimension](l[1]); #rank of the second module
   n:=LinearAlgebra[ColumnDimension](l[1]);
   for k to d do
       a[k]:=Matrix([seq(seq([l[k][i][j]],j=1..n),i=1..nprime)]); # rewrite every n' by n matrix into a n'*n by 1 matrix
   od;
   A:=Matrix([seq(a[i],i=1..d)]);  # combine d matrices to one matrix
   A;
end:


# Input:  an n'*n by d matrix A, n' and n.
# Output: a list of d n' by n matrix. Each matrix is a colomn of A.
MatBack:=proc(A,nprime,n)
   local r,d,i,k,s,m,B;
   r:=LinearAlgebra[RowDimension](A);
   d:=LinearAlgebra[ColumnDimension](A);
   if r<>nprime*n then return "error" fi;
   for i to d do
       for k to nprime do
           s:=1+(k-1)*n;
           m[i][k]:=Matrix([seq(A[j][i],j=s..s+n-1)]);
       od;
       B[i]:=Matrix([seq([m[i][k]],k=1..nprime)]);
   od;
   return [seq(B[i],i=1..d)];
end:



# Input: A and B are n' by n matrices. l1, l2, l are sets of variables, A is unique up to l1 and B is unique up to l2. l is the set of all variables.
# Output: a list of a new matrix D and a new variable list ll such that D is equal to A and B up to their own constants and D is unique up to the variable list ll.
# First check if A=B*C, i.e., check if all quotients(n*n') of entries in A dividing entries in B are equal.
# Then check if C=C1*C2 where the indets of elements in C1 is in l2 and indets of elements in C2 is in l1.
# return A/C2 and the variable set l1 intersect l2. 

combHom:=proc(A,B,l1,l2,l) 
 #options trace;
   local RN, CN, Zero, var1, var2, i, j, q, Numer, Denom, xn, yn, xd, yd, Q;
   RN:=LinearAlgebra[RowDimension](A);       # number of rows
   CN:=LinearAlgebra[ColumnDimension](A);    # number of columns
   Zero:=Matrix([seq([seq(i*j*0,i=1..CN)],j=1..RN)]);   
   var1:= l minus l1; var2:=l minus l2; 
   for i to RN do
       for j to CN do
           if simplify(B[i][j])<>0 then q:=normal(A[i][j]/B[i][j]); break; fi; 
       od;
   od;
   if not LinearAlgebra[Equal](simplify(A-q*B),Zero) then return {} fi;
   Numer:=numer(q);
   Denom:=denom(q);
   xn:=content(Numer,var1);
   yn:=content(Numer,var2);
   xd:=content(Denom,var1);
   yd:=content(Denom,var2);
   Q:=normal((q)*(xd)*(yd)/(xn)/(yn));
   if (var1 union var2) intersect indets(Q)<>{}
          then return {}; 
          else return [simplify(xd*A/xn), l1 intersect l2];   
   fi;
end:




# Input: A,B are n'*n by d matrices with d>1, n' is the dimension of the second module, n is the dimension of the first module and d is the demension of hx and #        hy. For every hx in univariate homomorphisms, we have one n' by n matrix, [hx(m'),Dx(hx(m')),...]*[m,Dx(m),...]^(-1), A is obtained by rewriteMat to      #        combine all such matrices. Likewise for B. l1, l2, l are sets of variables, A is unique up to l1 and B is unique up to l2. l is the set of all variables.
# Output: a list of a new matrix and a new variable list ll such that the matrix is equal to A and B up to their constants and unique up to the variable list ll.     # First use "findC" to find the d by d matrix C s.t. A=B*C.
# Then give some value to the variable in l1\l2 in C s.t. the new matrix from C is invertible, let D1 be it and let D2 be D1^(-1)*C.                          
# Return A*(D2^(-1)) and varable set l1 intersect l2.

combHom2:=proc(A,B,l1,l2,l)  
   #options trace;
   local C, ind, ind2, var, i, D1, D2;
   C:=findC(A,B);
   ind:=l1 minus l2; ind2:=l2 minus l1;
   if nops(ind)=1 then 
      var:=op(ind);
      i:=1; D1:=Matrix([[1,0],[1,0]]);
      while Determinant(D1)=0 do
          D1:=subs(var=i,C);
          i:=i+1;
      od;
      D2:=Multiply(MatrixInverse(D1),C);
   elif nops(ind2)=1 then
      var:=op(ind2);
      i:=1; D2:=Matrix([[1,0],[1,0]]);
      while Determinant(D2)=0 do
          D2:=subs(var=i,C);
          i:=i+1;
      od;
      D1:=Multiply(C,MatrixInverse(D2));
   fi;
   if `subset`(indets(D1),l2) and `subset`(indets(D2),l1) 
       then return [Multiply(A,MatrixInverse(D2)),l1 intersect l2]
       else return {}
   fi;
end:





# Input: two modules M1 and M2 with format [[variables],[derivatives],[matrices of derivatives]]. ([[x1,...,xp][Dx1,...,Dxp],[Mx1,...,Mxp]])
# Output: the homomorphism(s) between M1 and M2 if exist(s).
# 1. For each variable and each module, use "CycVec" to get the cyclic vector and its minimal operator w.r.t. the variable.
# 2. For each variable, use "Homomorphisms" in DEtools package to find the univariate homomorphisms between two modules.
# 3. For each variable, say x1, every element in the univariate homomorphisms, say hx1, gives one n' by n matrix [hx(m'),Dx(hx(m')),...]*[m,Dx(m),...]^(-1).
# 4. If d1=d2=...=dp=1, di is the dimension of the univariate homomorphism w.r.t. xi. Then use combHom to combine all of matrices by two matrices per time.
# 5. If the C(x1,...,xp)-vector spaces of univariate homomorphisms are equal and d1>1, then for each variable, rewrite the homomoprhisms into one n'*n by d matrix and #    then use combHom2 to combine all of them like 4.
# Comment. This "hom" can only be developed for two special cases.

hom := proc(M1, M2) 
 local varN,vars,nprime, n, i, var, Cyc1,Cyc2,m1,m2,j,m1prime,Dm1,Mat2,minop1,minop2,L,dim,k,hxm2,hxm2prime,Dhxm2,Mat1,H,Mat,G,d,checkdim,hh;
       #options trace;   
        varN:=nops(M1[1]);   # number of variables
        vars:={op(M1[1])};   # variables
        nprime:=LinearAlgebra[RowDimension](M2[3][1]);   # rank of M2 as a D-module
        n:=LinearAlgebra[RowDimension](M1[3][1]);        # rank of M1 as a D-module
        for i to varN do                                 # for each variable do
                 var:=M1[1][i];                          # i^th variable
                 Cyc1:=CycVec(M1,var);     # cyclic vector of M1 w.r.t. var and its minimal operator
                 Cyc2:=CycVec(M2,var);     # cyclic vector of M2 w.r.t. var and its minimal operator
                 m1:=Cyc1[1]; m2:=Cyc2[1]; # cyclic vectors
                 j:=1; 
                 m1prime:=m1; 
                 Dm1[0]:=m1; 
            while j<n do
                  Dm1[j]:=ApplyOp(M1,M1[2][i],m1prime);    
                  m1prime:=Dm1[j];
                  j:=j+1;
            od;
                  Mat2:=Matrix([seq(Dm1[j],j=0..n-1)]);       # [m,Dx(m),...]

                  minop1:=Cyc1[2]; minop2:=Cyc2[2];           # minimal operators w.r.t. m1 and m2
                  _Envdiffopdomain := [M1[2][i], M1[1][i]];
                  L := DEtools[Homomorphisms](minop2,minop1);               # homomorphisms between minop1 and minop2
                  if L=[] then return {} fi;
                  dim[i]:=nops(L);                                          # dimension of the univariate homomorphism
            for k to nops(L) do  # each univariate homomorphism corr. to a homomorphisms from M1 to M2 as Dxi-modules
                       hxm2:=ApplyOp(M2,L[k],m2);
                       j:=1; hxm2prime:=hxm2; Dhxm2[0]:=hxm2;
                  while j<n do
                       Dhxm2[j]:=ApplyOp(M2,M1[2][i],hxm2prime);
                       hxm2prime:=Dhxm2[j];          
                       j:=j+1;
                  od;
                       Mat1:=Matrix([seq(Dhxm2[j],j=0..n-1)]);         # (h(m2),Dxi(h(m2)),...)
                       H[k]:=simplify(LinearAlgebra[Multiply](Mat1,LinearAlgebra[MatrixInverse](Mat2))); 
           od;
           if dim[i]=1 then 
                  G[i]:=[H[1],{seq(M1[1][j],j=1..i-1),seq(M1[1][j],j=i+1..nops(M1[1]))}]; #d=1 just return the n' by n matrix and its constants
                  else 
                       Mat[i]:=rewriteMat([seq(H[k],k=1..nops(L))]);    # write all n' by n matrices into an n*n' by d matrix, d>1 is the dim(univariatehom)
                       G[i]:=[Mat[i],{seq(M1[1][j],j=1..i-1),seq(M1[1][j],j=i+1..nops (M1[1]))}]; # n'*n by d matrix and its constants
           fi;
        od;
        d:=dim[1]; checkdim:=1; 
        for i from 2 to varN do 
            if d<>dim[i] then checkdim:=0; break;fi;
        od; 
        if checkdim=1 and d=1 then                             # the easiest case when all dimensions of univariate homomorphisms equal 1.
          if varN=1 then return G[1][1]; else 
           hh:=combHom(G[1][1],G[2][1],G[1][2],G[2][2],vars);
           if varN=2 then return hh[1];
              else for i from 3 to varN do
                   hh:=combHom(hh[1],G[i][1],hh[2],G[i][2],vars);
              od;
              return hh[1];
           fi;
          fi;
        fi;
        if checkdim=1 and d>1 then    #and SPAN(Mat[i])=SPAN(Mat[j]) for any i,j. The case when all dimension of uni.hom equal but greater than 1. 
             hh:=combHom2(G[1][1],G[2][1],G[1][2],G[2][2],vars);
             if varN=2 then return MatBack(hh[1],nprime,n);
                else for i from 3 to varN do
                     hh:=combHom2(hh[1],G[i][1],h[2],G[i][2],vars);
                od;
                return MatBack(hh[1],nprime,n);
             fi;
        fi;
end: