Suppose $\tau\in\Bbb C$ and $\Im(\tau)>0$. Also, let $k\in\Bbb Z_{>2}$, and $A=\Bbb Z^2\setminus\{(0,0)\}$. Then the Eisenstein series $G_{2k}(\tau)$ is given by $$G_{2k}(\tau)=\sum_{(m,n)\in A}\frac{1}{(m+n\tau)^{2k}}.$$ It is then evident that $G_{2k}(\tau+1)=G_{2k}(\tau)$ for all $\tau$. Thus we may write $G_{2k}$ as a Fourier series $$G_{2k}(\tau)=\sum_{n\ge0}g_nq^n,$$ where $q=e^{2i\pi \tau}$. Supposedly, we can find these coefficients $g_n$ by computing the integral $$g_n=\int_0^1 e^{-2i\pi n\tau}G_{2k}(\tau)d\tau=\sum_{(u,v)\in A}\int_0^1 \frac{e^{-2\pi n\tau}}{(u+v\tau)^{2k}}d\tau.$$ Apparently, these coefficients have an explicit formula: $$\begin{align} g_0&=2\zeta(2k)\\ g_n&= (-1)^k\frac{2^{2k+1}\pi^{2k}}{(2k-1)!}\sigma_{2k-1}(n),\tag{1} \end{align}$$ but I have no clue how to prove this.
I thought about writing $(u+v\tau)^{-2k}$ as a power series $(u+v\tau)^{-2k}=\sum_{l\ge0}\alpha_l\tau^l$, but that seems overly complicated. This may be simplified by writing $$G_{2k}(\tau)=\sum_{n\ge1}\frac{1}{(a_n+b_n\tau)^{2k}},$$ where $(a_n,b_n)=($$\text{A174344}$$(n),$$\text{ A274923}$$(n))$.
Take note that I'm assuming the partial sums $\sum_{n=1}^{N}(a_n+b_n\tau)^{-2k}$ converge uniformly to $G_{2k}(\tau)$ as $N\to\infty$, hence I interchanged the sum and integral. Please correct me if this is not the case.
So, what's left is the integral $$j_n(a,b)=\int_0^1\frac{e^{-2i\pi n\tau}}{(a+b\tau)^{2k}}d\tau.$$ I'm fairly certain, since $\tau$ is a complex variable, that this integral will be taken over some path in the complex plane $\gamma$ that starts at $0$ and ends at $1$. This (potentially, but I'm not sure) being the case, I do not know which path to pick.
Could I have some help proving $(1)$? Is there a better approach than the evaluation of $j_n$?
Apostol in his Modular functions and Dirichlet series in Number Theory gives a very simple proof based on the partial fraction expansion of cotangent function.
We have $$\pi\cot\pi \tau=\frac{1}{\tau}+\sum_{m\in\mathbb {Z}, m\neq 0}\left(\frac{1}{\tau+m}-\frac{1}{m}\right)\tag{1}$$ If $\tau$ has a positive imaginary part and $q=\exp(2\pi i\tau) $ then $|q|<1$ and we have $$\pi\cot\pi\tau=-\pi i\left(1+2\sum_{r=1}^{\infty}q^r\right)$$ so that we have $$\frac{1}{\tau}+\sum_{m\in\mathbb {Z}, m\neq 0}\left(\frac{1}{\tau+m}-\frac{1}{m}\right)=-\pi i\left(1+2\sum_{r=1}^{\infty} q^r\right)$$ Differentiating the above with respect to $\tau$ we get $$-\sum_{m\in\mathbb {Z}} \frac{1}{(\tau+m)^2}=-(2\pi i) \sum_{r=1}^{\infty}rq^{r-1}\frac{dq}{d\tau}=-(2\pi i) ^2\sum_{r=1}^{\infty} r\exp(2\pi i r\tau) \tag{2}$$ Differentiating this as many times as needed we can get expression for the sum $\sum(m+\tau) ^{-2k}$. The job is complete if we replace $\tau$ by $n\tau$ and sum over $n$.
Thus for example $$\sum_{m\in\mathbb {Z}} \frac{1}{(m+n\tau)^{4}}=\frac{(2\pi i) ^4}{3!}\sum_{r=1}^{\infty} r^3\exp(2\pi irn\tau) $$ and more generally we have $$\sum_{m\in\mathbb {Z}} \frac{1}{(m+n\tau)^{2k}}=\frac{(2\pi i) ^{2k}}{(2k-1)!}\sum_{r=1}^{\infty} r^{2k-1}\exp(2\pi irn\tau) =\frac{(2\pi i) ^{2k}}{(2k-1)!}\sum_{r=1}^{\infty} r^{2k-1}q^{nr}\tag{3}$$ Now we have \begin{align} G_{2k}(\tau)&=\sum_{m,n=-\infty,(m,n)\neq (0,0)}^{\infty} \frac{1}{(m+n\tau)^{2k}}\notag\\ &=\sum_{m\in\mathbb {Z}, m\neq 0}\frac{1}{m^{2k}}+\sum_{n=1}^{\infty} \sum_{m\in\mathbb {Z}} \left(\frac{1}{(m+n\tau)^{2k}}+\frac{1}{(m-n\tau)^{2k}}\right)\notag\\ &=2\zeta(2k)+2\sum_{n=1}^{\infty} \sum_{m\in\mathbb {Z}} \frac{1}{(m+n\tau)^{2k}}\notag\\ &=2\zeta(2k)+2\cdot\frac{(2\pi i) ^{2k}}{(2k-1)!}\sum_{n=1}^{\infty} \sum_{r=1}^{\infty} r^{2k-1}q^{nr}\notag\\ &=2\left(\zeta(2k)+(-1)^k\frac{(2\pi)^{2k}}{(2k-1)!}\sum_{r=1}^{\infty} \frac{r^{2k-1}q^r}{1-q^r}\right)\notag\\ &=2\left(\zeta(2k)+(-1)^k\frac{(2\pi)^{2k}}{(2k-1)!}\sum_{n=1}^{\infty} \sigma_{2k-1}(n)q^n\right)\notag \end{align}
Elliptic function theory is difficult but not too difficult if one gets hold of good books. Essentially it requires a basic knowledge of calculus and a very deep skill in algebraic manipulation. At least thats how people like Legendre, Abel, Gauss, Jacobi and Ramanujan developed the subject. Then somehow Liouville and Weierstrass arrived on the scene and ditched algebraic manipulation altogether and exclusively used complex analysis making the topic widely inaccessible.