We study a second order ensemble method for fast computation of an ensemble of magnetohydrodynamics flows at small magnetic Reynolds number. Computing an ensemble of flow equations with different input parameters is a common procedure for uncertainty quantification in many engineering applications, for which the computational cost can be prohibitively expensive for nonlinear complex systems. We propose an ensemble algorithm that requires only solving one linear system with multiple right-hands instead of solving multiple different linear systems, which significantly reduces the computational cost and simulation time. Comprehensive stability and error analyses are presented proving conditional stability and second order in time convergent. Numerical tests are provided to illustrate theoretical results and demonstrate the efficiency of the proposed algorithm.