I am reading about adjoint systems of linear partial differential equations and appears the following example of a system over the domain $(a,b)\times (0,T)$, where $T>0$ and $a,\, b\in \mathbb{R}$:
$$
f_t + a_1 g_x =0,
$$
$$g_t - a_2 g_{xx} + a_3 f_x + a_4 h_x =0, $$
$$
h_t - a_5 h_{xx}+ a_6 g_x =0,
$$
where each $a_i$ is a real number. The adjoint system associated to the system above reads as
$$ - u_t - a_1 v_x =0, $$ $$ - v_t - a_2 v_{xx} - a_3 u_x - a_4 w_x =0, $$ $$ - w_t - a_5 w_{xx} - a_6 v_x =0.$$
How could I obtain this adjoint system ? What is the procedure to get it ?
You can write your original system as $A \vec{f} = \vec{0}$ where $\vec{f} = (f,g,h)$ and
\begin{align} A &= \begin{pmatrix} \partial_{t} & a_{1} \partial_{x} & 0 \\ a_{3} \partial_{x} & \partial_{t} - a_{2} \partial_{xx} & a_{4} \partial_{x} \\ 0 & a_{6} \partial_{x} & \partial_{t} - a_{5} \partial_{xx} \end{pmatrix} \end{align}
is a matrix of disembodied derivative operators. Recall that $M$ is adjoint to $N$ if $\langle M x, y \rangle = \langle x, N y \rangle$ for a suitable notion of $\langle \cdot, \cdot \rangle$. In this case, we want to use integration by parts and so we let $\langle \cdot, \cdot \rangle = L^{1}$, though we will also consider $\langle \cdot, \cdot \rangle$ to be a projection (dot product) so we can solve the entire system simultaneously.
Computing $\langle A \vec{f}, \vec{u} \rangle$ where $\vec{u} = (u, v, w)$, $\Omega = (a, b) \times (0, T)$ and IBP stands for integration by parts
\begin{align} 0 &= \langle A \vec{f}, \vec{u} \rangle \quad \text{(as $A \vec{f} = \vec{0}$)} \\\\ &= \int_{\Omega} (A \vec{f}) \cdot \vec{u} \ dx dt \\\\ &=\int_{\Omega} \left[ (f_{t} + a_{1} g_{x}) u + (a_{3} f_{x} + g_{t} - a_{2} g_{xx} + a_{4} h_{x}) v + (a_{6} g_{x} + h_{t} - a_{5} h_{xx}) w \right] \ dx dt \\\\ &\underbrace{=}_{\text{IBP}} \text{(terms that need to be evaluated)} \\ &- \int_{\Omega} \left[ f u_{t} + g (a_{1} u_{x}) + f (a_{3} v_{x}) + g (v_{t} + a_{2} v_{xx}) + h (a_{4} v_{x}) + g (a_{6} w_{x}) + h (w_{t} + a_{5} w_{xx}) \right] dx dt \\\\ &= \text{(terms that need to be evaluated)} \\ &- \int_{\Omega} \left[ f (u_{t} + a_{3} v_{x}) + g (a_{1} u_{x} + v_{t} + a_{2} v_{xx} + a_{6} w_{x}) + h (a_{4} v_{x} + w_{t} + a_{5} w_{xx}) \right] dx dt \\\\ &= \text{(terms that need to be evaluated)} + \int_{\Omega} \vec{f} \cdot (B \vec{u}) \ dx dt \\\\ &= \text{(terms that need to be evaluated)} + \langle \vec{f}, B \vec{u} \rangle \end{align}
where
\begin{equation} B = \begin{pmatrix} - \partial_{t} & - a_{3} \partial_{x} & 0 \\ -a_{1} \partial_{x} & - \partial_{t} - a_{2} \partial_{xx} & -a_{6} \partial_{x} \\ 0 & -a_{4} \partial_{x} & - \partial_{t} - a_{5} \partial_{xx} \end{pmatrix} \end{equation}
which is almost the result you have written down, except the coefficients $a_{1}$ and $a_{3}$ are swapped, as are the $a_{4}$ and $a_{6}$ coefficients, which is a typo on your end or the papers. You also need to specify boundary and initial conditions on the functions to determine what happens to the evaluated terms.