[guided]The purpose of this step is to isolate the only real computation in the proof. We have reduced the problem to a companion matrix $M$ and must show that the rank-one correction $-e_1e_n^\top q(M)$ changes the characteristic polynomial from $p$ to $q$.
First define the difference polynomial
\begin{align*}
h: \mathbb{R} \to \mathbb{R}, \qquad s \mapsto q(s)-p(s).
\end{align*}
Because both $p$ and $q$ are monic of degree $n$, the leading terms cancel, so $\deg h \le n-1$. Also $p(M)=0$ for this companion matrix: the identity $M^{j-1}e_1=e_j$ holds for $1 \le j \le n$, and the final-column relation gives exactly
\begin{align*}
M^n e_1 + c_{n-1}M^{n-1}e_1+\cdots+c_1Me_1+c_0e_1=0.
\end{align*}
Since $e_1,Me_1,\dots,M^{n-1}e_1$ form the standard basis, we must still check the identity on each vector of this basis. The matrix $p(M)$ commutes with $M$ because $p(M)$ is a polynomial in $M$. Hence, for every integer $k$ with $0\le k\le n-1$,
\begin{align*}
p(M)M^k e_1=M^k p(M)e_1=0.
\end{align*}
Thus $p(M)$ vanishes on the basis $e_1,Me_1,\dots,M^{n-1}e_1$, and therefore $p(M)=0$. Consequently $q(M)=h(M)$.
Let $D := \{r \in \mathbb{R}: p(r) \ne 0\}$. For $s \in D$, the condition $p(s)\ne 0$ is exactly the condition that $sI_n-M$ is invertible. Define the resolvent-vector map
\begin{align*}
x: D \to \mathbb{R}^n, \qquad s \mapsto (sI_n-M)^{-1}e_1.
\end{align*}
Fix $s \in D$.
The scalar quantity we need is $e_n^\top h(M)x(s)$. Since $h$ has degree at most $n-1$, it is enough to understand $e_n^\top M^j x(s)$ for $0 \le j \le n-1$.
Write
\begin{align*}
x(s)=x_1(s)e_1+x_2(s)e_2+\cdots+x_n(s)e_n,
\end{align*}
where each $x_i:D\to\mathbb R$ is the $i$-th coordinate function. Expanding $(sI_n-M)x(s)=e_1$ using the companion relations for $M$ gives
\begin{align*}
sx_1(s)+c_0x_n(s)=1,
\end{align*}
and, for $2\le i\le n$,
\begin{align*}
-x_{i-1}(s)+sx_i(s)+c_{i-1}x_n(s)=0.
\end{align*}
The second recurrence is best solved backward because it expresses $x_{i-1}(s)$ in terms of $x_i(s)$ and $x_n(s)$:
\begin{align*}
x_{i-1}(s)=sx_i(s)+c_{i-1}x_n(s).
\end{align*}
Iterating this identity gives, for every integer $k$ with $1\le k\le n-1$,
\begin{align*}
x_{n-k}(s)=\left(s^k+c_{n-1}s^{k-1}+\cdots+c_{n-k}\right)x_n(s).
\end{align*}
In particular,
\begin{align*}
x_1(s)=\left(s^{n-1}+c_{n-1}s^{n-2}+\cdots+c_1\right)x_n(s).
\end{align*}
Substituting this into the first coordinate equation gives
\begin{align*}
\left(s^n+c_{n-1}s^{n-1}+\cdots+c_1s+c_0\right)x_n(s)=1.
\end{align*}
The polynomial in parentheses is $p(s)$, so
\begin{align*}
e_n^\top x(s)=x_n(s)=\frac{1}{p(s)}.
\end{align*}
The remaining recurrence comes from rewriting the same equation as
\begin{align*}
Mx(s)=sx(s)-e_1.
\end{align*}
For $1 \le j \le n-1$, multiply by $M^{j-1}$ and then apply $e_n^\top$:
\begin{align*}
e_n^\top M^j x(s)=s\,e_n^\top M^{j-1}x(s)-e_n^\top M^{j-1}e_1.
\end{align*}
But $M^{j-1}e_1=e_j$, and $e_n^\top e_j=0$ for $j<n$. Thus
\begin{align*}
e_n^\top M^j x(s)=s\,e_n^\top M^{j-1}x(s).
\end{align*}
Starting from $e_n^\top x(s)=1/p(s)$, induction gives
\begin{align*}
e_n^\top M^j x(s)=\frac{s^j}{p(s)}
\end{align*}
for every $0 \le j \le n-1$.
Write
\begin{align*}
h(s)=h_0+h_1s+\cdots+h_{n-1}s^{n-1}.
\end{align*}
Using linearity and the preceding moment computation,
\begin{align*}
e_n^\top h(M)(sI_n-M)^{-1}e_1
=
\sum_{j=0}^{n-1} h_j e_n^\top M^j x(s)
=
\frac{h(s)}{p(s)}.
\end{align*}
Finally set $N:=sI_n-M$, $u:=e_1$, and $v:=h(M)^\top e_n$. Since $s\in D$, the matrix $N$ is invertible. We factor $N$ from the rank-one perturbation:
\begin{align*}
\det(N+uv^\top)=\det(N)\det(I_n+N^{-1}uv^\top).
\end{align*}
Set $a:=N^{-1}u \in \mathbb{R}^n$ and $b:=v \in \mathbb{R}^n$. We now compute the rank-one determinant directly. Write $a_j$ and $b_j$ for the $j$-th coordinates of $a$ and $b$. The $j$-th column of $I_n+ab^\top$ is $e_j+b_j a$. By multilinearity of the determinant in the columns, expanding the product produces the determinant of $I_n$, the $n$ terms with exactly one column replaced by $a$, and terms with at least two columns equal to scalar multiples of $a$. The latter terms vanish by alternation. Therefore
\begin{align*}
\det(I_n+ab^\top)=1+\sum_{j=1}^n b_j\det(e_1,\dots,e_{j-1},a,e_{j+1},\dots,e_n).
\end{align*}
The displayed determinant with $a$ in the $j$-th column equals $a_j$, so
\begin{align*}
\det(I_n+ab^\top)=1+\sum_{j=1}^n b_j a_j=1+b^\top a.
\end{align*}
Applying this with $a=N^{-1}u$ and $b=v$ yields
\begin{align*}
\det(sI_n-M+e_1e_n^\top h(M))
=
\det(sI_n-M)\left(1+e_n^\top h(M)(sI_n-M)^{-1}e_1\right).
\end{align*}
Since $\det(sI_n-M)=p(s)$, this becomes
\begin{align*}
\det(sI_n-M+e_1e_n^\top h(M))
=
p(s)\left(1+\frac{h(s)}{p(s)}\right)
=
p(s)+h(s)
=
q(s).
\end{align*}
This proves the identity for all $s$ with $p(s)\ne 0$. The two sides are polynomials in $s$, and a nonzero polynomial has only finitely many roots, so equality on the infinite set $\{s \in \mathbb{R}:p(s)\ne 0\}$ implies equality for every $s \in \mathbb{R}$.[/guided]
