An automated algorithm is presented which enables the harmonic balance equations for any polynomial type pure or cross-product non-linear differential system to be written down directly in terms of the coefficients of the governing equation and the complex amplitudes of general harmonic waveform. The system frequency response, in the form of multi-input, amplitude dependent describing function, is therefore readily computed. The method is illustrated by means of an example, and the results validated against detailed numeric simulation.