(PARI) /* By definition: */
{a(n)=local(A=1); A=sum(m=0, n, x^m/(1-x^2)^(2*m+1) * sum(k=0, m, binomial(m, k)^2 * x^k) * sum(k=0, m, binomial(m, k)^2 * x^(2*k)) +x*O(x^n)); polcoeff(A, n)}
for(n=0, 30, print1(a(n), ", "))
(PARI) /* By a binomial identity: */
{a(n)=polcoeff(sum(m=0, n, x^m*sum(k=0, m, binomial(m, k)^2 * x^(m-k) * sum(j=0, k, binomial(k, j)^2 * x^j )+x*O(x^n))), n)}
for(n=0, 30, print1(a(n), ", "))
(PARI) /* From a formula for a(n): */
{a(n)=sum(k=0, n\2, binomial(2*k, k)*binomial(n-k, k)^2)}
for(n=0, 30, print1(a(n), ", "))