This manuscript provides a comprehensive study of how fitting generalized Kubo lineshapes to multidimensional spectra compares to the CLS method in terms of the accuracy, precision and reliability of the resulting parameters. While Garret-Roe and coworkers have shown several examples of model fitting 2D IR waiting-time series using the fmincon function in MATLAB, a comparison of accuracy between model fitting and the CLS method by fitting to simulated spectra, where true values of parameters are known, has been missing. We find our model fitting routine improves precision over the CLS method by 8× to 15× on average for Kubo time constants and 8× to 50× for Kubo amplitudes and homogeneous dephasing, which is due in part to a novel figure of merit used in our fitting algorithm that we refer to as the scale invariant gradient norm (SIGN). We find that numerical instabilities associated with some fitting parameters appearing nearly indistinguishable at certain points in parameter space are the primary cause of sudden ceasing, which could be mistaken for local minima. Importantly, the SIGN readily identifies these events, which enables swift correction by our fitting algorithm.