Receiving indexing errors on curve fitting code which has been working previously before changing the fitting equation. In desperate need of advice..

I'm trying to fit some equations to some data. I am doing this by using the curve fitting tool, generating the code, and then I call the function inside a for loop to fit the equation to 'k' sets of data (the number of data sets I import into Matlab). Below is the way I've been defining the equations and calling them. I've been working with these fits for months, and I've just found out I have to fit some different equations to the data. So literally all I've changed in the code is the fitting equation 'ft' and hence the expression for 'yfitExp', however now it's giving me an error stating that brace indexing does not support my 'yfitExp' variable. When I use parentheses (so 'yfitExp(k)') it give me the error message "Unable to perform assignment because the left and right sides have a different number of elements.".

I will say my gamma fit seems to be working fine so I've included the code for this below the exponential fit. All of the other fits apart from gamma are producing the same error messages all of the sudden but I really can't pick where I've gone wrong. I know its obviously having an issue with the type of variable 'yfitExp' is (its a double) but I can't pick what I've done differently on the exponential fit to make it produce errors.

If anyone can provide some input I would be really grateful as I am very desperate to solve this issue (I clearly have a lot of work to do on my understanding of indexing and loops because I'm hopeless at them haha).

%==============================
%   EXPONENTIAL DISTRIBUTION
%==============================
for k = 1:length(theFiles)
    txt = [theFiles(k).name];
    col = lines(length(theFiles));
  
    [yfitExp{k}, gof_exp(k), a_exp(k), b_exp(k)] = expFit(dose{k}, sf{k});
    
    % Plot fit with data.
    figure(14);
    hold on
    plot(dose{k}, sf{k},'Color',col(k,:),'Marker','x','LineStyle','none','DisplayName',txt);
    plot(dose{k}, yfitExp{k},'Color',col(k,:),'HandleVisibility','off')
    h = errorbar(dose{k},sf{k},errs{k},'Marker','none','LineStyle','none','Color',[col(k,1), col(k,2), col(k,3)],'HandleVisibility','off');
    legend('Location','southoutside');
    % Label axes
    xlabel( 'Dose (Gy)' ); ylabel( 'Surviving Fraction' );
    title('Exponential Distribution Fit - Nonlinear Least Squares');
    grid on; grid minor; set(gcf,'color','w');
end

%==================================
%   EXPONENTIAL DISTRIBUTION FIT
%==================================
function [yfitExp, gof, a, b] = expFit(x, y)

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

% Set up fittype and options.
ft = fittype( 'a*exp(-a*x) + b', 'independent', 'x', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [0.983170757650345 0.830932322999125];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

coeffs = coeffvalues(fitresult);
a = coeffs(1);
b = coeffs(2);

yfitExp = a*exp(-a*x) + b;
end

%=========================
%   GAMMA DISTRIBUTION
%=========================
for k = 1:length(theFiles)
    txt = [theFiles(k).name];
    col = lines(length(theFiles));
    
    [yfitGamma{k}, gof_gamma(k), a_gamma(k), b_gamma(k)] = gammaFit(dose{k}, sf{k});
    
    % Plot fit with data.
    figure(13);
    hold on
    plot(dose{k}, sf{k},'Color',col(k,:),'Marker','x','LineStyle','none','DisplayName',txt);
    plot(dose{k}, yfitGamma{k},'Color',col(k,:),'HandleVisibility','off')
    h = errorbar(dose{k},sf{k},errs{k},'Marker','none','LineStyle','none','Color',[col(k,1), col(k,2), col(k,3)],'HandleVisibility','off');
    legend('Location','southoutside');
    % Label axes
    xlabel( 'Dose (Gy)' ); ylabel( 'Surviving Fraction' );
    title('Gamma Distribution Fit - Nonlinear Least Squares');
    grid on; grid minor; set(gcf,'color','w');
end

%=============================
%   GAMMA DISTRIBUTION FIT
%=============================
function [yfitGamma, gof, a, b] = gammaFit(x, y)

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

% Set up fittype and options.
ft = fittype( '(b^a/gamma(a))*x.^(a-1).*exp(-b*x)', 'independent', 'x', 'dependent', 'y' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.StartPoint = [0.983170757650345 0.830932322999125];
opts.Upper = [2 2];
opts.Lower = [0 0];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

coeffs = coeffvalues(fitresult);
a = coeffs(1);
b = coeffs(2);

yfitGamma = (b^a/gamma(a))*x.^(a-1).*exp(-b*x);
end