On the numerical solution of coupled eigenvalue differential equations arising in molecular spectroscopy