pt = {{1}}; Table[rhs = CoefficientList[(k^2 + x)^k, x]; qt = Join[pt, {vars = Array[Subscript[a, #] &, k + 1]}]; b = MatrixPower[PadRight[qt], k] ; {out} = vars /. Solve[Thread[Reverse[Last[b]] == Reverse[rhs]], vars]; pt = Join[pt, {out}]; out, {k, 20}] not the fastest of codes : the first 20 rows take 43 seconds on my old machine. W.