There is an established result that the Fourier expansion of a particular ratio of Jacobi elliptic theta functions: $$\frac{\theta_1(x+y)\theta_1'(0)}{\theta_1(x)\theta_1(y)} = \cot(x)+\cot(y)+4\sum_{n=1}^{\infty}\sum_{m=1}^{\infty}q^{2mn}\sin(2my+2nx)$$ where $|\Im(x)| <|\Im(\pi\tau)|$ and $|\Im(y)| <|\Im(\pi\tau)|$ c.f. Whittaker and Watson 5th edition page 515 Example 21.13.
One can evaluate the integrals in $x$ and $y$ for the Fourier coefficient by considering contour integrals about the parallelogram in the complex plane with vertices at $\left\{\pm \dfrac{\pi}{2}, \pm\dfrac{\pi}{2}+\pi\tau\right\}$ and arrive at the result by computing the sum.
I am interested in evaluating the Fourier expansion for the following ratio of Jacobi elliptic theta functions: $$\frac{\theta_4(x+y)\theta_1'(0)}{\theta_4(x)\theta_4(y)}=\sum_{m,n}e^{2inx+2imy}A_{m,n}$$ where $$A_{m,n} = \frac{1}{\pi^2}\int_{-\pi/2}^{\pi/2}\mathrm{d}ye^{-2imy}\int_{-\pi/2}^{\pi/2}\mathrm{d}x e^{-2inx } \frac{\theta_4(x+y)\theta_1'(0)}{\theta_4(x)\theta_4(y)}$$ I can evaluate the integral with respect to $x$ by considering the parallelogram contour as before, giving $$A_{m,n} = \frac{2i}{\pi}\int_{-\pi/2}^{\pi/2}\mathrm{d}ye^{-2imy} \frac{q^{-n}e^{-iy}}{1-q^{-2n}e^{-2iy}} \frac{\theta_1(y)}{\theta_4(y)}$$ However, if I try to compute this integral using the same method I instead generate the following recurrence relation:
$$ A_{m,n} - q^{-2m}A_{m,n+1} = -4\frac{\theta_4(0)}{\theta_1'(0)} \frac{q^{-m}q^{-(n+\frac{1}{2})}}{1-q^{-(2n+1)}}$$
This seems simple, but I do not have an explicit boundary condition to solve this difference relation. I can see that the original function is symmetric under interchange of $x\leftrightarrow y$ and $x\rightarrow -x, y\rightarrow-y$ simultaneously, which may imply $A_{m,n}=A_{n,m}$ and $A_{m,n}=A_{-m,-n}$ respectively, but I can't seem to reconcile this with the difference equation.
Any help would be greatly appreciated!