For structures with non-proportional damping, complex eigenvectors or mode shapes must be used in order to decouple the equations of motion. The resulting equations can then be solved in a systematic way. The necessity of solving a complex eigenvalue problem of a large system remains an obstacle for the practical application of the method. This study utilkes the fact that in practice only a small number of the complex modes are needed. Therefore, these complex modes can be approximated by a linear combination of a small number of the undamped modes, which can be obtained by well established methods with less cost. An additional eigenvalue problem is then solved in a subspace with a much smaller dimension to provide the best combination coefficient for each complex mode. The method of solution for the decoupled equations is then carried over, using the approximate complex modes expressed in undamped mode shapes, to result in simple formulas for the time-and frequency-domain solution. Thus, an efficient modal superposition method is developed for non-proportionally damped systems. The accuracy of this approximate method is studied through an example. By comparing the frequency response result using the approximate method with that using the exact complex modes, it is found that the error is negligible.