This study investigated $\eta$-Hermitian solutions to a two-sided matrix equation and a system of matrix equations over the skew-field of quaternions. We employed solvability conditions to identify the general solution for this system. This method deduces the Moore-Penrose inverse, and the rank equalities for the system's coefficients. We utilized these strategies to develop an algorithm that can compute general solutions. We also used these algorithms in numerical instances to verify the theoretical results.