The options structure to lsqcurvefit must be the 
argument to it, since (except for name-value pair arguments), arguments to MATLAB functions are positional.
argument to it, since (except for name-value pair arguments), arguments to MATLAB functions are positional. So in your call to lsqcurvefit, use:
[y, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@MB, y0, E_fit, delta_Abs_norm_fit, lb, ub, [], [], [], [], [], optim_lsq);
I tested that and it worked. (I also added a ga call to get a better initial estimate of ‘y0’, and then let lsqcurvefit tweak it.)
%% Preparation
clear;clc
data = importdata("FCPIB-293K-2.5...bg -chirp.csv"); % insert file path within parenthesis
%% Preamble
% Fundamental constants
h = 4.0135667696*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K
% Clean up of data to select range of values
wavelength = data(1:end, 1);
delay_t = data(1, 1:end); % conatains all of the delay times
E = (h*c)./(wavelength*10^-9); % contains all of the probe energies
Range_E = E>=1.5 & E<=2.2;
Range_T = delay_t>=0.5 & delay_t<=1000;
% for one delay time
T = find(Range_T);
T_min = min(T);
T_max = max(T);
t = 57; % choose an integer b/w T_min and T_max
delaytime = delay_t(1, t);
disp(delaytime)
% Initial parameter guess and bounds
lb = [0, 293, -1]; ub = [Inf, 1200, 1];
y0 = [2*10^9, 1000, 0.5];
% Data for fitting
E_p = E(Range_E); % selected probe energies
delta_Abs = -1*data(Range_E,t);
delta_Abs_norm = delta_Abs./max(abs(delta_Abs)); % normalised delta_Abs
Range_Efit = E_p>=1.62 & E_p<=max(E_p);
E_fit = E_p(Range_Efit);
delta_Abs_norm_fit = delta_Abs_norm(Range_Efit);
% Fitting function
function F = MB(y, E_fit)
F = y(1).*exp(-(E_fit./(8.617333268*10^-5.*y(2)))) + y(3);
end
%% Curve fitting options
% % Initial parameter guess and bounds
% lb = [0, 293, -1]; ub = [Inf, 800, 1];
% y0 = [1.2*10^9, 700, 0.5];
% lsqcurvefit and choose between different algorithm that lsqcurvefit employs (3C1, comment those lines that are not choosen and uncomment the line that is choosen, if not, matlab will take the last line of "optim_lsq" by default)
optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'MaxFunctionEvaluations',10^10, 'MaxIterations', 10^10, 'FunctionTolerance',10^-10, 'StepTolerance', 10^-10);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxFunctionEvaluations',10^10, 'MaxIterations',10^10, 'FunctionTolerance',10^-20, 'StepTolerance', 10^-20);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point', 'MaxFunctionEvaluations',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-20, 'StepTolerance', 10^-20);
ftns = @(y) norm(delta_Abs_norm_fit - MB(y, E_fit));
y0 = ga(ftns, 3, [], [], [], [], lb, ub)
% Solver for lsqcurvefit
[y, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@MB, y0, E_fit, delta_Abs_norm_fit, lb, ub, [], [], [], [], [], optim_lsq);
y
%% Plot command
plot(E_p, delta_Abs_norm,'Black')
hold on
plot(E_fit, MB(y, E_fit), 'LineWidth', 1.0, 'Color', 'red')
xlabel('Probe Photon Energy (eV)')
ylabel('Normalised \Delta A (a.u.)')
legend('Experimental Data', 'Fitted Curve')
disp(y(1,1))
disp(y(1,2))
disp(y(1,3))
