For problems with very fine surface meshes, typically the most time-consuming step of a boundary element method (BEM, also called a panel method) is solving the final linear system of equations. Many have already studied how to efficiently solve the dense, asymmetric systems which arise in elliptic BEMs. However, this has not been studied for a supersonic aerodynamic BEM, for which the governing PDE is hyperbolic. Due to this hyperbolic character, the matrix equation which arises from a supersonic BEM has a large number of identically zero elements. But the resulting linear system of equations is also not sparse in the standard sense. Hence, the efficient solution of the linear system of equations arising in a supersonic BEM is here considered. A novel sorting algorithm is developed whereby the non-zero elements may be arranged into a useful structure with minimal cost. A novel direct solution method is developed here based on fast Givens rotations and the QR decomposition. This novel solver leverages the unique matrix structure to solve the supersonic system of equations more quickly than traditional direct methods. This novel method is then compared to other direct and iterative matrix solvers and is shown to be more robust than iterative solvers and more efficient than other direct solvers, with a computational time complexity of approximately O(N2.5).