The simulation of time-periodic unsteady flows is a central problem in aeronautical applications, especially in turbomachinery. The so-called harmonic balance (HB) method which uses a spectral discretisation of the time derivative has been shown to be a highly efficient approach for applications in unsteady aerodynamics and non-linear aeroelasticity. Unlike linearised frequency-domain methods, HB takes the non-linear interaction between harmonics into account. In contrast to other disciplines (e.g. electrical circuit analysis or structural dynamics), all HB solvers in Computational Fluid Dynamics (CFD) seem to use pseudotime stepping, thereby adopting the traditional approach to achieve steady solutions. In the authors’ experience, HB together with pseudotime stepping can give unsteady solutions of high accuracy at moderate costs, provided the solver converges. There are, however, occasionally configurations where, at least for some operating conditions, it seems extremely hard to achieve converged HB solutions, which raises the question of the optimal solution technique. In this paper, we give a physical motivation for pseudotime stepping. We show that, even for highly non-linear problems, pseudotime marching HB solvers inherit important properties from the standard time-integration approach. Roughly speaking, we show that along certain lines in the pseudotime-time plane the pseudotime HB solution corresponds to a discrete solution of the original ordinary differential equation. This shows that, given sufficiently many harmonics and small pseudotime steps, the HB solver should converge to asymptotically periodic solutions provided the initial solution is appropriate. On the other hand, we see that self-sustained flow instabilities can prevent the HB solver from converging. We illustrate our results by means of the van der Pol oscillator as well as unsteady flow problems for a NACA profile.