  We introduce matrix-valued weight functions of arbitrary size, which are
analogues of the weight function for the Gegenbauer or ultraspherical
polynomials for the parameter $\nu>0$. The LDU-decomposition of the weight is
explicitly given in terms of Gegenbauer polynomials. We establish a
matrix-valued Pearson equation for these matrix weights leading to explicit
shift operators relating the weights for parameter $\nu$ and $\nu+1$. The
matrix coefficients of the Pearson equation are obtained using a special
matrix-valued differential operator in a commutative algebra of symmetric
differential operators.
  The corresponding orthogonal polynomials are the matrix-valued Gegenbauer
polynomials which are eigenfunctions for the symmetric matrix-valued
differential operators. Using the shift operators we find the squared norm and
we establish a simple Rodrigues formula. The three-term recurrence relation is
obtained explicitly using the shift operators as well. We give an explicit
non-trivial expression for the matrix entries of the matrix-valued Gegenbauer
polynomials in terms of scalar-valued Gegenbauer and Racah polynomials using
the LDU-decomposition and differential operators. The case $\nu=1$ reduces to
the case of matrix-valued Chebyshev polynomials previously obtained using group
theoretic considerations.
