F = 1+3*y^2
G = y^2*(1-y^2)^2

a={a1,a2,a3,a4,a5,a6,a7,a8,a9,a10,a11,a12,a13,a14,a15,a16,a17,a18,a19,a20,a21,a22,a23,a24,a25,a26,a27,a28,a29,a30,a31,a32,a33,a34,a35,a36,a37};

For[m=1,m<38,m++,
  n=6*m;
  f = (F^(n/2))*(G^0)+Sum[(a[[j]])*F^(n/2-3*j)*G^(j),{j,1,n/6}];
  ta = Table[a[[i]],{i,1,(n-6)/6}];
  (* Find a_1,...,a_{m-1} *)
  sol = Solve[Table[Coefficient[f,y^(2*i)]==0,{i,1,(n-6)/6}],ta];
  (* Substitute a_1,...,a_{m-1} for f*)
  fsol = f /. sol[[1]];
  sol0 = Solve[Coefficient[fsol,y^(2*m)]==alpha,{a[[m]]}];
  fsol0 = fsol /. sol0[[1]];
  file = ToString[m]<>".txt";
  (* Output a weight enumerator *)
  Export[file,"W_{4,"<>ToString[n]<>"}="<>ToString[InputForm[Collect[Expand[fsol0],y]]]]
]
