I have a general form equation
$$\ddot{y} + a_1 \dot{y} + a_2 y = b_0 \ddot{u} + b_1 \dot{u} + b_2 u$$ in which I am trying to find the state space representation.
$$\dot{X} = AX + BU$$
$$Y=CX +DU$$
to find the matrix $B$ I have found online that many use this method
$\beta_0 = b_0$
$\beta_1 = b_1 - a_1 \beta_0$
$\beta_2 = b_2 - a_1 \beta_1 - a_2 \beta_0$
where $B = \begin{bmatrix}\beta_1 \\ \beta_2\end{bmatrix}$ and $D = \beta_0$
I can compute this of course, but I do not understand how this works. I am looking for an answer that shows how this works or derives this method.
First of all, keep in mind that you can have several state-space representations. That is, you have an infinite number of ways to define the state variables $x$ and you will get a different $B$ matrix for each. But, I will attempt to answer your question.
First, find the transfer function of the system by taking Laplace Transform (assuming zero initial conditions): $$ \frac{Y(s)}{U(s)} = \frac{b_0s^2+b_1s+b_0}{s^2+a_1s + a_2}$$ Since this is an improper transfer function, divide by the denominator to simplify: $$ \frac{Y(s)}{U(s)} = \frac{(b_1-a_1b_0)s+(b_2-b_0a_2)}{s^2+a_1s + a_2}+b_0$$
From this transfer function, you can build a state-space representation by using visual inspection. Since the system is a second order, your $A$ matrix is a 2x2, $B$ is a 2x1, $C$ is a 1x2, and $D$ is scalar.
Now, I do not know which representation are they using and without the $A$ and $C$ formulas I think is impossible to tell. But, for instance, if you follow the observable canonical form, you will get that your $B$ is the coefficients $[(b_1-a_1b_0), (b_2-a_2b_0)]^T$ of the numerator. The $D$ matrix is the constant coefficient $b_0$ (always the case, regardless of the state-space form you choose).