(PARI) {a(n)=local(V=Vec((1+x+x^2+x^(2*n+3))^n)); V*V~}
(PARI) {a(n)=sum(k=0, n, binomial(n, k)^2*sum(j=0, k, binomial(2*k-j, j)*binomial(2*k, j)))}
(PARI) {a(n)=n!^2*polcoeff(sum(m=0, n, x^m/m!^2+x*O(x^n)) *sum(m=0, n, sum(k=0, m, binomial(2*m-k, k)*binomial(2*m, k))*x^m/m!^2), n)}