System Identification -- Transfer Function from Frequency Response Data

I have some frequency response data (amplitude, phase and frequency)
of a system that I collected using a frequency response analyzer and
am trying to obtain its transfer function (poles and zeros).

This general problem is addressed in
http://www.mathworks.com/support/solutions/data/7423.shtml.  In short,
I tried to use the invfreqs approach suggested in this link, but got
horrible results.  The other suggestion mentioned is to use the
3rd-party Frequency Domain Toolbox ($950).  Other than buying a new
toolbox, does anyone have any suggestions for obtaining the transfer
function from frequency domain data?