<s>[INST] Maclaurin expansion of $\exp(x)/(1-\exp(x))^2$ I am trying to expand following function $f(x)=\frac{e^x}{(1-e^x)^2}$ using the Maclaurin series.
According to wolframalpha it should be:
$\frac{1}{x^2}-\frac{1}{12}+\frac{x^2}{240}-\frac{x^4}{6048}+O(x^6)$
However, I am struggling to arrive at this result following a standard formula for Maclaurin expansion:
$f(x) = f(0)+f'(0)x+\frac{f''(0)}{2!}x^2+...$
In our case $\lim(x \to 0)\frac{e^x}{ (1 - e^x)^2 } = \infty$.
0 is a singular point for the $f(x)$ and zero appears in the denominator of every n'th order derivative as well. 
So, how can I expand this function?
 [/INST] Bernoulli numbers $B_n$ are defined in a way that:
$$\frac{x}{e^x-1}= \sum_{n=0}^{\infty} \frac{B_n x^n}{n!}\tag{1}\label{eq1}$$
Divide \eqref{eq1} by $x$:
$$\frac{1}{e^x-1}= \sum_{n=0}^{\infty} \frac{B_n x^{n-1}}{n!}\tag{2}\label{eq2}$$
Calculate the first derivative of both sides of the \eqref{eq2}:
$$-\frac{e^x}{(1-e^x)^2}= \sum_{n=0}^{\infty} B_n \frac{(n-1)x^{n-2}}{n!}\tag{3}\label{eq3}$$
Here are first few Bernoulli numbers: $B_0 = 1, B_1 = -\frac{1}{2}, B_2 = \frac{1}{6}, B_3 = 0, B_4 = -\frac{1}{30}, B_5 = 0, B_6 = \frac{1}{42}, B_7=0$
Substitute these in \eqref{eq3}:
$$\frac{e^x}{(1-e^x)^2}= -B_0x^{-2}-B_2\frac{1}{2!}x^0-B_4\frac{3}{4!}x^2-B_6\frac{5}{6!}x^4+O(x^6)=\frac{1}{x^2}-\frac{1}{12}+\frac{x^2}{240}-\frac{x^4}{6048}+O(x^6)$$
</s>