How does Curve Fitting Toolbox's fit() function works with non-differentiable expressions?

As far as I understand, fit() uses Levenberg-Marquardt's method of function optimization, which is a gradient-based method (modified to incorporate a second-order term).This essentially means that we need a way to compute a gradient of the objective function to find the direction of the next step.(With respect to the parameter variables.)

If I'm trying to optimize a complex Matlab expression, possibly including classes, how does Matlab compute this gradient?

I can imaging computing the gradient numerically, or using the backpropagation (chain rule). But how does Matlab actually do it?

I have a relatively difficult function to fit (20 parameters), but curiously, if I set the initial guess not too far from the true value (when I'm doing the testing), the optimization does converge to the true values of the parameters.

I'm using Matlab 2013a.

Thanks a lot in advance.