Three different methods have been used for computing numerical solutions of Eq. (1) or (2). The first two methods generate solutions forward in time from an arbitrary initial condition; all the steady states to be presented were obtained by both methods. The third method computes a Fourier series expansion of a limit cycle by successive approximation. All reported cases of stable limit cycles were checked with this method as well; furthermore some unstable limit cycles were found with this method only.

A. Taylor series method

The solution is advanced in steps of size h by a two-stage procedure. First the derivatives up to sixth order at t = 0 are obtained via the differential equation in terms of the derivatives to fifth order at t = -L. Then a Taylor series expansion at t = 0 is used to obtain the solution value at t = h in terms of the derivatives up to sixth order at t = 0. This two-stage procedure is then repeated with time t incremented by h, and so on.

B. Validated numerical integration method

Numerical integration can be applied to the integral form Eq. (2); this approach has the advantage that it is possible to obtain rigorous bounds on the error at each step. The trapezoidal rule is the simplest numerical integration. Then, Eq. (2) gives the following form:

,

(3)

where e relates to the error of the numerical integration and can be evaluated as follows:

(4)

or

.

(5)

C. Harmonic balance method of finding a periodic solution

If a solution θ(t) is a limit cycle, its derivative dθ(t)/dt is a periodic function of t; this is true even for a limit cycle of the second kind where θ(t) itself is not a periodic function. Thus a harmonic balance method may be employed. A computer-simulated periodic derivative with period T can be expanded in a trigonometric series as follows:

.

(6)

By integrating this series, θ(t) is obtained as:

.

(7)

From Eq. (7), it follows that a₀ = 0 for a limit cycle of the first kind and that a₀ = 2kπ/T for a limit cycle of the second kind with rotation k.