Expressing non-homogeneous differential equation into a system of first order differential equations The differential equation that describes my system is given as 
$y^{(n)} + a_{n-1}y^{(n-1)}+\cdots +a_1\dot{y} +a_0y = b_{n-1}u^{(n-1)} + \cdots +b_1\dot{u} + b_0u + g(y(t),u(t))$
I want to express the above differential equation into a system of linear differential equations of the form
$\dot{x} = Ax+ Bu + B_pg$
$y = Cx$
The matrices are given as follows: However, I am not able to prove, how to get them
$A = \begin{bmatrix} 0 & 1 & 0&\cdots &0\\
0&0&1&\cdots&0\\
\vdots &\vdots & \vdots & \vdots & \vdots \\
-a_0 & -a_1&-a_2&\cdots & -a_{n-1}
\end{bmatrix}$ 
$B = \begin{bmatrix} 0 \\0\\ \vdots\\b_0
\end{bmatrix}$
$C =\begin{bmatrix}
1 &b_1/b_0 & b_2/b_0 & \cdots &b_{n-1}/b_0
\end{bmatrix}$
How do I get the above matrices from the differential equation form as shown above?
 A: Despite the fact that the $(n-1)$-th order ODE is not "simple looking" as the ones usually shown in the examples of transformation into a first order system, the process is absolutely the same. Defining the $n$-dimensional vectors $\mathbf{u} $ and $\mathbf{y}$ as
$$
\begin{split}
\mathbf{u} &=(u_{n-1},u_{n-2}\ldots, u_1, u_0)= \Big(u^{(n-1)}, u^{(n-2)}, \ldots, \dot{u} , u\Big)\\
\mathbf{y} &=(y_{n-1},y_{n-2}\ldots, y_1, y_0)= \Big(y^{(n-1)}, y^{(n-2)}, \ldots, \dot{y\,} , y\Big)\\
\end{split}
$$
we have that the first order ODE system equivalent to the 
$$
\begin{align}
\dot{y_0} & = y_1\\
\dot{y_1} & = y_2\\
\vdots\,&=\vdots\\
\dot{y\,}_{\!n-2} & = y_{n-1}\\
\dot{y\,}_{\!n-1} & = -a_{n-1}y_{n-1} -\ldots -a_1y_1 +a_0y_0 + b_{n-1}u_{n-1} + \ldots +b_1u_1 + b_0u_0 + g\big(\mathbf{y}(t),\mathbf{u}(t)\big)
\end{align}
$$
i.e.
$$
\begin{split}
\dot{\mathbf{y}} &=
\begin{pmatrix}
0 & 1 & 0 & \ldots & 0\\
0 & 0 & 1 & \ldots & 0\\
\vdots & \vdots & \vdots & \ldots & \vdots\\
0 & 0 & 0 & \ldots & 1\\
-a_0 & -a_1 & -a_2 & \ldots & -a_{n-1}\\
\end{pmatrix}\mathbf{y}\\
\\
&\qquad +
\begin{pmatrix}
0 & 0 & 0 & \ldots & 0\\
0 & 0 & 0 & \ldots & 0\\
\vdots & \vdots & \vdots & \ldots & \vdots\\
0 & 0 & 0 & \ldots & 0\\
b_0 & b_1 & b_2 & \ldots & b_{n-1}\\
\end{pmatrix}\mathbf{u} + 
\begin{pmatrix}
0 \\
0 \\
\vdots \\
0 \\
g\big(\mathbf{y}(t),\mathbf{u}(t)\big)
\end{pmatrix}.\\
\end{split}\label{1}\tag{1}
$$
and the matrices $A$, $B$, are immediately acknowledged: I have kept the last term of the system in a vector form since its structure is clearer. However, you can express it as a product of the vector 
$$
\mathbf{g}=\underbrace{(0,0,\ldots,g)}_{n}
$$ 
and the matrix $B_p\in M^{n\times n}$ whose all entries are $0$ except the entry $b^p_{nn}$ with is set to $1$. 
A: The aim of such a transformation is to eliminate the derivatives of $u$, as the control input is usually just a stream of $u$ values, possibly with noise, and numerical differentiation is unreliable.
In a first step, move all derivatives on the left side
$$
\frac{d}{dt}[\underbrace{y^{(n-1)}+a_{n-1}y^{(n-2)}+a_1y-b_{n-1}u^{(n-2)}-...-b_1u}_{=x_{n-1}}]=-a_0y+b_0u+g(y,u)
$$
Call the term inside the derivative on the left side $x_{n-1}$ (index $=$ highest $y$ derivative order). Then in the same way
$$
\frac{d}{dt}[\underbrace{y^{(n-2)}+a_{n-1}y^{(n-3)}+a_2y-b_{n-1}u^{(n-3)}-...-b_2u}_{=x_{n-2}}]=x_{n-1}-a_1y+b_1u
$$
etc. until
$$
\frac{d}{dt}y=x_1-a_{n-1}y+b_{n-1}u
$$
so that  the resulting system is
$$
\frac{d}{dt}
\left[\begin{alignedat}{1}
y\\x_1\\x_2\\\vdots\\x_{n-2}\\x_{n-1}
\end{alignedat}\right]
=
\left[\begin{alignedat}{1}
x_1&-a_{n-1}y&+b_{n-1}u&\\
x_2&-a_{n-2}y&+b_{n-2}u&\\
x_3&-a_{n-3}y&+b_{n-3}u\\
\vdots~~&\\
x_{n-1}&-a_1y&+b_1u\\
&-a_0y&+b_0u&+g(y,u)
\end{alignedat}\right]
\\~\\
=
\begin{bmatrix}
-a_{n-1}&1&0&\\
-a_{n-2}&0&1&0&\\
-a_{n-3}&0&0&1&0\\
\vdots&&&\ddots&\ddots&\\
-a_1&0&&&0&1\\
-a_0&0&&&&0
\end{bmatrix}
\begin{bmatrix}y\\x_1\\x_2\\\vdots\\x_{n-2}\\x_{n-1}\end{bmatrix}
+
\begin{bmatrix}
b_{n-1}\\b_{n-2}\\b_{n-3}\\\vdots\\b_1\\b_0
\end{bmatrix}u
+
\begin{bmatrix}
0\\0\\ \\\vdots\\0\\g(y,u)
\end{bmatrix}
$$
which has more-or-less the structure in question, however with the transpose of the companion matrix and different use of the $b$ vectors.

This is different from the proposed solution. There, the idea seem to be to set $y=\tilde b(D)x$, where $D$ is the differentiation operator and $\tilde b_k=b_k/b_0$. Then the equation reads as
$$
\tilde b(D)[a(D)x-b_0u]=g(y,u)
$$
Now it is easy to transform $a(D)x-b_0u=v$ into the usual matrix form using $\vec x=[x_0,...,x_{n-1}]^T$ with $x_k=x^{(k)}$. It is less clear how $\tilde b(D)v=g(y,u)$ results in a form with some matrix $B_p$.
