We perform nonperturbative calculations of resonant multiwave mixing of two high-intensity laser beams of different frequencies both for a general three-level system and for the 3 Siy2 -3 Piy2 electronic resonance of the sodium atom. Our calculations proceed by direct numerical integration of the densitymatrix equations and subsequent Fourier analysis of the calculated time-dependent density-matrix elements. We examine the case where two nearly degenerate input beams are each detuned many collisional linewidths from electronic resonance, and the frequency difference between the input beams is tuned through a ground-state resonance. Unlike previous high-intensity analyses, each input beam can be resonant with numerous transitions simultaneously; for a three-level system, we refer to this as an "M-type" system, in analogy with Aand V-type systems. Calculated four-, six-, and eight-wave-mixing spectra exhibit subharmouic resonances of the type observed experimentally by Trebino and Rahu [Opt. Lett. 12, 912 (1987)]. The three-level system is studied to gain insight into the physics of these subharmouic resonances. For more quantitative comparison with experiment, we use a model that includes the 16 Zeeman and hyper5ne states in the 3 S&~2 and the 3 P&&2 levels, and accounts for arbitrary linear polarization of the input laser radiation. The relative intensities of the resonance features in the calculated multiwave-mixing spectra show good agreement with experiment. These resonance features broaden and shift as laser intensity increases, as predicted for four-wave-mixing processes in previous, less general high-intensity analyses. Both theory and experiment also show additional resonant features in the eight-wave-mixing spectra that appear at very high intensity.