You can find a fully worked-out derivation of $P(z)$ on an elliptic curve in Appendix A of Feynman Graph Integrals and Almost Modular Forms 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 physics.stackexchange 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.$$