You can find a fully worked-out derivation of $P(z)$ on an elliptic curve in Appendix A of <A HREF="http://arxiv.org/abs/1112.4015">Feynman Graph Integrals and Almost Modular Forms</A> by S. Li (2011). Basically, the result you quote appears upon evaluation of $\sum_{n,m=-\infty}^{\infty}(z-m-n)^{-2}$. This gives the Weierstrass elliptic function $\wp$ and Eisenstein series $E_2$ in the second equation of your posting.

You ask "What am I doing wrong?"

There is some confusing terminology here (mentioned also in response to your earlier <A HREF="http://physics.stackexchange.com/questions/110655/how-to-compute-the-propogator-for-chern-simons-on-a-torus">physics.stackexchange</A> posting). The object $P(z)$ is called the "propagator" in the conformal field theory context, but it is actually the second derivative of your propagator. More specifically, the function $P(z)$ is not the inverse of the Laplacian $\Delta_z$, which is the function you are computing in your posting, but instead

$$P(z)=\int_0^\infty \frac{d^2}{dz^2}e^{-t\Delta_z}\,dt.$$