A systematic algorithm for generating the polyharmonic balance equations for any system within a broad class of time-delayed differential, or non-linear difference equation models, is presented. The method, which is readily automated, enables the balance equations to be written down directly in terms of the coefficients of the governing equation, and the complex amplitudes of a general harmonic waveform. The system frequency response or amplitude dependent describing function is then readily computed. The method is illustrated by means of examples including both a time-delayed differential system example and a discrete time NARX model application. The results are validated against detailed numeric simulation which confirms the accuracy and efficacy of the approach.