Clear[statevariables, timeflux,  spaceflux,  Dop,  LB,  IntCoef, 
      TimeFluxMx, SpaceFluxMx,  dimofsystem,  FMx,  EVal,  REVMx,  
      LEVMx,  REV,  LEV,  LBNum,  colnumber, rownumber,  LBMx, LBCMx, 
      IntCoefNum]
\(Print["These are the state variables, time/space flux functions and the dimension of the system"]

 (*P - 
        System*) 
 (*statevariables = \({u, 
              v}
timeflux = \({u, 
                v}
spaceflux = {P[
                v], \(-u\)}\)\)*) 

 (*ahat\ system*) 
 (*statevariables = \({u, v, 
              w}
timeflux = \({u, v, 
                w}
spaceflux = {a[w]*u + b[w]*v, 
              c[w]*u + d[w]*v, 
              f[w]}\)\)*) 

 (*Euler\ Gas\ Dynamics*) 
 (*statevariables = \({p, u, 
              e}
timeflux = \({p, p*u, 
                p*u^2/2 + p*e}
spaceflux = {p*u, \ 
              p*u^2 + p, 
              p*u^3/2 + p*u*e + 
                u*p}\)\)*) 

 (*Maxwell' 
        s\ system*) 
 (*statevariables = \({B1, B2, B3, E1, 
              E2, E3}
timeflux = \({B1, B2, B3, E1, E2, 
                E3}
spaceflux = {0, \(-E3\)/ϵ, 
              E2/ϵ, 0, 
              B3/μ, \(-B2\)/μ}\)\)*) 

 (*String\ System*) 
 (*statevariables = \({U1, U2, U3, V1, V2, 
              V3}
timeflux = \({U1, U2, U3, V1, V2, 
                V3}
\((1 - 
                  1/\((Sqrt[
                        U1^2 + U2^2 + 
                          U3^2])\))\)
spaceflux = {V1, V2, 
              V3, \((1 - 1/\((Sqrt[U1^2 + U2^2 + U3^2])\))\) 
                U1, \((1 - 1/\((Sqrt[U1^2 + U2^2 + U3^2])\))\) 
                U2, \((1 - 1/\((Sqrt[U1^2 + U2^2 + U3^2])\))\) 
                U3}\)\)*) 

 (*Template*) 
 (*
statevariables\  = \ \({}
timeflux\  = \ \({}
spaceflux\  = \ {}\)\)
*) 
\)
dimofsystem = Length[statevariables]
Dop[f_] := 
      Transpose[
        Table[(∂_(Part[statevariables, i])) f, {i, 1, 
            dimofsystem}]]
LB[U_, V_] := Dop[U] . V - Dop[V] . U
IntCoef[U_, V_, W_] := 
      U . LB[V, W] 

Print["Time-Flux-matrix, Space-Flux-matrix:"]
TimeFluxMx = Dop[timeflux];
SpaceFluxMx = Dop[spaceflux];
MatrixForm[TimeFluxMx]
MatrixForm[SpaceFluxMx]

Print["This is the quasi-linear flux matrix:"]
FMx = 
        Simplify[
          Inverse[TimeFluxMx] . SpaceFluxMx];
MatrixForm[FMx]

Print["Eigenvalues, Right Eigenvectors, Left Eigenvectors:"]
EValVec = Eigenvalues[FMx];
REVMx = Transpose[Eigenvectors[FMx]];
LEVMx = Inverse[REVMx];
MatrixForm[Simplify[EValVec]]
MatrixForm[Simplify[REVMx]]
MatrixForm[Simplify[LEVMx]]
EVal[n_] := Part[EValVec, n]
REV[n_] := Part[Transpose[REVMx], n]
LEV[n_] := Part[LEVMx, n]

Print["Lie Bracket Matrix, Lie Bracket Coefficient Matrix"]
LBMx = 
        Table[Table[
            Simplify[LB[REV[colnumber], REV[rownumber]]], {colnumber, 1,  
              dimofsystem}], {rownumber, 1, 
            dimofsystem}];
LBNum[n1_, n2_] := Part[Part[LBMx, n1], n2]
LBCMx = 
        Table[Table[
            Simplify[
              LinearSolve[REVMx, LBNum[colnumber, rownumber]]], {colnumber, 
              1,  dimofsystem}], {rownumber, 1, 
            dimofsystem}];
MatrixForm[LBMx]
MatrixForm[LBCMx]

Print["Tensor of interaction coefficients"]
IntCoefNum[n1_, n2_, n3_] := 
      IntCoef[LEV[n1], REV[n2], REV[n3]]
MatrixForm[
      Table[Table[
          Table[Simplify[IntCoefNum[lev, rev1, rev2]], {lev, 1, 
              dimofsystem}], {rev1, 1,  dimofsystem}], {rev2, 1, 
          dimofsystem}]]