In this article, an efficient time–space multigrid (TS-MG) method for solving the harmonic balance (HB) equation system is proposed. The principle of the time–space multigrid method is to coarsen grids in both space and time dimensions simultaneously when coarse grids are formed. The inclusion of time in the time–space multigrid is to address the instability issue or diminished convergence speedup of the spatial multigrid (S-MG) due to larger grid reduced frequencies on coarse grids. With the proposed method, the unsteady flow governing equation will still be solved on all grid levels. Comparing to the finest grid, fewer harmonics and thus fewer equations will be solved consequently on coarse grids. Discrete Fourier transform (DFT) and inverse discrete Fourier transform (IDFT) are used to achieve solution prolongation and restriction between different time grid levels. Results from the proposed method are compared with those obtained from the traditional S-MG and time domain methods. It is found that the TS-MG method can increase solution stability, reduce analysis central processing unit (CPU) time cost required for convergence, save memory usage and has no adverse effect on solution accuracy.