Summary: We provide a solution to the \(\beta\)-Jacobi matrix model problem posed by Dumitriu and the first author. The random matrix distribution introduced here, called a matrix model, is related to the model of Killip and Nenciu, but the development is quite different. We start by introducing a new matrix decomposition and an algorithm for computing this decomposition. Then we run the algorithm on a Haar-distributed random matrix to produce the \(\beta\)-Jacobi matrix model. The Jacobi ensemble on \({\mathbb R}^{n}\), parametrized by \(\beta > 0\), \(a > -1\),and \(b > -1\), is the probability distribution whose density is proportional to \(\prod_{i}\lambda_{i}^{({\beta}/{2})(a+1)-1}(1-\lambda_{i})^{({\beta}/{2})(b+1)-1}\prod_{i<j}|\lambda_{i}-\lambda_{j}|^{\beta}\). The matrix model introduced in this paper is a probability distribution on structured orthogonal matrices. If \(J\) is a random matrix drawnfrom this distribution, then a CS decomposition can be taken, \[ J=\left[\begin{matrix} U_1& \\ &U_2\end{matrix} \right]\left[\begin{matrix} C & S \\ -S & C \end{matrix} \right]\left[\begin{matrix} V_1 & \\ & V_2 \end{matrix} \right]^T, \] in which \(C\) and \(S\) are diagonal matrices with entries in \([0,1]. J\) is designed so that the diagonal entries of \(C\), squared, follow the law of the Jacobi ensemble. When \(\beta = 1\) (resp., \(\beta = 2\)), the matrix model is derived by running a numerically inspired algorithm on a Haar-distributed random matrix from the orthogonal (resp., unitary) group. Hence, the matrix model generalizes certain features of the orthogonal and unitary groups beyond \(\beta = 1\) and \(\beta = 2\) to general \(\beta > 0\). Observing a connection between Haar measure on the orthogonal (resp., unitary) group and pairs of real (resp., complex) Gaussian matrices, we find a direct connection between multivariate analysis of variance (MANOVA) and the new matrix model.