A new method of modelling molecular orientational dynamics is presented and illustrated by application to optical second harmonic generation (SHG). The method highlights the intricate dependence of the harmonic signal on the form and evolution of the molecular orientational distribution, and reference is made to some examples taken from the recent literature. Specifically, the cases considered relate to SHG in (a) poled polymers, (b) molecules oriented within molecular sieves and (c) molecules in films or adsorbed on surfaces. The technique invokes rotational averages weighted by static or time-dependent distribution functions expressed in terms of Legendre polynomials. Rotational diffusion is used to model the decay of second harmonic intensity associated with a growth or recovery of bulk isotropy, the results allowing several discrete contributions associated with different time constants to be characterised. For polarisation studies using surfaces and films, it is shown that caution is required in how the results are interpreted to reflect molecular orientation, approximations in commonly employed theory readily producing spurious results.