I'm trying to understand the method exposed in Dorf's "Modern Control Systems", for approximating a high order system with a lower order one. The method goes as follows (using the exact same notation): Let $G_H(s)$ and $G_L(s)$, be the transfer functions of the original system and the approximate one, respectively. Take $G_H(s)/G_L(s)$, and call $M(s)$ and $\Delta(s)$ the numerator and denominator polynomials of the resulting quotient, respectively. Define
$$M_{2q} = \sum_{k=0}^{2q} \frac{(-1)^{k+q}M^{(k)}(0)M^{(2q-k)}(0)}{k!(2q-k)!}$$
and
$$\Delta_{2q} = \sum_{k=0}^{2q} \frac{(-1)^{k+q}\Delta^{(k)}(0)\Delta^{(2q-k)}(0)}{k!(2q-k)!}$$
and make $M_{2q}=\Delta_{2q}$, for $q=0,1,2\dots$, as needed to solve for all the unknowns in $G_L(s)$.
The questions is, why this works, how'd they derived it. I tried to undersand it by looking at the book's references and found a paper from the 60s by E. Davison but couldn't access it (paywall). Any help to see how it works or another source to delve in.