Your solution is basically correct, there are just some special cases that need to be handled. It is not always true that $\wp(z) - \wp(a_k)$ has a simple zero in $a_k$. If $a_k$ is a zero of $\wp'$, then $\wp(z)-\wp(a_k)$ has a double zero in $a_k$, and similar for the poles $b_k$ of course. If none of the zeros or poles of $f$ coincides with a zero of $\wp'$, then the construction goes through without any problems, and you have
$$f(z) = C\prod_{k=1}^n \frac{\wp(z)-\wp(a_k)}{\wp(z)-\wp(b_k)}$$
where the $a_k$ resp. $b_k$ are the zeros resp. poles of $f$ in the fundamental parallelogram for an even elliptic $f$ that has neither a pole nor a zero in $0$.
What if one (or more) of the $a_k$ resp. $b_k$ is a zero of $\wp'$?
In this question, we saw that $\wp'$ has the three distinct zeros
$$\rho_1 = \frac{\omega_1}{2},\; \rho_2 = \frac{\omega_1+\omega_2}{2},\; \rho_3 = \frac{\omega_2}{2},$$
and since the order of $\wp'$ is three, these are all simple zeros, and $\wp'$ has no other zeros (modulo the lattice $\Omega = \langle \omega_1,\omega_2\rangle$). The argument used the oddness and periodicity of $\wp'$, but of course $f'$ is also an odd elliptic function for the lattice $\Omega$, so the same argument yields
$$-f'(\rho_1) = f'(-\rho_1) = f'(-\rho_1+\omega_1) = f'(\rho_1),$$
hence $f'(\rho_1) = 0$, if $f$ doesn't have a pole in $\rho_1$,and similar for $\rho_2$, $\rho_3$. Thus if any of the $\rho_i$ is a zero of $f$, it is a zero of even order (if the order is greater than $2$, divide out one factor $\wp(z)-\wp(\rho_i)$ and repeat the argument), and you include the factor $\wp(z)-\wp(\rho_i)$ only half as often in the product. If one of the $\rho_i$ is a pole of $f$, the same argument for $1/f$ shows that the pole must have even order, and then you include the factor $\dfrac{1}{\wp(z)-\wp(\rho_i)}$ only half as often as the multiplicity of the pole would indicate.
Now, if $f$ has a zero or a pole in any of the $\rho_i$, it may happen that the halving of the factors $\wp(z) - \wp(\rho_i)$ produces a different number of factors in the numerator than in the denominator. But that means that $f$ then must have either a zero or a pole in $0$, so this cannot happen for an even elliptic function that has neither a zero nor a pole in $0$ (sorry, I'd rather have a more elegant proof of that fact, but this will have to do for now).
Best Answer
$$\wp_{\tau}(z) = \frac1{z^2} + \sum_{(n,m) \ne (0,0)}\frac1{(z+n\tau+m)^2}-\frac1{(n\tau+m)^2}$$ is $1$ periodic and analytic on $\Im(z) \in (0,|\Im(\tau)|)$ so it has a Fourier series valid for $\Im(z) \in (0,|\Im(\tau)|)$ which is a Laurent series for $\wp_{\tau}(\frac{\log u}{2i\pi})$ valid for $|e^{2i \pi u}| \in (1,e^{|\Im(\tau)|})$
The Fourier coefficients are found from the pole expansion of inverse of trigonometric functions. For $\Im(z) > 0$ $$\sum_{m}\frac{1}{(z+m)^2}= \frac{d}{dz} \frac1{1-e^{2i \pi z}}= \sum_{k \ge 0} (2i\pi k) e^{2 i \pi kz}$$ For $\Im(z) < 0$ $$\sum_{m}\frac{1}{(z+m)^2}= \frac{d}{dz} \frac1{1-e^{2i \pi z}}= \sum_{k \ge 1} (-2i\pi k) e^{-2 i \pi kz}$$