A pair of regular matrix functionals \(\{\textbf{u}_0, \textbf{u}_1\}\) is said to be a matrix left coherent pair if their corresponding sequences of matrix orthogonal polynomials \(\{P_n \}_{\ge 0}\) and \(\{T_n \}_{\ge 0}\) satisfy the following structure relation: \(\begin{aligned} T_n(x)=\frac{1}{n+1}P^{\prime }_{n+1}(x)+\frac{1}{n}\sigma _nP^{\prime }_{n}(x), \ n\ge 1, \end{aligned}\) where \(\sigma _n\) is a complex matrix for every natural number n. In this contribution, we build and implement algorithms, in the framework of matrix coherence, that allow the explicit computation of sequences of matrix orthogonal polynomials associated with matrix left coherent pairs and Sobolev matrix orthogonal polynomials. Moreover, we apply these algorithms to determine the Fourier–Sobolev matrix-valued coefficients associated with a matrix function living in an appropriate function space in terms of the Sobolev matrix orthogonal polynomials with respect to a Hermitian matrix weight. We illustrate the obtained results and algorithms with some numerical examples.