© 2018 Society for Industrial and Applied Mathematics. We consider a set of linear hyperbolic equations coupled by a linear source term and introduce a surface hopping Gaussian beam method as its numerical solver. The Gaussian beam part is basically classic, while the surface hopping part is derived from the equations. The whole algorithm shows high efficiency and good parallelizability. An application on the quantum-classical Liouville equations is presented to show its potential use in practice.