Solution (source code)

= Solution

For a nonzero lattice point $\omega$ and $|z|<|\omega|$,
$$
\frac1{(z-\omega)^2}-\frac1{\omega^2}
=\sum_{r=1}^{\infty}(r+1)\frac{z^r}{\omega^{r+2}}.
$$
Pairing $\omega$ with $-\omega$ cancels odd powers. There is no constant term, so $a_0=0$, and the <laurent coefficients of the Weierstrass elliptic function> are
$$
\boxed{a_{2k}=(2k+1)
\sum_{(m,n)\ne(0,0)}
\frac1{(m\omega_1+n\omega_2)^{2k+2}},
\qquad k\geq1.}
$$