$\zeta(0)$ and the cotangent function

This is not a completely satisfactory answer. I would like a simpler one. Nevertheless still probably a good exercise in Complex variables. I will only sketch it.

What we want to show is equivalent to $$\zeta(2n)=-\frac{1}{2\pi i}\int_{C_r}\frac{\pi z \cot(\pi z)}{2z^{2n+1}}\,dz,\qquad n\ge 0,\quad n\in{\bf Z}.\tag{1}$$ In fact this will be true for all $n\in{\bf Z}$. For $z=ix$ with $x>0$ we have $$\cot(\pi z)=\cot(\pi i x)=-i-i\frac{2}{e^{2\pi x}-1}$$ it is convenient to write (1) as $$\zeta(2n)=-\frac{1}{2 i}\int_{C_r}\frac{ \cot(\pi z)+i}{2z^{2n}}\,dz,\qquad n\in{\bf Z}.\tag{2}$$

Consider now the region $\Omega$ equal to ${\bf C}$ with a cut along the positive imaginary axis. Let $\log z$ denote the determination of the logarithm in $\Omega$ with $-\frac{3\pi}{2}<\arg(z)<\frac{\pi}{2}$, and let $C'_r$ be the path of integration that start at $i\infty$ to $ir$ (left border of the imaginary positive axis), then follows the circumference $C_r$ from $ir$ to $ir$ and then go from $ir$ to $i\infty$ (right border of the imaginary positive axis). It is easy to show that (2) is equivalent to (3) $$\zeta(2n)=-\frac{1}{2 i}\int_{C'_r}\frac{ \cot(\pi z)+i}{2z^{2n}}\,dz,\qquad n\in{\bf Z}. \tag{3}$$

The integral defines an entire function $$f(s)=-\frac{1}{4 i}\int_{C'_r}\bigl(\cot(\pi z)+i\bigr)e^{-s\log z}\,dz.\tag{4}$$

Expanding the integral we get $$f(s)= \frac{i}{2}\bigl(e^{-\pi i s/2}-e^{3\pi i s/2}\bigr)\int_r^{\infty}\frac{x^{-s}}{e^{2\pi x}-1}\,dx-\frac{1}{4 i}\int_{C_r}\bigl(\cot(\pi z)+i\bigr)e^{-s\log z}\,dz.$$ When $r\to0$ the last integral tends to $0$ if we have $\Re(s)=\sigma<0$. So in this case we get $$f(s)= \frac{i}{2}\bigl(e^{-\pi i s/2}-e^{3\pi i s/2}\bigr)\int_0^{\infty}\frac{x^{-s}}{e^{2\pi x}-1}\,dx,\qquad \sigma=\Re(s)<0.$$ from which we get (5)

$$f(s)= e^{\pi i s/2}\sin(\pi s)(2\pi)^{s-1}\int_0^{\infty}\frac{x^{-s}}{e^{x}-1}\,dx,\qquad \sigma<0.\tag{5}$$ Applying Titchmarsh (2.4.1) $$\zeta(s)=\frac{1}{\Gamma(s)}\int_0^\infty\frac{x^{s-1}}{e^x-1}\,dx,\qquad \sigma>1$$ and the functional equation, yields that for $\sigma<0$ we have $$f(s)=e^{\pi i s/2}\cos(\pi s/2)\zeta(s)\tag{6}$$ Therefore for all $s$ we have $$e^{\pi i s/2}\cos(\pi s/2)\zeta(s)=-\frac{1}{4 i}\int_{C'_r}\bigl(\cot(\pi z)+i\bigr)e^{-s\log z}\,dz.\tag{7}$$ Since $f(2n)=\zeta(2n)$ for all $n\in{\bf Z}$ we have proved that the coefficient of $z^{2n}$ in the Laurent series for $-\frac{\pi z}{2}\cot(\pi z)$ is equal to $\zeta(2n)$ for all $n\in{\bf Z}$.

Here is an explanation based on the Euler-Maclaurin summation formula. (Or rather, since we'll only ever need two terms of the Euler-Maclaurin summation, it's really more or less just "the trapezoid rule".) I think it's a good explanation because it sticks to the general structure of the argument outlined in the question, and its "18th-century-friendly" spirit.

First, let's review how to apply Euler-Maclaurin to $\zeta$. We have $$ \begin{align} \zeta(s) &= \sum_{n=1}^N \frac{1}{n^s} + \sum_{n=N+1}^{\infty} \frac{1}{n^s} \\ &= \sum_{n=1}^N \frac{1}{n^s} + \int_N^{\infty} \frac{1}{x^s} \, dx - \frac12 \frac{1}{N^s} + \mathrm{Error} \\ &= \sum_{n=1}^N \frac{1}{n^s} + \frac{1}{s-1} \frac{1}{N^{s-1}} - \frac12 \frac{1}{N^s} + \mathrm{Error} \tag{1}. \\ \end{align} $$ We can say more about the error term later - for now, all we need to know is that $\mathrm{Error} = O ( \frac{1}{N^{\operatorname{Re}(s)+1}} )$ as $N \to \infty$, for each $s$. The key thing to note is that, even though the computation initially required $\operatorname{Re}(s)>1$ in order to be valid, in fact the "equation" $$ \zeta(s) = \sum_{n=1}^N \frac{1}{n^s} + \frac{1}{s-1} \frac{1}{N^{s-1}} - \frac12 \frac{1}{N^s} + O \left( \frac{1}{N^{\operatorname{Re}(s)+1}} \right) \tag{2} $$ has a unique constant solution $\zeta(s)$, for each $s$ with $\operatorname{Re}(s) > -1$ and $s \ne 1$. This is one way to define $\zeta(s)$ for all such $s$. In particular, we can plug in $0$ and immediately find that $\zeta(0) = -\frac12$.

What about the function $f(z) = -\frac{\pi z}{2} \cot(\pi z)$? We can also apply Euler-Maclaurin to $f$ in the same way: $$ \begin{align} f(z) &= -\frac12 + \sum_{n=1}^N \frac{z^2}{n^2-z^2} + \sum_{n=N+1}^{\infty} \frac{z^2}{n^2-z^2} \\ &= -\frac12 + \sum_{n=1}^N \frac{z^2}{n^2-z^2} + \int_N^{\infty} \frac{z^2}{x^2-z^2} \, dx - \frac12 \frac{z^2}{N^2-z^2} + O \left( \frac{1}{N^3} \right) \tag{3} \\ \end{align} $$ as $N \to \infty$, for each $z$. The right-hand side of (3), without the error term, is what we'll call $f_N(z)$; we can simplify it to $$ f_N(z) = \sum_{n=1}^N \frac{z^2}{n^2-z^2} + z \cdot \frac12 \left( \log\left(1+\frac{z}{N}\right) - \log\left(1-\frac{z}{N}\right) \right) - \frac12 \frac{N^2}{N^2-z^2} \tag{4}. $$ I should emphasize that (4) is just a somewhat more elaborate variant of $-\frac12 + \sum_{n=1}^{\infty} \frac{z^2}{n^2-z^2}$, much like how (2) is just a somewhat more elaborate variant of $\sum_{n=1}^{\infty} \frac{1}{n^s}$.

Now when we expand (4) into power series, we recognize the expression from (2) as the coefficients, and we recognize that we can make the sum run from $k=0$ to $\infty$, not just $k=1$ to $\infty$: $$ \begin{align} f_N(z) &= \sum_{n=1}^N \sum_{k=1}^{\infty} \frac{z^{2k}}{n^{2k}} + \sum_{k=1}^{\infty} \frac{1}{2k-1} \frac{1}{N^{2k-1}} z^{2k} - \frac12 \sum_{k=0}^{\infty} \frac{z^{2k}}{N^{2k}} \\ &= \sum_{k=0}^{\infty} \left( \sum_{n=1}^N \frac{1}{n^{2k}} + \frac{1}{2k-1} \frac{1}{N^{2k-1}} - \frac12 \frac{1}{N^{2k}} \right) z^{2k}. \tag{5} \\ \end{align} $$ This is exactly what we want: take the limit as $N \to \infty$ to conclude that $f(z) = \sum_{k=0}^{\infty} \zeta(2k) z^{2k}$.

(To be rigorous in the final step, we'd need to be more precise about the error term in (1). The actual bound we get from Euler-Maclaurin is $ \lvert \mathrm{Error} \rvert \le \mathrm{constant} \cdot \int_N^{\infty} \left\lvert \frac{s(s+1)}{x^{s+2}} \right\rvert \, dx $ for all $N$ and all $s$ such that $\operatorname{Re}(s) > -1$ and $s \ne 1$. This lets us control the size of the difference $\sum_{k=0}^{\infty} \zeta(2k) z^{2k} - f_N(z)$.)

This proves that all the even-power coefficients of the power series of $-\frac{\pi z}{2} \cot(\pi z)$ are the corresponding values of $\zeta$, including $\zeta(0)$ as the constant term.