> restart;with(linalg):
M1      :=   1:                #  No. of non-zero real eigenvalues
M2      :=   2:                #  No. of complex conjugate eigenvalues
omg    :=    1:                #  omega_c of linearized system 
Ord     :=   2:                 
N       :=   2 + M1 + M2*2:     # In my case N=7

s:=sqrt(1-c^2):#c:=0.5:b:=1:             
func[1] := (1/2)*(c+s)*x[1]+s*x[3]+(s-c)*x[1]*x[2];
func[2] := (1/2)*(c-s)*x[1]+s*x[3]+(c+s)*x[1]*x[2];
func[3] := -b*x[2]-x[1]^2;  
Solve := proc(k,Y,X)
   local XX:
   if  k=1 then
       if  coeff(Y,X)<>0 then
           XX := solve(Y,X):
       else
           XX := 0:
       fi:
   else
       XX := Y:
   fi:
   XX := rationalize(XX):
   XX := expand(XX):
   X  := XX:
end:

solution := proc(n)
   local  i, j, k, A, B, C, D, E, f, myA, myb, temp:
   global Ord, func, F1, F2, FF, X, Dr, Dphi, DDr, DDphi:

# Solve for solution  x_n1.

   for i from 1 to N do
       f[i] := 0:
       for j from 1 to n do
           f[i] := f[i]
           + combine(subs(eps=0,diff(FF[i,j],eps$(n+1-j))/(n+1-j)!),trig):
       od:
       f[i] := f[i]+ combine(subs(eps=0,FF[i,n+1]),trig):
   od:
   for i from 1 to N do
       for j from 1 to n do
           f[i] := f[i]-(r^(j+1)*Dr[j]*diff(X[i,n-j],r)
                       + r^j*Dphi[j]* combine(diff(X[i,n-j],theta),trig)):
       od:
   od:
   X[1,n] := A[0]:
   for i from 2 to n+1 do
       X[1,n] := X[1,n] + A[i]*cos(i*theta) + B[i]*sin(i*theta):
   od:
   X[1,n] := r^(n+1)*X[1,n]:
   F1 := diff(X[1,n],theta$2) + X[1,n] - diff(f[1],theta) - f[2]:

# Solve for the coefficients of the normal forms, D_n a  and  D_n phi.

   Solve(1, coeff(F1, sin(theta)), Dr[n]):
   Solve(1, coeff(F1, cos(theta)), Dphi[n]):

   for i from 2 to n+1 do
       Solve(1, coeff(F1, cos(i*theta)), A[i]):
       Solve(1, coeff(F1, sin(i*theta)), B[i]):
   od:
   Solve(1, F1, A[0]):

# Find solution x_n2.

   X[1,n] := X[1,n]:
   X[2,n] := diff(X[1,n], theta) - f[1]:

# Solve for solutions  x_np, p = 3, 4, ... M1+2.

   if M1 <> 0 then
      for i from 3 to M1+2 do
          A := 'A':
          B := 'B':
          X[i,n] := A[0]:
          for j from 1 to n+1 do
              X[i,n] := X[i,n] + A[j]*cos(j*theta) + B[j]*sin(j*theta):
          od:
          F1 := combine(diff(X[i,n], theta) - f[i], trig):
          for j from 1 to n+1 do
              for k from 1 to 2 do
                  if   k=1 then
                       temp := coeff(F1, cos(j*theta)):
                  elif k=2 then
                       temp := coeff(F1, sin(j*theta)):
                  fi:
                  C[k,1] := coeff(temp, A[j]):
                  C[k,2] := coeff(temp, B[j]):
                  C[k,3] := C[k,1]*A[j] + C[k,2]*B[j] - temp:
              od:
              myA  := array([[C[1,1], C[1,2]], [C[2,1], C[2,2]]]):
              myb  := array( [C[1,3], C[2,3]]):
              myb  := linsolve(myA, myb):
              Solve(2, myb[1], A[j]):
              Solve(2, myb[2], B[j]):
          od:
          F1 := combine(F1, trig):
          if F1 <> 0 then
             Solve(1, F1, A[0]):
          fi:
          X[i,n] := X[i,n]:
      od:
   fi:

# Solve for solutions  x_nq  and  x_n(q+1), q = M1+3, M1+5, ... N-1.

   if M2 <> 0 then
      for i from M1+3 by 2 to N do
          A := 'A':
          B := 'B':
          X[i,n]   := A[0]:
          X[i+1,n] := D[0]:
          for j from 1 to n+1 do
              X[i,n]   := X[i,n]   + A[j]*cos(j*theta) + B[j]*sin(j*theta):
              X[i+1,n] := X[i+1,n] + D[j]*cos(j*theta) + E[j]*sin(j*theta):
          od:
          F1:= combine(diff(X[i,n],theta) - f[i], trig):
          F2:= combine(diff(X[i+1,n],theta) - f[i+1], trig):
          for j from 1 to n+1 do
              for k from 1 to 4 do
                  if   k=1 then
                       temp := coeff(F1, cos(j*theta)):
                  elif k=2 then
                       temp := coeff(F1, sin(j*theta)):
                  elif k=3 then
                       temp := coeff(F2, cos(j*theta)):
                  elif k=4 then
                       temp := coeff(F2, sin(j*theta)):
                  fi:
                  C[k,1] := coeff(temp, A[j]):
                  C[k,2] := coeff(temp, B[j]):
                  C[k,3] := coeff(temp, D[j]):
                 # C[k,4] := coeff(temp, E[j]):
                 # C[k,5] := C[k,1]*A[j] + C[k,2]*B[j]
                           +C[k,3]*D[j] + C[k,4]*E[j] - temp:
              od:
              myA  := array([[C[1,1], C[1,2], C[1,3], C[1,4]],
                             [C[2,1], C[2,2], C[2,3], C[2,4]],
                             [C[3,1], C[3,2], C[3,3], C[3,4]],
                             [C[4,1], C[4,2], C[4,3], C[4,4]]]):
              myb  := array( [C[1,5], C[2,5], C[3,5], C[4,5]] ):
              myb  := linsolve(myA, myb):
              Solve(2, myb[1], A[j]):
              Solve(2, myb[2], B[j]):
              Solve(2, myb[3], D[j]):
              Solve(2, myb[4], E[j]):
          od:
          F1 := combine(F1, trig):
          if F1 <> 0 then
             Solve(1, F1, A[0]):
          fi:
          F2 := combine(F2, trig):
          if F2 <> 0 then
             Solve(1, F2, D[0]):
          fi:
      od:
   fi:
end:

func[1] := func[1]-omg*x[2]:
func[2] := func[2]+omg*x[1]:
for i from 1 to N do
    func[i] := func[i]/omg:
od:
for i from 1 to N do
    for j from 1 to N do
        func[i] := subs(x[j]=eps*x[j],func[i]):
    od:
    func[i] := rationalize(func[i]):
    func[i] := expand(func[i]):
    for j from 1 to Ord+1 do
        FF[i,j] := subs(eps=0,diff(func[i],eps$j)/j!):
    od:
od:

# Main program; write the formal ordered perturbation solution  x_i.

for i from 1 to N do
    x[i] := 0:
    for j from 0 to Ord do
        x[i] := x[i] + X[i,j]*eps^j:
    od:
od:

# Zero order solution  x_0.

X[1,0] :=   r*cos(theta);
X[2,0] := - r*sin(theta);
if N > 2 then
   for i from 3 to N do
       X[i,0] := 0:
   od:
fi:

# Call procedure "solution" to obtain the asymptotic solutions and the
# coefficients of the normal forms for the i-step perturbation equations.

for i from 1 to Ord do
    solution(i):
    print(` Order = `, i):
    for j from 1 to i do
        DDr[j]   := simplify(subs(r=1,omg*Dr[j])):
        DDphi[j] := simplify(subs(r=1,omg*Dphi[j])):
    od:
od:
for i from 1 to N do
    x[i] := x[i]:
od: 
print(` a1= `,  DDr[2], ` a2 = `, DDr[4]);  
print(` t1= `,  Dphi[2], ` t2 = `, Dphi[4]);  
save i, Ord, DDr, DDphi, output:       # Store results in the file "output"


          /            (1/2)\                (1/2)     
        1 |    /     2\     |        /     2\          
        - \c + \1 - c /     / x[1] + \1 - c /      x[3]
        2                                              

             /        (1/2)    \          
             |/     2\         |          
           + \\1 - c /      - c/ x[1] x[2]
          /            (1/2)\                (1/2)     
        1 |    /     2\     |        /     2\          
        - \c - \1 - c /     / x[1] + \1 - c /      x[3]
        2                                              

             /            (1/2)\          
             |    /     2\     |          
           + \c + \1 - c /     / x[1] x[2]
                                      2
                        -b x[2] - x[1] 
                          r cos(theta)
                         -r sin(theta)
Error, (in combine) too many levels of recursion
Error, too many levels of recursion
                  a1= , DDr[2],  a2 = , DDr[4]
                 t1= , Dphi[2],  t2 = , Dphi[4]
Warning, unassigned variable `DDr` in save statement
Warning, unassigned variable `DDphi` in save statement