We solve the van der Pol oscillator only up to order 2. This method can be continued indefinitely in the same way, where the order-n term \epsilon^n x_n consists of a harmonic term a_n\cos(t) + b_n\cos(t), plus some super-harmonic terms a_{n, 2}\cos(2t) + b_{n, 2}\cos(2t) + \cdots. The coefficients of the super-harmonic terms are solved directly, and the coefficients of the harmonic term are determined by expanding down to order-(n+1), and eliminating its secular term.
