Matrix spherical functions associated to the symmetric pair $(G, K)=$ $\left(\mathrm{SU}(m+2), \mathrm{S}(\mathrm{U}(2) \times \mathrm{U}(m))\right.$, having reduced root system of type $\mathrm{BC}_{2}$ are studied. We consider a $K$-representation $\left(\pi, V_{\pi}\right)$ arising from the $\mathrm{U}(2)$-part of $K$, then the induced representation $\operatorname{Ind}_{K}^{G} \pi$ is multiplicity free. The corresponding spherical functions, i.e. $\Phi: G \rightarrow \operatorname{End}\left(V_{\pi}\right)$ satisfying $\Phi\left(k_{1} g k_{2}\right)=\pi\left(k_{1}\right) \Phi(g) \pi\left(k_{2}\right)$ for all $g \in G, k_{1}, k_{2} \in K$, are studied by studying certain leading coefficients. This is done explicitly using the action of the radial part of the Casimir operator on these functions and their leading coefficients. To suitably grouped matrix spherical functions we associate two-variable matrix orthogonal polynomials giving a matrix analogue of Koornwinder's 1970 s two-variable orthogonal polynomials, which are Heckman-Opdam polynomials for $\mathrm{BC}_{2}$. In particular, we find explicit orthogonality relations and the polynomials being eigenfunctions to a second order matrix partial differential operator. This is joint work with Jie Liu (Radboud $\mathrm{U}$ ).[-]