Chi Square Statistic on Curve Fitting

After looking through the forum for a bit, I discovered there has been a question asked about how to switch out the r^2 goodness of fit on the curve fitting toolbox for a chi^2. The answer was that it is not possible and that you would need to code it yourself. My question is... how do you actually do that?

For some background on my specific problem, I have around 10,000 sets of data, each of which I am curve fitting using the MATLAB curve fitting toolbox. My code is below. As I have so much data, I've painstakingly put it into several cell-type matrices that I have indexed carefully. The equation that I put into the curve fitting function is basically a special one used for light curves of point-sources. I don't think it will matter when addressing my problem.

My current code gives me two matrices: bv (holding the coefficients found by the curve fitting) and goodness_of_fit (holding the gof values that MATLAB also provides).

I know that MATLAB has a built-in chi^2 function, but I am not sure where or how to use it. Do I put it into the curve fitting function of my current code? Or do I calculate it separately later? Some help would be very much appreciated.

function [bv, goodness_of_fit] = createFit(year, dataN, t_0_baseline, bv, goodness_of_fit, coeff_set, goodness_set)

for num = 1:length(dataN{1,year})
    x = dataN{1,year}{1,num}{1,1}(1:end,1);
    y = dataN{1,year}{1,num}{1,1}(1:end,2);

      [xData, yData] = prepareCurveData(x,y);

      %Establish function.
      t_0 = t_0_baseline{1,year}(num,2);
      baseline = t_0_baseline{1,year}(num,3);
      t_0 = num2str(t_0);
      baseline = num2str(baseline);
      name = strcat('-(2.5*log((sqrt((v*(x-',t_0,'))^2+b^2))^2+2)/((sqrt((v*(x-',t_0,'))^2+b^2))*sqrt((sqrt((v*(x-',t_0,'))^2+b^2))^2+4)) -',baseline,')');

      % Set up fittype and options.
      ft = fittype( name, 'independent', 'x', 'dependent', 'y' );                 
      excludedPoints = excludedata( xData, yData, 'range', [0 50]);
      opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
      opts.Display = 'Off';
      opts.StartPoint = [0.972740854003014 0.192028349427775];
      opts.Exclude = excludedPoints;

      % Fit model to data.
      [fitresult, gof] = fit( xData, yData, ft, opts );
      % Get coefficients.
      coeff = coeffvalues(fitresult);
      coeff_set{num} = coeff;

      % Get good of fit.
      goodness = gof;
      goodness_set{num} = goodness;

end
bv{1,year} = coeff_set;
goodness_of_fit{1,year} = goodness_set;
    %{
    % Plot fit with data.
    figure( 'Name', 'untitled fit' );
    h = plot( fitresult, xData, yData, excludedPoints );
    legend( h, 'y vs. x', 'Excluded y vs. x', 'untitled fit', 'Location', 'NorthEast' );
    % Label axes
    xlabel x
    ylabel y
    grid on
    %}
end