A general method for computing harmonic generation in reflection and transmission from planar nonmagnetic multilayer structures is described. The method assumes plane waves and treats harmonic generation in the parametric approximation. The method is applied in studying the second- and third-harmonic generation properties of thin crystal silicon layers surrounded by thermal oxide. Most independent components of the nonlinear susceptibility tensor have unique signatures with silicon layer thickness $d$, allowing their strength to be determined in principle by measuring harmonic generation as a function of $d$. Surface and bulk contributions to third-harmonic generation are cleanly distinguished, with the bulk signal dominating. Four of six nonvanishing components of ${\ensuremath{\chi}}^{(2)}$ are independent. An approximate value for the bulk susceptibility component ${\ensuremath{\delta}}^{\ensuremath{'}}$, which is accessible only in multibeam experiments and has not previously been measured, is obtained.