Finding the coefficients $p_0,p_1,p_2,q_1,q_2,q_3$ of Padé approximation Determine the Padé approximation of degree $5$ with $ n =2 $ and $ m= 3$ for $f(x) = e^{-x}$.

Suppose $r$ is a rational function of degree $N$.$$
r(x) = \frac{p(x)}{q(x)} = \frac{p_0 +p_1x + \cdots + p_nx^n}{q_0 +q_1x+ \cdots + q_mx^m}.
$$
The Padé approximation technique is$$
f(x)-r(x) = f(x)-\frac{p(x)}{q(x)} = \frac{f(x)q(x)-p(x)}{q(x)} = \frac{\left(\sum\limits_{i=0}^{\infty} a_ix^i \right)\left(\sum\limits^m_{i=0}q_{i}x^i \right)-\sum\limits^n_{i=0} p_ix^i}{q(x)}.
$$
To find the Padé approximation we need to choose $p_0,p_1,p_2,q_1$ and $q_2,q_3$. The Maclaurin expansion of $e^{-x}$ is $\sum\limits_{i = 0}^{\infty} \dfrac{(-1)^i}{i!} x^i$. Since $q_0 \neq 0$, thus $q_0= 1$ if $x= 0$, and $p_0 =1$. How does one find $p_1,p_2,q_1,q_2,q_3$ ao that$$
r_{2,3}(x) = \frac{p_0+p_1x+p_2x^2}{q_0+q_1x+q_2x^2+q_{3}x^3}?
$$
Does one have to transform this into a linear equation to solve?$$
\begin{bmatrix}
a_0 &  &  & &\\
a_1 & a_0 & & &\\
a_2 & a_1 & a_0&& \\
a_3 & a_2 & a_1& a_0&\\
a_4 & a_3 & a_2 & a_1&a_0\\
\end{bmatrix} \begin{bmatrix}
q_1\\
q_2\\
q_3\\
0\\
0
\end{bmatrix} - \begin{bmatrix}
p_1\\
p_2\\
0\\
0\\
0
\end{bmatrix} = -\begin{bmatrix}
a_1\\
a_2\\
a_3\\
a_4\\
a_5
\end{bmatrix}
$$
How do I solve this linear system so I may obtain the solution? Can someone demonstrate a concise method of Gaussian elimination to find $r_{2,3}$ and that can also work for $r_{1,4}$ or any $r_{n,m}$ for Padé approximation?
 A: $\newenvironment{gmatrix}{\left\lgroup\begin{matrix}}{\end{matrix}\right\rgroup}$In general, to find the Padé approximation $r(x) = \dfrac{\sum\limits_{k = 0}^n p_k x^k}{\sum\limits_{k = 0}^m q_k x^k}$ for $f(x) = \sum\limits_{k = 1}^\infty a_k x^k$, it reduces to find the solution to$$
\begin{gmatrix}
a_0 &&&\\
a_1 & a_0 &&\\
\vdots & \ddots & \ddots &\\
a_{m + n} & \cdots & a_1 & a_0
\end{gmatrix} \begin{gmatrix}
q_0 \\ q_1 \\ \vdots \\ q_{m + n}
\end{gmatrix} = \begin{gmatrix}
p_0 \\ p_1 \\ \vdots \\ p_{m + n}
\end{gmatrix}, \tag{1}
$$
where $p_{n + 1} = \cdots = p_{m + n} = 0$, $q_{m + 1} = \cdots = q_{m + n} = 0$. In order to solve this system of equations, usually it is assumed that $q_0 = 1$ and (1) becomes$$
\begin{gmatrix}
a_0 &&\\
\vdots & \ddots &\\
a_{m + n - 1} & \cdots & a_0
\end{gmatrix} \begin{gmatrix}
q_1 \\ \vdots \\ q_{m + n}
\end{gmatrix} - \begin{gmatrix}
p_1 \\ \vdots \\ p_{m + n}
\end{gmatrix} = \begin{gmatrix}
a_1 \\ \vdots \\ a_{m + n}
\end{gmatrix},\ p_0 = a_0.
$$
Since\begin{align*}
&\mathrel{\phantom{=}}{} \begin{gmatrix}
a_0 &&\\
\vdots & \ddots &\\
a_{m + n - 1} & \cdots & a_0
\end{gmatrix} \begin{gmatrix}
q_1 \\ \vdots \\ q_{m + n}
\end{gmatrix} = \begin{gmatrix}
a_0 &&\\
\vdots & \ddots &\\
a_{m + n - 1} & \cdots & a_0
\end{gmatrix} \begin{gmatrix}
q_1 \\ \vdots \\ q_m \\ 0 \\ \vdots \\ 0
\end{gmatrix}\\
&= \begin{gmatrix}
a_0\\
\vdots & \ddots\\
a_{m - 1} & \cdots & a_0\\
a_m & \cdots & a_1 & 0\\
\vdots & \ddots & \vdots & \vdots & \ddots\\
a_{m + n - 1} & \cdots & a_n & 0 & \cdots & 0
\end{gmatrix} \begin{gmatrix}
q_1 \\ \vdots \\ q_m \\ 0 \\ \vdots \\ 0
\end{gmatrix} = \begin{gmatrix}
A & B\\
C & O
\end{gmatrix} \begin{gmatrix}
q_1 \\ \vdots \\ q_m \\ 0 \\ \vdots \\ 0
\end{gmatrix},
\end{align*}
where$$
A = \begin{gmatrix}
a_0\\
\vdots & \ddots\\
a_{m - 1} & \cdots & a_0
\end{gmatrix},\ B = -I_n,\ C = \begin{gmatrix}
a_m & \cdots & a_1\\
\vdots & \ddots & \vdots\\
a_{m + n - 1} & \cdots & a_n
\end{gmatrix},
$$
and $O$ is an $m × n$ zero matrix, and$$
-\begin{gmatrix}
p_1 \\ \vdots \\ p_{m + n}
\end{gmatrix} = -\begin{gmatrix}
p_1 \\ \vdots \\ p_n \\ 0 \\ \vdots \\ 0
\end{gmatrix} = \begin{gmatrix}
A & B\\
C & O
\end{gmatrix} \begin{gmatrix}
0 \\ \vdots \\ 0 \\ p_1 \\ \vdots \\ p_n
\end{gmatrix},
$$
then (1) is equivalent to$$
\begin{gmatrix}
A & B\\
C & O
\end{gmatrix} \begin{gmatrix}
q_1 \\ \vdots \\ q_m \\ p_1 \\ \vdots \\ p_n
\end{gmatrix} = \begin{gmatrix}
a_1 \\ \vdots \\ a_{m + n}
\end{gmatrix},\ p_0 = a_0. \tag{2}
$$
Now, to find the Padé approximation $r(x) = \dfrac{\sum\limits_{k = 0}^2 p_k x^k}{\sum\limits_{k = 0}^3 q_k x^k}$ for $f(x) = \mathrm{e}^{-x} = \sum\limits_{k = 1}^\infty \dfrac{(-1)^k}{k!} x^k$, from (2) there is $q_0 = 1$, $p_0 = 1$, and$$
\begin{gmatrix}
1 &&& -1\\
-1 & 1 &&& -1\\
\dfrac{1}{2} & -1 & 1\\
-\dfrac{1}{6} & \dfrac{1}{2} & -1\\
\dfrac{1}{24} & -\dfrac{1}{6} & \dfrac{1}{2}
\end{gmatrix} \begin{gmatrix}
q_1 \\ q_2 \\q_3 \\ p_1 \\ p_2
\end{gmatrix} = \begin{gmatrix}
-1 \\ \dfrac{1}{2} \\ -\dfrac{1}{6} \\ \dfrac{1}{24} \\ -\dfrac{1}{120}
\end{gmatrix},
$$
thus $q_1 = -\dfrac{3}{5}$, $q_2 = -\dfrac{3}{20}$, $q_3 = -\dfrac{1}{60}$, $p_1 = \dfrac{2}{5}$, $p_2 = -\dfrac{1}{20}$. Therefore,$$
r_{2, 3}(x) = \frac{1 + \dfrac{2}{5}x - \dfrac{1}{20} x^2}{1 - \dfrac{3}{5} x - \dfrac{3}{20} x^2 - \dfrac{1}{60} x^3}.
$$
