`type/mytype`:= proc (v) if type(v, indexed) and op(0, v) = E or op(0, v) = F or op(0, v) = Theta then return true else return false end if end proc:

`type/mutype`:= proc (v) if type(v, indexed) and op(0, v) = mu then return true else return false end if end proc:

 Arrange := proc (v) local L, m; if v = 0 then return [0, 0, 0] else L := indets(v, mytype); m := nops(L) end if; if 1 < m then return [m, L[1], add(i, i = normal([op(collect(v, L, distributed))]))] elif m = 1 then return 1, L[1], coeff(v, L[1])*L[1] end if end proc:
 
 grade1 := proc (V) local L; L := Arrange(V); if L[1] = 1 and L[2] <> E[10000000000, 10000000000] and L[2] <> F[10000000000, 10000000000] then return 2*ldegree(coeff(L[3], L[2]), indets(L[3], mutype))+op(2, L[2]) end if end proc:
 
 GRADE := proc (V) local J, w; return max(map(grade1, {op(V+E[10000000000, 10000000000]+F[10000000000, 10000000000])})) end proc:
 
 ELKTRASFORMATION := proc (V1, l, k, a0) local c; c := 0; if l < k then c := E[l+1, k]*coeff(V1, E[l, k])/((2*l+2)*a0) end if; return c end proc:
 
 FLKTRASFORMATION := proc (V1, l, k, a0) local d; d := 0; if l < k then d := (1/2)*F[l+1, k]*coeff(V1, F[l, k])/(a0*(l+2)) end if; return d end proc:
 
 ThetaLKTRASFORMATION:= proc (V1, l, k, a0) local e; e := 0; if l < k then e := (1/2)*Theta[l+1, k]*coeff(V1, Theta[l, k])/(a0*(l+1)) end if; return e end proc:
 
TRANSFORMATION := proc (n, VF) local transformation, trans, i, W, INDETS, NEWVF, a0, a1; transformation := 0; trans := 0; NEWVF := VF+E[20000, 20000]+F[2000, 2000]; a0 := coeff(VF, F[-1, 0]); a1 := subs(mu[1] = 0, mu[2] = 0, mu[3] = 0, mu[4] = 0, mu[5] = 0, mu[6] = 0, mu[7] = 0, mu[8] = 0, mu[9] = 0, mu[10] = 0, mu[11] = 0, a0); for i to nops(NEWVF) do if E[20000, 20000] = op(i, NEWVF) or F[2000, 2000] = op(i, NEWVF) then trans := 0; transformation := transformation+trans else W := op(i, NEWVF); INDETS := indets(W, mytype); if n = GRADE(W) then if op(0, INDETS[1]) = E then trans := ELKTRASFORMATION(op(i, NEWVF), op(1, INDETS[1]), op(2, INDETS[1]), a1); transformation := transformation+trans elif op(0, INDETS[1]) = F then trans := FLKTRASFORMATION(op(i, NEWVF), op(1, INDETS[1]), op(2, INDETS[1]), a1); transformation := transformation+trans elif op(0, INDETS[1]) = Theta then trans := ThetaLKTRASFORMATION(op(i, NEWVF), op(1, INDETS[1]), op(2, INDETS[1]), a1); transformation := transformation+trans end if else trans := 0; transformation := transformation+trans end if end if end do; return transformation end proc:
 
 
 LIEBRACKETEE := proc (l, k, m, n) return (n-k)*E[l+m, k+n] end proc:
 
 LIEBRACKETFF := proc (l, k, m, n) return ((m+1)*(k+2)-(l+1)*(n+2))*F[l+m, k+n] end proc:
 
 LIEBRACKETThetaTheta:= proc (l, k, m, n) return 0 end proc:
 
 LIEBRACKETFE := proc (l, k, m, n) return (n+2)*(m*(k+2)-n*(l+1))*E[l+m, k+n]/(k+n+2)-k*(k+2)*F[l+m, k+n]/(k+n+2) end proc:
 
 LIEBRACKETEF := proc (m, n, l, k) return (n+2)*(-m*(k+2)+n*(l+1))*E[m+l, n+k]/(n+k+2)+k*(k+2)*F[m+l, n+k]/(n+k+2) end proc:
 
LIEBRACKETFTheta:= proc (l, k, m, n) return (m*(k+2)-n*(l+1))*Theta[l+m, k+n] end proc:

LIEBRACKETThetaF:= proc (m, n, l, k) return (-m*(k+2)+n*(l+1))*Theta[m+l, n+k] end proc:

 LIEBRACKETETheta:= proc (l, k, m, n) return n*Theta[l+m, k+n] end proc:
 
LIEBRACKETThetaE:= proc (m, n, l, k) return -n*Theta[m+l, n+k] end proc:
 
 LIEPRODUCT := proc (V1, V2) local INDETS1, INDETS2; 
INDETS1 := indets(V1, mytype); INDETS2 := indets(V2, mytype);
 if V1 = 0 or V2 = 0 then return 0 
 elif op(0, INDETS1[1]) = E and op(0, INDETS2[1]) = E then 
 return coeff(V1, INDETS1[1], 1)*coeff(V2, INDETS2[1], 1)*LIEBRACKETEE(op(1, INDETS1[1]), op(2, INDETS1[1]), op(1, INDETS2[1]), op(2, INDETS2[1])) 
 elif op(0, INDETS1[1]) = F and op(0, INDETS2[1]) = F then
  return coeff(V1, INDETS1[1], 1)*coeff(V2, INDETS2[1], 1)*LIEBRACKETFF(op(1, INDETS1[1]), op(2, INDETS1[1]), op(1, INDETS2[1]), op(2, INDETS2[1])) 
  elif op(0, INDETS1[1]) = Theta and op(0, INDETS2[1]) = Theta then return coeff(V1, INDETS1[1], 1)*coeff(V2, INDETS2[1], 1)*LIEBRACKETThetaTheta(op(1, INDETS1[1]), op(2, INDETS1[1]), op(1, INDETS2[1]), op(2, INDETS2[1])) 
  elif op(0, INDETS1[1]) = E and op(0, INDETS2[1]) = F then
   return coeff(V1, INDETS1[1], 1)*coeff(V2, INDETS2[1], 1)*LIEBRACKETEF(op(1, INDETS1[1]), op(2, INDETS1[1]), op(1, INDETS2[1]), op(2, INDETS2[1])) 
   elif op(0, INDETS1[1]) = F and op(0, INDETS2[1]) = E then
    return coeff(V1, INDETS1[1], 1)*coeff(V2, INDETS2[1], 1)*LIEBRACKETFE(op(1, INDETS1[1]), op(2, INDETS1[1]), op(1, INDETS2[1]), op(2, INDETS2[1])) elif op(0, INDETS1[1]) = E and op(0, INDETS2[1]) = Theta then return coeff(V1, INDETS1[1], 1)*coeff(V2, INDETS2[1], 1)*LIEBRACKETETheta(op(1, INDETS1[1]), op(2, INDETS1[1]), op(1, INDETS2[1]), op(2, INDETS2[1])) elif op(0, INDETS1[1]) = Theta and op(0, INDETS2[1]) = E then
     return coeff(V1, INDETS1[1], 1)*coeff(V2, INDETS2[1], 1)*LIEBRACKETThetaE(op(1, INDETS1[1]), op(2, INDETS1[1]), op(1, INDETS2[1]), op(2, INDETS2[1])) elif op(0, INDETS1[1]) = F and op(0, INDETS2[1]) = Theta then
      return coeff(V1, INDETS1[1], 1)*coeff(V2, INDETS2[1], 1)*LIEBRACKETFTheta(op(1, INDETS1[1]), op(2, INDETS1[1]), op(1, INDETS2[1]), op(2, INDETS2[1])) elif op(0, INDETS1[1]) = Theta and op(0, INDETS2[1]) = F then
       return coeff(V1, INDETS1[1], 1)*coeff(V2, INDETS2[1], 1)*LIEBRACKETThetaF(op(1, INDETS1[1]), op(2, INDETS1[1]), op(1, INDETS2[1]), op(2, INDETS2[1])) end if end proc:
 
 LIEPRODUCTS := proc (f, g) local NEWF, NEWG, i, j, lie, lieproducts; lie := 0; lieproducts := 0; if f = 0 or g = 0 then return 0 else NEWF := f+E[20000, 20000]+F[2000, 2000]; NEWG := g+E[20000, 20000]+F[2000, 2000]; for i to nops(NEWF) do for j to nops(NEWG) do if E[20000, 20000] = op(i, NEWF) or E[20000, 20000] = op(j, NEWG) or F[2000, 2000] = op(i, NEWF) or F[2000, 2000] = op(j, NEWG) then lie := 0 else lie := LIEPRODUCT(op(i, NEWF), op(j, NEWG)) end if; lieproducts := lieproducts+lie end do end do; return lieproducts end if end proc:
 
 ADJOINT := proc (N, Y, VF) local i, n, N10, NEWVF, adjoint; if Y = 0 then adjoint := VF else N10 := GRADE(Y); n := floor(N/N10); NEWVF := expand(simplify(VF)); adjoint := VF; for i to n do NEWVF := LIEPRODUCTS(Y, NEWVF); adjoint := adjoint+NEWVF/factorial(i) end do end if; adjoint := expand(simplify(adjoint)); return adjoint end proc:


TRUNCATE := proc (VF, n) local i, NEWVF, W, trun, truncate; trun := 0; truncate := 0; NEWVF := expand(VF+E[20000, 20000]+F[2000, 2000]); for i to nops(NEWVF) do W := op(i, NEWVF); if W <> E[20000, 20000] and W <> F[2000, 2000] and GRADE(W) <= n then trun := op(i, NEWVF) else trun := 0 end if; truncate := truncate+trun end do; return truncate end proc:

 SECONDNORMALFORM := proc (N, VF) local newvf1, newvf2, i, secondnormalform; secondnormalform := expand(simplify(TRUNCATE(VF, N))); for i to N do newvf1 := TRANSFORMATION(i, expand(secondnormalform)); newvf2 := TRUNCATE(ADJOINT(N, newvf1, expand(secondnormalform)), N); secondnormalform := expand(newvf2) end do; return secondnormalform end proc:
 
 gradenn := proc (V, s) local L, grade, grad, W, INDETS; grade := 0; grad := 0; L := Arrange(V); if L[1] = 1 and L[2] <> E[10000000000, 10000000000] and L[2] <> F[10000000000, 10000000000] then W := L[2]; INDETS := indets(W, mytype); if op(0, INDETS[1]) = E or op(0, INDETS[1]) = F then 
 grad := s*(op(2, INDETS[1])-op(1, INDETS[1]))+op(2, INDETS[1])+ldegree(coeff(L[3], L[2]), indets(L[3], mutype)); grade := max(grade, grad) else grad := s*(op(2, INDETS[1])-op(1, INDETS[1]))+op(2, INDETS[1])+s+1+ldegree(coeff(L[3], L[2]), indets(L[3], mutype)); grade := max(grade, grad) end if end if; return grade end proc:
 
 GRADE2n := proc (V, s) local J, w; return max(map(gradenn, {op(V+E[10000000000, 10000000000]+F[10000000000, 10000000000])}, s)) end proc:

 ADJOINT2 := proc (N, Y, VF, s) local i, n, N10, NEWVF, adjoint; if Y = 0 then adjoint := VF else N10 := GRADE2n(Y, s); n := floor(N/N10); NEWVF := expand(VF); adjoint := VF; for i to n do NEWVF := LIEPRODUCTS(Y, NEWVF); adjoint := adjoint+NEWVF/factorial(i) end do end if; adjoint := expand(adjoint); return adjoint end proc:
 

 CoeffE := proc (V) local W, i, INDETS, A, VF; 
 A := 10^100; 
 VF := subs(seq(mu[i] = 0, i = 1 .. 100), expand(V)+E[1000, 1000]); 
 for i to nops(VF) do 
 W := op(i, VF); 
 INDETS := indets(W, mytype);
  if op(0, INDETS[1]) = E and op(2, INDETS[1]) <> 0 and W <> E[1000, 1000] then
   if op(1, INDETS[1]) <= A
    then A := op(1, INDETS[1]) end if end if end do; return A end proc:
 
     
 CoeffF := proc (V) local W, i, INDETS, A, VF;
  A := 10^100; 
  VF := subs(seq(mu[i] = 0, i = 1 .. 100), expand(V)+E[1000, 1000]);
   for i to nops(VF) do 
   W := op(i, VF);
    INDETS := indets(W, mytype);
     if op(0, INDETS[1]) = F and op(2, INDETS[1]) <> 0 and op(2, INDETS[1]) <> -1 and W <> E[1000, 1000] then if op(1, INDETS[1]) <= A then A := op(1, INDETS[1]) end if end if end do; return A end proc:
 
Coefftheta:= proc (VF) local W, i, INDETS, A; A := 10^100; for i to nops(VF) do W := op(i, VF); INDETS := indets(W, mytype); if op(0, INDETS[1]) = Theta then if op(1, INDETS[1]) <= A then A := op(1, INDETS[1]) end if end if end do; return A end proc:
 
 ZFPRODUCT := proc (m, n, k, l) return (l+2)*F[k+m, l+n]/(l+n+2)+(l*m+2*m-n-k*n)*E[k+m, l+n]/(l+n+2) end proc:
 
 ZEPRODUCT := proc (m, n, k, l) return E[k+m, l+n] end proc:
 
ZThetaPRODUCT:= proc (m, n, k, l) return Theta[k+m, l+n] end proc:
 
 ZPRODUCT := proc (m, n, VF) local zproduct, zp, i, W, INDETS, NEWVF, a0; zproduct := 0; zp := 0; NEWVF := Arrange(expand(VF+E[20000, 20000]+F[2000, 2000]))[3]; for i to nops(NEWVF) do if E[20000, 20000] = op(i, NEWVF) or F[2000, 2000] = op(i, NEWVF) then zp := 0; zproduct := zproduct+zp else W := op(i, NEWVF); INDETS := indets(W, mytype); if op(0, INDETS[1]) = E then zp := coeff(NEWVF, INDETS[1], 1)*ZEPRODUCT(m, n, op(1, INDETS[1]), op(2, INDETS[1])); zproduct := zproduct+zp elif op(0, INDETS[1]) = F then zp := coeff(NEWVF, INDETS[1], 1)*ZFPRODUCT(m, n, op(1, INDETS[1]), op(2, INDETS[1])); zproduct := zproduct+zp elif op(0, INDETS[1]) = Theta then zp := coeff(NEWVF, INDETS[1], 1)*ZThetaPRODUCT(m, n, op(1, INDETS[1]), op(2, INDETS[1])); zproduct := zproduct+zp end if end if end do; return zproduct end proc:
 
 Ytransformationmm := proc (VF, m) return F[m, m]/(2*m+2)+(1/2)*E[m, m]/m end proc:
 
Ytransformationtheta:= proc (VF, m) return Theta[m+1, m+1]/(2*m+2) end proc:

 TRUNCATE2 := proc (VF, n, s) local i, NEWVF, W, trun, truncate; trun := 0; truncate := 0; NEWVF := expand(VF+E[20000, 20000]+F[2000, 2000]); for i to nops(NEWVF) do W := op(i, NEWVF); if W <> E[20000, 20000] and W <> F[2000, 2000] and GRADE2n(W, s) <= n then trun := op(i, NEWVF) else trun := 0 end if; truncate := truncate+trun end do; return truncate end proc:
 
 WM := proc (VF, k, m) local t, UM, s, a1, r, b1, WM; 
 s := CoeffE(VF); r := CoeffF(VF); a1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), F[r, r])); b1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), E[s, s])); WM[0] := 0; UM[0] := 1; for t to m do WM[t] := (1/2)*a1*((t-1)*s+k+2)*((k-2*t+3)*(s+1)+1)*WM[t-1]/(b1*(s*t+k+2)*((t-1)*(s+1)+s))-(s*(t-2)+k)*WM[t-1]/(2*(t-1)*(s+1)+2*s)+s*(s+2)*(k-t+3)*UM[t-1]/((2*(t-1)*(s+1)+2*s)*(s*t+k+2)); UM[t] := -(1/2)*((t-1)*s+k)*((t-1)*s+k+2)*UM[t-1]/(t*(s+1)*(s*t+k+2))-(1/2)*a1*(2*t-4-k)*UM[t-1]/(b1*t)+(1/2)*a1*s*(s+2)*WM[t-1]/(b1*t*(s+1)*(s*t+k+2)) end do; return WM[m] end proc:
 
 UM := proc (VF, k, m) local t, UM, s, a1, r, b1, WM; s := CoeffE(VF); r := CoeffF(VF); a1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), F[r, r])); b1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), E[s, s])); UM[0] := 1; WM[0] := 0; for t to m do UM[t] := -(1/2)*((t-1)*s+k)*((t-1)*s+k+2)*UM[t-1]/(t*(s+1)*(s*t+k+2))-(1/2)*a1*(2*t-4-k)*UM[t-1]/(b1*t)+(1/2)*a1*s*(s+2)*WM[t-1]/(b1*t*(s+1)*(s*t+k+2)); WM[t] := (1/2)*a1*((t-1)*s+k+2)*((k-2*t+3)*(s+1)+1)*WM[t-1]/(b1*(s*t+k+2)*((t-1)*(s+1)+s))-(s*(t-2)+k)*WM[t-1]/(2*(t-1)*(s+1)+2*s)+s*(s+2)*(k-t+3)*UM[t-1]/((2*(t-1)*(s+1)+2*s)*(s*t+k+2)) end do; return UM[m] end proc:
 
 transformationf := proc (VF, k) local a1, A, m, r, s, b1; A := 0; a1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), F[r, r])); b1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), E[s, s])); for m from 0 to k+1 do A := A+WM(VF, k, m)*b1^m*E[m*(s+1)-1, m*s+k]+b1^m*UM(VF, k, m)*F[m*(s+1)-1, m*s+k] end do; return A end proc:
 
 UM0 := proc (VF, k, m) local t, UM, s, a1, r, b1, WM0, UM0; s := CoeffE(VF); r := CoeffF(VF); a1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), F[r, r])); b1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), E[s, s])); WM0[0] := 1; UM0[0] := 0; for t to m do WM0[t] := (1/2)*a1*((t-1)*s+k+2)*(k-2*t+2)*WM0[t-1]/(b1*t*(s*t+k+2))-(1/2)*(s*(t-2)+k)*WM0[t-1]/(t*(s+1))+(1/2)*s*(s+2)*(k-t+2)*UM0[t-1]/(t*(s+1)*(s*t+k+2)); UM0[t] := -((t-1)*s+k)*((t-1)*s+k+2)*UM0[t-1]/((2*t*(s+1)+2)*(s*t+k+2))-(1/2)*a1*((2*t-2-k)*(s+1)-s)*UM0[t-1]/(b1*(t*(s+1)+1))+(1/2)*a1*s*(s+2)*WM0[t-1]/(b1*(t*(s+1)+1)*(s*t+k+2)) end do; return UM0[m] end proc:
 
 WM0 := proc (VF, k, m) local t, UM, s, a1, r, b1, WM0, UM0; s := CoeffE(VF); r := CoeffF(VF); a1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), F[r, r])); b1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), E[s, s])); WM0[0] := 1; UM0[0] := 0; for t to m do UM0[t] := -((t-1)*s+k)*((t-1)*s+k+2)*UM0[t-1]/((2*t*(s+1)+2)*(s*t+k+2))-(1/2)*a1*((2*t-2-k)*(s+1)-s)*UM0[t-1]/(b1*(t*(s+1)+1))+(1/2)*a1*s*(s+2)*WM0[t-1]/(b1*(t*(s+1)+1)*(s*t+k+2)); WM0[t] := (1/2)*a1*((t-1)*s+k+2)*(k-2*t+2)*WM0[t-1]/(b1*t*(s*t+k+2))-(1/2)*(s*(t-2)+k)*WM0[t-1]/(t*(s+1))+(1/2)*s*(s+2)*(k-t+2)*UM0[t-1]/(t*(s+1)*(s*t+k+2)) end do; return WM0[m] end proc:
 
transformationvarepsilon:= proc (VF, k) local a1, A, m, r, s, b1; A := 0; a1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), F[r, r])); b1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), E[s, s])); for m from 0 to k do A := A+WM0(VF, k, m)*b1^m*E[m*(s+1), m*s+k]+b1^m*UM0(VF, k, m)*F[m*(s+1), m*s+k] end do; return A end proc:
 
 CHANGE := proc (VF, A, n) local i, l, b, f, e, d, W, k, c, INDETS, VF1, s; s := CoeffE(VF); l := 0; b := 0; f := 0; VF1 := Arrange(VF); VF1 := VF1[3]; for i to nops(VF1) do W := op(i, VF1); INDETS := indets(W, mytype); k := op(1, INDETS[1]); if op(0, INDETS[1]) = F then c := coeff(VF1, INDETS[1], 1); l := expand(l+c*add(binomial(k+1, j)*F[k-j, k-j]*A^j, j = 0 .. k)) elif op(0, INDETS[1]) = E then d := coeff(VF1, INDETS[1], 1); b := expand(b+d*add((binomial(k+1, j)-(binomial(k+1, j)-binomial(k, j))/(k+2-j))*E[k-j, k-j]*A^j, j = 0 .. k)+d*add((binomial(k+1, j)-binomial(k, j))*F[k-j, k-j]*A^j/(k+2-j), j = 1 .. k)+d*F[-1, -1]*A^(k+1)) else e := coeff(VF1, INDETS[1], 1); f := expand(f+e*add(Theta[k-j, k-j]*A^j, j = 0 .. k)) end if end do; return TRUNCATE2(l+b+f+coeff(VF, F[-1, -1])*F[-1, -1], n, s) end proc:
 
 SECOND := proc (VF, n) local r, i, j, t, W, INDETS, k, a, b, c, d, f, e, A, l, VF1, NASRIN, s, a1, b1; s := CoeffE(VF); a1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), F[s, s])); b1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), E[s, s])); VF1 := Arrange(VF); VF1 := VF1[3]; for t to n do NASRIN := VF1; for i to nops(VF1) do W := op(i, VF1); INDETS := indets(W, mytype); if op(0, INDETS[1]) = F and op(1, INDETS[1]) = s-1 and GRADE2n(W, s) = t then A := expand(-coeff(W, F[s-1, s-1])/(a1*(s+1)+b1/(s+1))); NASRIN := Arrange(CHANGE(VF1, A, n)); NASRIN := NASRIN[3] end if end do; VF1 := NASRIN end do; return VF1+coeff(VF, F[-1, 0])*F[-1, 0] end proc:
 
 progFs1 := proc (NVF, N, N1, s) local CO, l1, VF; CO := coeff(NVF, F[s-1, s-1]); l1 := ldegree(CO, {seq(mu[i], i = 1 .. 100)}); CO := mtaylor(CO, [seq(mu[i], i = 1 .. 100)], l1+1); if l1 < 0 then if s-1+s*l1 = N then print(secondary, 1); VF := SECOND(NVF, N1) else return NVF end if; return VF else return NVF end if end proc:
 
 progFs := proc (NVF, N, N1, s, a1) local l1, q1, gamma2, VF1, Y, TimeStateTrans, VF, CO; CO := simplify(coeff(NVF, F[s, s])-a1); l1 := ldegree(CO, {seq(mu[i], i = 1 .. 100)}); CO := mtaylor(CO, [seq(mu[i], i = 1 .. 100)], l1+1); if l1 < 0 then if s+s*l1 = N then gamma2 := -CO/(a1*((1/2)*s+1)); VF1 := expand(NVF+gamma2*ZPRODUCT(0, 0, NVF)); Y := expand((1/2)*gamma2*F[0, 0]); TimeStateTrans := ADJOINT2(N1, Y, VF1, s); VF := SECONDNORMALFORM(N1, TimeStateTrans) else return NVF end if; return VF else return NVF end if end proc:
 
 progFk := proc (NVF, N1, q1, s, b1, l) local k, beta2, gamma, Y, TimeStateTrans, VF; k := (q1-2*s)/(s+1); if 0 <= k and type(k, integer) then beta2 := 2*b1^(k+2)*((s+1)*(k+1)+s)*WM(NVF, k, k+2); if beta2 < 0 then gamma := expand(-mtaylor(coeff(NVF, E[q1, q1]), [seq(mu[i], i = 1 .. 100)], l+1)/beta2); if gamma < 0 then Y := expand(gamma*transformationf(NVF, k)); TimeStateTrans := ADJOINT2(N1, Y, NVF, s); VF := SECONDNORMALFORM(N1, TimeStateTrans) else return NVF end if else return NVF end if else return NVF end if; return VF end proc:
 
 progFT := proc (NVF, N, N1, l) local s, b1, a1, q1, m, A, gamma, VF1, Y, VF, TimeStateTrans; 
 s := CoeffE(expand(NVF)); 
 b1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(NVF), E[s, s]));
  a1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(NVF), F[s, s])); 
  q1 := N-l*s; 
  if l = ldegree(coeff(NVF, F[q1, q1]), {seq(mu[i], i = 1 .. 100)}) then
   m := q1-s; 
   if 0 < m then
 A := 2*a1*(s-m)/((2*m+2)*(m+2))+(-2)*b1*m*(m+2)/((2*m+2)*(m+s+2)*(m+2))+m*a1*s*(s+2)/(m*(m+s+2)*(m+2))+a1*(s+2)/(m+s+2); if A <> 0 then 
 gamma := -mtaylor(coeff(NVF, F[q1, q1]), [seq(mu[i], i = 1 .. 100)], l+1)/A; 
   VF1 := expand(NVF+gamma*ZPRODUCT(m, m, NVF)); 
   Y := expand(gamma*Ytransformationmm(VF, m));
    TimeStateTrans := expand(ADJOINT2(N1, Y, VF1, s)); 
    VF := SECONDNORMALFORM(N1, TimeStateTrans) 
    else return NVF end if else return NVF end if; return VF else return NVF end if end proc:
 
 ORBITALNF1 := proc (N1, VF) local m, gamma, VF1, Y, TimeStateTrans, NVF, s, b1, VF2, r, a1, A, alpha1, alpha2, beta3, beta4, E1, F1, alpha, beta, N, k, Y2, VF3, B, DET, c1, l, q1, CO, l1, gamma2, beta1, beta2;
 r := CoeffF(VF); print(R=r); 
 s := CoeffE(VF); print(S=s); 
 a1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), F[r, r])); 
 b1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), E[s, s]));
 c1 := subs(seq(mu[i] = 0, i = 1 .. 100), coeff(expand(VF), Theta[0, 0]));
   NVF := VF; if r < s then print(''it*isnt*at*our*case'') else for N to N1 do for l from 0 to floor(N1/s) do NVF := NVF-c1*Theta[0, 0]; NVF := progFs1(NVF, N, N1, s); NVF := progFs(NVF, N, N1, s, a1); q1 := N-l*s; k := (q1-s)/(s+1); if 0 <= k and type(k, integer) then alpha1 := simplify(2*b1^(k+1)*((s+1)*(k+1)+1)*UM0(NVF, k, k+1)); alpha2 := simplify(2*b1^(k+1)*(s+1)*(k+1)*WM0(NVF, k, k+1)); beta3 := simplify(a1*(s-k*(s+1))/((k*(s+1)+1)*(k*(s+1)+2))-b1*k*(s+1)/((k*(s+1)+1)*(k*(s+1)+s+2))+a1*s*(s+2)/((k*(s+1)+s+2)*(k*(s+1)+2))+a1*(s+2)/(k*(s+1)+s+2)); beta4 := simplify(b1*s*(s+2)/((k*(s+1)+1)*(k*(s+1)+s+2)*(k*(s+1)+2))+b1*(s-k*(s+1))/(k*(s+1)+2)+b1); E1 := simplify(-mtaylor(coeff(NVF, E[q1, q1]), [seq(mu[i], i = 1 .. 100)], l+1)); F1 := simplify(-mtaylor(coeff(NVF, F[q1, q1]), [seq(mu[i], i = 1 .. 100)], l+1)); DET := simplify(alpha1*beta4-beta3*alpha2); if DET < 0 then alpha := simplify((F1*beta4-E1*beta3)/(alpha1*beta4-beta3*alpha2)); beta := simplify((alpha1*E1-alpha2*F1)/(alpha1*beta4-beta3*alpha2)); Y2 := expand(alpha*`transformation&varepsilon;`(NVF, k)); TimeStateTrans := expand(ADJOINT2(N1, Y2, NVF, r)); VF2 := simplify(TimeStateTrans); VF3 := expand(simplify(VF2+beta*ZPRODUCT(k*(s+1), k*(s+1), VF2))); Y := expand(beta*Ytransformationmm(NVF, k*(s+1), k*(s+1))); TimeStateTrans := expand(ADJOINT2(N1, Y, VF3, s)); NVF := SECONDNORMALFORM(N1, TimeStateTrans) end if end if; NVF := progFT(NVF, N, N1, l) end do end do; VF2 := TRUNCATE2(NVF, N1, s); return Arrange(VF2)[3] end if end proc:
 
 CONVERTR := proc (VF1) local i, W, L, K, XX, rhorho, ThetaTheta, INDETS, V11, VF; XX := 0; rhorho := 0; ThetaTheta := 0; VF := expand(VF1); for i to nops(VF) do W := op(i, VF); INDETS := indets(W, mytype); L := op(1, INDETS[1]); K := op(2, INDETS[1]); if op(0, INDETS[1]) = F then XX := XX+expand(coeff(W, INDETS[1], 1)*(K-L+1)*x^(L+1)*rho^(2*K-2*L)); rhorho := rhorho-expand(coeff(W, INDETS[1], 1)*((1/2)*L+1/2)*x^L*rho^(2*K-2*L+1)); ThetaTheta := ThetaTheta elif op(0, INDETS[1]) = E then XX := XX+expand(coeff(W, INDETS[1], 1)*x^(L+1)*rho^(2*K-2*L)); rhorho := rhorho+expand((1/2)*coeff(W, INDETS[1], 1)*x^L*rho^(2*K-2*L+1)); ThetaTheta := ThetaTheta else XX := XX; rhorho := rhorho; ThetaTheta := ThetaTheta+expand(coeff(W, INDETS[1], 1)*x^L*rho^(2*K-2*L)) end if end do; return [XX, rhorho, ThetaTheta] end proc:
