Flatten[Table[2^(2*n+1) + 1 + 2^(n+1) {-1, 1}, {n, 0, 40}]] (* J. Mulder (jasper.mulder(AT)planet.nl), Jan 28 2010 *)
PROG
(PARI) my(x='x+O('x^40)); Vec((1+4*x-6*x^2-16*x^3+20*x^4)/((1-x)*(1- 2*x^2)*(1-4*x^2))) \\ G. C. Greubel, Jun 01 2019
(Magma) R<x>:=PowerSeriesRing(Integers(), 40); Coefficients(R!( (1+4*x-6*x^2-16*x^3+20*x^4)/((1-x)*(1- 2*x^2)*(1-4*x^2)) )); // G. C. Greubel, Jun 01 2019
(SageMath) ((1+4*x-6*x^2-16*x^3+20*x^4)/((1-x)*(1- 2*x^2)*(1-4*x^2))).series(x, 40).coefficients(x, sparse=False) # G. C. Greubel, Jun 01 2019