PolyHarmonic Distortion modeling in MATLAB?

PolyHarmonic Distortion modeling approach is a black-box frequency-domain model that provides a foundation for measurement, modeling, and simulation of driven nonlinear systems. Basically, the model uses the concept of harmonic superposition to analytically describe the large-signal B waves' response of a nonlinear system, linearized around a large-signal point, in terms of a linear mapping of the stimulus A waves, similar to classic S parameters. Since it relies heavily on arrays and matrixes, it seems to be a good candidate for implementation in MATLAB. Agilent has a trademark on their own name for PolyHarmonic Distortion modeling; however they do not have intellectual property rights for the technique itself. Has anyone done any MATLAB coding using PolyHarmonic Distortion modeling?