Difference in computational time

Hi all, I am working on a curve fitting project and I have a couple of curves to fit. For 2 of the curves, the fitting function involves an integration and I was told that I will need to make each of the data loop through the function (using the for loop). The below is my first code, which takes less than a second to complete.

tic
% Preparation

clear; clc 

load Steady_State_Data.mat

%% Preamble

% Fundamental constants 
h = 4.0135667696*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units

% Clean up of data to select range of values 
Absorbance = log10(T_substrate./T_sample);
e = E >= 0 & E <= 2.0; % returns boolean value of the indices that satisfies the logical condition defined above
A = find(e); % gives indices that are non-zero
N = length(A); % no. of elements that are non-zero 
n = length(E) + 1; 

% Data for fitting
E_p = E(n-N:n-1);
Abs = Absorbance(n-N:n-1) - min(Absorbance);

% Fitting function: Elliott's Model for Steady State with non-parabolic factor 
function F = EM_SS_wR(p, e_p)
  for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*...
    (integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...
    126*(p(5))^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), p(3), Inf, 'ArrayValued', 1)) + ...
    p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(...
    (1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + ...
    (1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + ...
    (1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + ...
    (1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + ...
    (1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + ...
    (1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + ...
    (1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6))); 
  end
  F = F(:);
end

% Carrier Contribution
function F = EM_SS_wR_CC(p, e_p)
  for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*...
    (integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...
    126*(p(5))^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), p(3), Inf, 'ArrayValued', 1));
  end
  F = F(:);
end

% Excitoninc Contribution 
EM_SS_wR_EC = @(p, E_p) p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(...
    (1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + ...
    (1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + ...
    (1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + ...
    (1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + ...
    (1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + ...
    (1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + ...
    (1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6)));    



%% Curve fitting options

% Initial parameter guess and bounds
lb = [-Inf, -Inf, 0, 0, -Inf, -Inf]; ub = [Inf, Inf, 55, 0.030, Inf, Inf];
p0 = [0.13, 0.11, 1.4, 0.025, 0.12, 0.035]; % refer to the next line for their order
% p0 = [A1, A2, Eg, Eb, R, g]


% 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',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxFunctionEvaluations',1000, 'MaxIterations', 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point', 'MaxFunctionEvaluations',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);


% Solver for lsqcurvefit 
[p, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@EM_SS_wR, p0, E_p, Abs, lb, ub);


%% Error bars calculation 
ci = nlparci(p, residual, 'Jacobian', jacobian);
PCI = table(ci(:,1), p(:), ci(:,2),'VariableNames',{'CI 0.025','p','CI 0.975'});
Parameter_CI = table2array(PCI);
upper_bar = (Parameter_CI(:,3) - Parameter_CI(:,2))./2;
lower_bar = (Parameter_CI(:,2) - Parameter_CI(:,1))./2;


%% Parameter values (refer to command window)
p1 = p(1,1);
p2 = p(1,2);
p3 = p(1,3);
p4 = p(1,4);
p5 = p(1,5);
p6 = p(1,6);
e1 = lower_bar(1,1);
e2 = lower_bar(2,1);
e3 = lower_bar(3,1);
e4 = lower_bar(4,1);
e5 = lower_bar(5,1);
e6 = lower_bar(6,1);
X1 = [' A1 = ', num2str(p1), char(177), num2str(e1)];
X2 = [' A2 = ', num2str(p2), char(177), num2str(e2)];
X3 = [' Eg = ', num2str(p3), char(177), num2str(e3)];
X4 = [' Eb = ', num2str(p4), char(177), num2str(e4)];
X5 = [' R = ', num2str(p5), char(177), num2str(e5)];
X6 = [' g = ', num2str(p6), char(177), num2str(e6)];
disp(X1);
disp(X2);
disp(X3);
disp(X4);
disp(X5);
disp(X6);

% Table of parameter values
parameter_name = {'A1'; 'A2'; 'Eg'; 'Eb'; 'R'; 'g'};
parameter_values = [p1; p2; p3; p4; p5; p6];
error = [e1; e2; e3; e4; e5; e6];
T = table(parameter_values, error, 'RowNames',parameter_name);
disp(T)


%% Plot command

plot(E_p, Abs, 'o', 'Color','b') % scatter plot
hold on
plot(E_p, EM_SS_wR(p, E_p), 'LineWidth', 1.0, 'Color','red') % curve fitting 
plot(E_p, EM_SS_wR_CC(p, E_p), 'LineWidth', 1.0,'Color','black') % carrier contribution
plot(E_p, EM_SS_wR_EC(p, E_p), 'LineWidth', 1.0, 'Color','green') % excitonic contribution
% Table of parameter values
TString = evalc('disp(T)');
TString = strrep(TString,'<strong>','\bf');
TString = strrep(TString,'</strong>','\rm');
TString = strrep(TString,'_','\_');
FixedWidth = get(0,'FixedWidthFontName');
annotation(gcf,'Textbox','String',TString,'Interpreter','tex',...
           'FontName',FixedWidth,'Units','Normalized','Position',[0.15, 0.8, 0.1, 0.1]);
xlabel('Probe energy (eV)')
ylabel('Absorbance.O.D')
legend('Experimental Data', 'Fitted Curve', 'Carrier Contribution','Excitonic Contribution', 'Location','southeast')
hold off 


%% Miscellaneous

% For a rectangular pic: 
% drag horizontally until it just covers the letter 'c' in the word 'col'

% For a square pic: 
% drag horizontall until it just covers the letter 'o' in the word 'Contribution' in the next line)
% legend('Experimental Data', 'Fitted Curve', 'Carrier Contribution','Excitonic Contribution')
toc

However, when using my second code (shown below) it takes around 5 minutes to complete running. May I please ask for the difference in the computational time?

tic
%% Preparation 
clc; clear

data = importdata("FCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv"); % insert file path within parenthesis 
load Steady_State_Parameter_Values.mat

%% Preamble 

% Fundamental constants
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K

% Data 
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
delay_t = data(1, :);

% Parameter values from Steady state fiting 
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);

t1 = 0.5342;
t2 = 1.0257;
t3 = 1.9673;
t4 = 3.9554;

carrier_T = [1198.8, 816.7, 446.8, 328.7]; 

col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;


% Data for fitting 
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
delta_Abs1 = data(Range_W,col1);
delta_Abs2 = data(Range_W,col2);
delta_Abs3 = data(Range_W,col3);
delta_Abs4 = data(Range_W,col4);

% Fitting function: Elliott's Model for Transient Absorption
function F = EM_TA_wR1(x, e_p)
load Steady_State_Parameter_Values.mat
kB = 8.617333268*10^-5; % units: eV/ K
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6); 
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 1;
    for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
    126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
    A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
    (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
    126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
    x(6)*A2*(4*pi*(x(2))^3/2)./(E_p).*1/x(3)*(...
    (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
    (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
    (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
    (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
    (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
    (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
    (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
    A2*(4*pi*Eb^3/2)./E_p.*1/g*(...
    (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
    (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
    (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
    (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
    (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
    (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
    (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
    end
F = F(:);
end 

function F = EM_TA_wR2(x, e_p)
data = importdata("Experimental data\Transient Absorption\FCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv"); % insert file path within parenthesis 
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K 
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6); 
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 2;
    for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
    126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
    A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
    (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
    126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
    x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(...
    (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
    (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
    (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
    (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
    (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
    (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
    (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
    A2*(4*pi*Eb^3/2)*1/g*(...
    (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
    (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
    (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
    (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
    (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
    (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
    (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
    end
F = F(:);
end 

function F = EM_TA_wR3(x, e_p) 
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K 
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6); 
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 3;
    for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
    126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
    A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
    (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
    126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
    x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(...
    (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
    (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
    (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
    (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
    (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
    (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
    (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
    A2*(4*pi*Eb^3/2)*1/g*(...
    (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
    (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
    (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
    (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
    (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
    (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
    (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
    end
F = F(:);
end 

function F = EM_TA_wR4(x, e_p)
load Steady_State_Parameter_Values.mat
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K 
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6); 
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 4;
    for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
    126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
    A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
    (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
    126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
    x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(...
    (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
    (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
    (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
    (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
    (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
    (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
    (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
    A2*(4*pi*Eb^3/2)*1/g*(...
    (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
    (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
    (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
    (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
    (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
    (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
    (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
    end
F = F(:);
end 

% Curve fitting
lb = [Eg, Eb, g, 0.3, 0.5, 0.2]; ub = [55, 0.05, 0.05, 20, 2, 1];
x0 = [1.65, 0.03, 0.03, 1.3, 1, 0.3]; 
optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'MaxFunctionEvaluations',10^5, 'MaxIterations', 10^5, 'FunctionTolerance',10^-20, 'StepTolerance', 10^-20);
delta_Abs(:,1)= data_new2(:,1);
x1 = lsqcurvefit(@EM_TA_wR1, x0, E_p, delta_Abs1, lb, ub, optim_lsq);
x2 = lsqcurvefit(@EM_TA_wR2, x1, E_p, delta_Abs2, lb, ub, optim_lsq);
x3 = lsqcurvefit(@EM_TA_wR3, x2, E_p, delta_Abs3, lb, ub, optim_lsq);
x4 = lsqcurvefit(@EM_TA_wR4, x3, E_p, delta_Abs4, lb, ub, optim_lsq);


plot(E_p, delta_Abs1, 'o','color','blue')
hold on 
plot(E_p, delta_Abs2, 'o', 'color', 'yellow')
plot(E_p, delta_Abs3, 'o', 'color', 'green')
plot(E_p, delta_Abs4, 'o', 'color', 'magenta')
plot(E_p, EM_TA_wR1(x1, E_p), 'color', 'blue', 'LineWidth', 4.0)
plot(E_p, EM_TA_wR2(x2, E_p), 'color', 'yellow', 'LineWidth', 4.0)
plot(E_p, EM_TA_wR3(x3, E_p), 'color', 'green', 'LineWidth', 4.0)
plot(E_p, EM_TA_wR4(x4, E_p), 'color', 'magenta', 'LineWidth', 4.0)
xlabel('Probe Photon Energy (eV)')
ylabel('\Delta A (O.D.)')
legend('0.5 ps', '1.0 ps', '2.0 ps', '4.0 ps', 'Location', 'southeast')
hold off

disp(x1)
disp(x2)
disp(x3)
disp(x4)

toc

[Steven Lord]

What happens when you profile the two codes? What sections of each code take the most time? Is it the equivalent section in the two codes or is the hotspot different between the two codes?

Are you confident that both the codes are solving the same problem? If so, double check. If I had a dollar for every time I've written "the same code" two different ways only to realize that they weren't actually doing the same thing ... I would have a non-zero number of dollars :)

Answer 1

[Star Strider]

It appears to work correctly.

What else needs to be done? (W.R.T. ‘Miscellaneous’, you can use the daspect funciton or different axis options.)

tic
% Preparation

clear; clc 

load Steady_State_Data.mat

%% Preamble

% Fundamental constants 
h = 4.0135667696*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units

% Clean up of data to select range of values 
Absorbance = log10(T_substrate./T_sample);
e = E >= 0 & E <= 2.0; % returns boolean value of the indices that satisfies the logical condition defined above
A = find(e); % gives indices that are non-zero
N = length(A); % no. of elements that are non-zero 
n = length(E) + 1; 

% Data for fitting
E_p = E(n-N:n-1);
Abs = Absorbance(n-N:n-1) - min(Absorbance);

% Fitting function: Elliott's Model for Steady State with non-parabolic factor 
function F = EM_SS_wR(p, e_p)
  for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*...
    (integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...
    126*(p(5))^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), p(3), Inf, 'ArrayValued', 1)) + ...
    p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(...
    (1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + ...
    (1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + ...
    (1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + ...
    (1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + ...
    (1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + ...
    (1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + ...
    (1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6))); 
  end
  F = F(:);
end

% Carrier Contribution
function F = EM_SS_wR_CC(p, e_p)
  for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = p(1)*(2*pi*sqrt(p(4))/E_p)*(1/p(6))*...
    (integral(@(E)sech(((E_p - E)./p(6)))*(1 + 10*p(5)*(E - p(3)) + ...
    126*(p(5))^2*(E - p(3))^2)/(1 - exp(-2*pi*sqrt(p(4)/(E - p(3))))), p(3), Inf, 'ArrayValued', 1));
  end
  F = F(:);
end

% Excitoninc Contribution 
EM_SS_wR_EC = @(p, E_p) p(2)*((2*pi*p(4)^3/2)/E_p)*1/p(6)*(...
    (1/1^3)*sech((E_p - p(3) + p(4)/1^2)./p(6)) + ...
    (1/2^3)*sech((E_p - p(3) + p(4)/2^2)./p(6)) + ...
    (1/3^3)*sech((E_p - p(3) + p(4)/3^2)./p(6)) + ...
    (1/4^3)*sech((E_p - p(3) + p(4)/4^2)./p(6)) + ...
    (1/5^3)*sech((E_p - p(3) + p(4)/5^2)./p(6)) + ...
    (1/6^3)*sech((E_p - p(3) + p(4)/6^2)./p(6)) + ...
    (1/7^3)*sech((E_p - p(3) + p(4)/7^2)./p(6)));    



%% Curve fitting options

% Initial parameter guess and bounds
lb = [-Inf, -Inf, 0, 0, -Inf, -Inf]; ub = [Inf, Inf, 55, 0.030, Inf, Inf];
p0 = [0.13, 0.11, 1.4, 0.025, 0.12, 0.035]; % refer to the next line for their order
% p0 = [A1, A2, Eg, Eb, R, g]


% 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',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'trust-region-reflective', 'MaxFunctionEvaluations',1000, 'MaxIterations', 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'interior-point', 'MaxFunctionEvaluations',1000, 'MaxIterations', 1000, 'FunctionTolerance',10^-9, 'StepTolerance', 10^-9);


% Solver for lsqcurvefit 
[p, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@EM_SS_wR, p0, E_p, Abs, lb, ub);

Local minimum possible.

lsqcurvefit stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

%% Error bars calculation 
ci = nlparci(p, residual, 'Jacobian', jacobian);
PCI = table(ci(:,1), p(:), ci(:,2),'VariableNames',{'CI 0.025','p','CI 0.975'});
Parameter_CI = table2array(PCI);
upper_bar = (Parameter_CI(:,3) - Parameter_CI(:,2))./2;
lower_bar = (Parameter_CI(:,2) - Parameter_CI(:,1))./2;


%% Parameter values (refer to command window)
p1 = p(1,1);
p2 = p(1,2);
p3 = p(1,3);
p4 = p(1,4);
p5 = p(1,5);
p6 = p(1,6);
e1 = lower_bar(1,1);
e2 = lower_bar(2,1);
e3 = lower_bar(3,1);
e4 = lower_bar(4,1);
e5 = lower_bar(5,1);
e6 = lower_bar(6,1);
X1 = [' A1 = ', num2str(p1), char(177), num2str(e1)];
X2 = [' A2 = ', num2str(p2), char(177), num2str(e2)];
X3 = [' Eg = ', num2str(p3), char(177), num2str(e3)];
X4 = [' Eb = ', num2str(p4), char(177), num2str(e4)];
X5 = [' R = ', num2str(p5), char(177), num2str(e5)];
X6 = [' g = ', num2str(p6), char(177), num2str(e6)];
disp(X1);

 A1 = 0.041666±0.0026577

disp(X2);

 A2 = 36.9262±16.039

disp(X3);

 Eg = 1.6333±0.006288

disp(X4);

 Eb = 0.024755±0.004483

disp(X5);

 R = 0.22164±0.0071215

disp(X6);

 g = 0.024537±0.00067958

% Table of parameter values
parameter_name = {'A1'; 'A2'; 'Eg'; 'Eb'; 'R'; 'g'};
parameter_values = [p1; p2; p3; p4; p5; p6];
error = [e1; e2; e3; e4; e5; e6];
T = table(parameter_values, error, 'RowNames',parameter_name);
disp(T)

          parameter_values      error   
          ________________    __________

    A1        0.041666         0.0026577
    A2          36.926            16.039
    Eg          1.6333          0.006288
    Eb        0.024755          0.004483
    R          0.22164         0.0071215
    g         0.024537        0.00067958

%% Plot command

plot(E_p, Abs, 'o', 'Color','b') % scatter plot
hold on
plot(E_p, EM_SS_wR(p, E_p), 'LineWidth', 1.0, 'Color','red') % curve fitting 
plot(E_p, EM_SS_wR_CC(p, E_p), 'LineWidth', 1.0,'Color','black') % carrier contribution
plot(E_p, EM_SS_wR_EC(p, E_p), 'LineWidth', 1.0, 'Color','green') % excitonic contribution
% Table of parameter values
TString = evalc('disp(T)');
TString = strrep(TString,'<strong>','\bf');
TString = strrep(TString,'</strong>','\rm');
TString = strrep(TString,'_','\_');
FixedWidth = get(0,'FixedWidthFontName');
annotation(gcf,'Textbox','String',TString,'Interpreter','tex',...
           'FontName',FixedWidth,'Units','Normalized','Position',[0.15, 0.8, 0.1, 0.1]);
xlabel('Probe energy (eV)')
ylabel('Absorbance.O.D')
legend('Experimental Data', 'Fitted Curve', 'Carrier Contribution','Excitonic Contribution', 'Location','southeast')
hold off 

%% Miscellaneous

% For a rectangular pic: 
% drag horizontally until it just covers the letter 'c' in the word 'col'

% For a square pic: 
% drag horizontall until it just covers the letter 'o' in the word 'Contribution' in the next line)
% legend('Experimental Data', 'Fitted Curve', 'Carrier Contribution','Excitonic Contribution')
toc

Elapsed time is 3.191509 seconds.

.

[Jack]

Hi Star, thanks for the tips on "Miscellaneous"! The second code just takes way too long as comapred to the first one (the one that you've ran), and I am wondering what is the thing that is taking up so much time.

Also, is there a more efficient way for the second code (i.e. not calling the function 4 time)? Because each y-data has a associated fitting function (these fitting functions are the same, but that they differ in the value of a constant named carrier_T in my code) and I am wondering how can I code it to avoid having to write out the same function 4 times for each value of the constant. Thank you!

[Star Strider]

My pleasure!

I’m slightly mystified. What code versions do you want to compare? Are they posted here?

[Jack]

Hi Star, the two versions of the code are posted in my original post. Please do let me know if you have trouble seeing them!

one of them is at the star of the post, and the other is seperated by a few lines of words

[Torsten]

The functions "EM_TA_wRx" are usually called thousands of times from "lsqcurvefit". So don't load the same data in these functions again and again.

Load them once before you call "lsqcurvefit" for the first time and pass them to the functions afterwards.

[Jack]

I see. Do I have to define a global variable such that the function recognises the variables outside of it? The reason for me loading the data into the function is due to the fact that I am unable to get the function to recognise the variables I need it to. May I ask what I should do other than loading the data into these functions? Thank you!

[Torsten]

Z = load(...)
[p, residualnorm, residual, exitflag, output, lambda, jacobian] = lsqcurvefit(@(p,E_p)EM_SS_wRx(p,E_p,Z),...)
    
function f = EM_SS_wRx(p,E_p,Z)
...
end

[Star Strider]

Correct. I didn’t catch that initially because I didn’t copy it and run it (thhinking everything was captured in the same ‘Copy’ operattion.

tic
%% Preparation 
clc; clear

data = importdata("FCPIB-293K-2.5...bg -chirp.csv"); % insert file path within parenthesis 
load('Steady_State_P...er_Values.mat') 

%% Preamble 

% Fundamental constants
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K

% Data 
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
delay_t = data(1, :);

% Parameter values from Steady state fiting 
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6);

t1 = 0.5342;
t2 = 1.0257;
t3 = 1.9673;
t4 = 3.9554;

carrier_T = [1198.8, 816.7, 446.8, 328.7]; 

col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;


% Data for fitting 
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
delta_Abs1 = data(Range_W,col1);
delta_Abs2 = data(Range_W,col2);
delta_Abs3 = data(Range_W,col3);
delta_Abs4 = data(Range_W,col4);

% Fitting function: Elliott's Model for Transient Absorption
function F = EM_TA_wR1(x, e_p)
load('Steady_State_P...er_Values.mat') 
kB = 8.617333268*10^-5; % units: eV/ K
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6); 
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 1;
    for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
    126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
    A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
    (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
    126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
    x(6)*A2*(4*pi*(x(2))^3/2)./(E_p).*1/x(3)*(...
    (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
    (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
    (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
    (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
    (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
    (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
    (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
    A2*(4*pi*Eb^3/2)./E_p.*1/g*(...
    (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
    (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
    (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
    (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
    (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
    (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
    (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
    end
F = F(:);
end 

function F = EM_TA_wR2(x, e_p)
data = importdata("FCPIB-293K-2.5...bg -chirp.csv"); % insert file path within parenthesis 
load('Steady_State_P...er_Values.mat') 
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K 
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6); 
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 2;
    for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
    126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
    A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
    (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
    126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
    x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(...
    (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
    (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
    (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
    (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
    (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
    (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
    (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
    A2*(4*pi*Eb^3/2)*1/g*(...
    (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
    (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
    (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
    (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
    (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
    (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
    (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
    end
F = F(:);
end 

function F = EM_TA_wR3(x, e_p) 
data = importdata("FCPIB-293K-2.5...bg -chirp.csv"); % insert file path within parenthesis 
load('Steady_State_P...er_Values.mat') 
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K 
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6); 
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 3;
    for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
    126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
    A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
    (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
    126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
    x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(...
    (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
    (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
    (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
    (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
    (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
    (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
    (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
    A2*(4*pi*Eb^3/2)*1/g*(...
    (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
    (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
    (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
    (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
    (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
    (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
    (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
    end
F = F(:);
end 

function F = EM_TA_wR4(x, e_p)
data = importdata("FCPIB-293K-2.5...bg -chirp.csv"); % insert file path within parenthesis 
load('Steady_State_P...er_Values.mat') 
h = 4.1356677*10^-15; % units: eV/ Hz
c = 3*10^8; % SI units
kB = 8.617333268*10^-5; % units: eV/ K 
Wavelength = data(:, 1);% units: nm
E = (h*c)./(Wavelength*10^-9);
A1 = p(1,1);
A2 = p(1,2);
Eg = p(1,3);
Eb = p(1,4);
R = p(1,5);
g = p(1,6); 
col1 = 56;
col2 = 63;
col3 = 74;
col4 = 87;
Range_E = E >= 1.5 & E <= 2.0;
Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
E_p = E(Range_E); % selected probe energies
data_new2 = data(Range_W, [col1,col2,col3,col4]);
carrier_T = [1198.8, 816.7, 446.8, 328.7];
n = 4;
    for i = 1:numel(e_p)
    E_p = e_p(i);
    F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
    126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
    A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
    (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
    126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
    x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(...
    (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
    (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
    (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
    (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
    (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
    (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
    (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
    A2*(4*pi*Eb^3/2)*1/g*(...
    (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
    (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
    (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
    (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
    (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
    (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
    (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
    end
F = F(:);
end 

% Curve fitting
lb = [Eg, Eb, g, 0.3, 0.5, 0.2]; ub = [55, 0.05, 0.05, 20, 2, 1];
x0 = [1.65, 0.03, 0.03, 1.3, 1, 0.3]; 
optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'MaxFunctionEvaluations',10^5, 'MaxIterations', 10^5, 'FunctionTolerance',10^-20, 'StepTolerance', 10^-20);
delta_Abs(:,1)= data_new2(:,1);
x1 = lsqcurvefit(@EM_TA_wR1, x0, E_p, delta_Abs1, lb, ub, optim_lsq);

Local minimum possible.
lsqcurvefit stopped because the relative size of the current step is less than
the value of the step size tolerance.

x2 = lsqcurvefit(@EM_TA_wR2, x1, E_p, delta_Abs2, lb, ub, optim_lsq);

Local minimum possible.
lsqcurvefit stopped because the relative size of the current step is less than
the value of the step size tolerance.

x3 = lsqcurvefit(@EM_TA_wR3, x2, E_p, delta_Abs3, lb, ub, optim_lsq);

Local minimum possible.
lsqcurvefit stopped because the relative size of the current step is less than
the value of the step size tolerance.

x4 = lsqcurvefit(@EM_TA_wR4, x3, E_p, delta_Abs4, lb, ub, optim_lsq);

Local minimum possible.
lsqcurvefit stopped because the relative size of the current step is less than
the value of the step size tolerance.

plot(E_p, delta_Abs1, 'o','color','blue')
hold on 
plot(E_p, delta_Abs2, 'o', 'color', 'yellow')
plot(E_p, delta_Abs3, 'o', 'color', 'green')
plot(E_p, delta_Abs4, 'o', 'color', 'magenta')
plot(E_p, EM_TA_wR1(x1, E_p), 'color', 'blue', 'LineWidth', 4.0)
plot(E_p, EM_TA_wR2(x2, E_p), 'color', 'yellow', 'LineWidth', 4.0)
plot(E_p, EM_TA_wR3(x3, E_p), 'color', 'green', 'LineWidth', 4.0)
plot(E_p, EM_TA_wR4(x4, E_p), 'color', 'magenta', 'LineWidth', 4.0)
xlabel('Probe Photon Energy (eV)')
ylabel('\Delta A (O.D.)')
legend('0.5 ps', '1.0 ps', '2.0 ps', '4.0 ps', 'Location', 'southeast')
hold off

disp(x1)

    1.6506    0.0378    0.0397    1.4603    0.9383    0.3780

disp(x2)

    1.6533    0.0405    0.0301    1.5501    0.9295    0.2000

disp(x3)

    1.6544    0.0347    0.0251    1.5752    1.0028    0.2000

disp(x4)

    1.6551    0.0345    0.0251    1.5903    1.0065    0.2000

toc

Elapsed time is 153.813683 seconds.

Loading the files and then passing their outputs as an additional parameter —

% tic
% %% Preparation 
% clc; clear
% 
% data = importdata("FCPIB-293K-2.5...bg -chirp.csv"); % insert file path within parenthesis 
% LD = load('Steady_State_P...er_Values.mat')
% p = LD.p
% 
% %% Preamble 
% 
% % Fundamental constants
% h = 4.1356677*10^-15; % units: eV/ Hz
% c = 3*10^8; % SI units
% kB = 8.617333268*10^-5; % units: eV/ K
% 
% % Data 
% Wavelength = data(:, 1);% units: nm
% E = (h*c)./(Wavelength*10^-9);
% delay_t = data(1, :);
% 
% % Parameter values from Steady state fiting 
% A1 = p(1,1);
% A2 = p(1,2);
% Eg = p(1,3);
% Eb = p(1,4);
% R = p(1,5);
% g = p(1,6);
% 
% t1 = 0.5342;
% t2 = 1.0257;
% t3 = 1.9673;
% t4 = 3.9554;
% 
% carrier_T = [1198.8, 816.7, 446.8, 328.7]; 
% 
% col1 = 56;
% col2 = 63;
% col3 = 74;
% col4 = 87;
% 
% 
% % Data for fitting 
% Range_E = E >= 1.5 & E <= 2.0;
% Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
% E_p = E(Range_E); % selected probe energies
% data_new2 = data(Range_W, [col1,col2,col3,col4]);
% delta_Abs1 = data(Range_W,col1);
% delta_Abs2 = data(Range_W,col2);
% delta_Abs3 = data(Range_W,col3);
% delta_Abs4 = data(Range_W,col4);
% 
% % Fitting function: Elliott's Model for Transient Absorption
% function F = EM_TA_wR1(x, e_p, p)
% % load Steady_State_Parameter_Values.mat
% kB = 8.617333268*10^-5; % units: eV/ K
% A1 = p(1,1);
% A2 = p(1,2);
% Eg = p(1,3);
% Eb = p(1,4);
% R = p(1,5);
% g = p(1,6); 
% carrier_T = [1198.8, 816.7, 446.8, 328.7];
% n = 1;
%     for i = 1:numel(e_p)
%     E_p = e_p(i);
%     F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
%     126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
%     A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
%     (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
%     126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
%     x(6)*A2*(4*pi*(x(2))^3/2)./(E_p).*1/x(3)*(...
%     (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
%     (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
%     (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
%     (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
%     (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
%     (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
%     (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
%     A2*(4*pi*Eb^3/2)./E_p.*1/g*(...
%     (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
%     (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
%     (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
%     (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
%     (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
%     (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
%     (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
%     end
% F = F(:);
% end 
% 
% function F = EM_TA_wR2(x, e_p, p, data)
% % data = importdata("Experimental data\Transient Absorption\FCPIB-293K-2.5mW-400nm-Jan072021 -ibg -bg -chirp.csv"); % insert file path within parenthesis 
% % load Steady_State_Parameter_Values.mat
% h = 4.1356677*10^-15; % units: eV/ Hz
% c = 3*10^8; % SI units
% kB = 8.617333268*10^-5; % units: eV/ K 
% Wavelength = data(:, 1);% units: nm
% E = (h*c)./(Wavelength*10^-9);
% A1 = p(1,1);
% A2 = p(1,2);
% Eg = p(1,3);
% Eb = p(1,4);
% R = p(1,5);
% g = p(1,6); 
% col1 = 56;
% col2 = 63;
% col3 = 74;
% col4 = 87;
% Range_E = E >= 1.5 & E <= 2.0;
% Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
% E_p = E(Range_E); % selected probe energies
% data_new2 = data(Range_W, [col1,col2,col3,col4]);
% carrier_T = [1198.8, 816.7, 446.8, 328.7];
% n = 2;
%     for i = 1:numel(e_p)
%     E_p = e_p(i);
%     F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
%     126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
%     A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
%     (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
%     126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
%     x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(...
%     (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
%     (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
%     (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
%     (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
%     (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
%     (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
%     (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
%     A2*(4*pi*Eb^3/2)*1/g*(...
%     (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
%     (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
%     (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
%     (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
%     (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
%     (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
%     (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
%     end
% F = F(:);
% end 
% 
% function F = EM_TA_wR3(x, e_p, p, data) 
% % load Steady_State_Parameter_Values.mat
% h = 4.1356677*10^-15; % units: eV/ Hz
% c = 3*10^8; % SI units
% kB = 8.617333268*10^-5; % units: eV/ K 
% Wavelength = data(:, 1);% units: nm
% E = (h*c)./(Wavelength*10^-9);
% A1 = p(1,1);
% A2 = p(1,2);
% Eg = p(1,3);
% Eb = p(1,4);
% R = p(1,5);
% g = p(1,6); 
% col1 = 56;
% col2 = 63;
% col3 = 74;
% col4 = 87;
% Range_E = E >= 1.5 & E <= 2.0;
% Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
% E_p = E(Range_E); % selected probe energies
% data_new2 = data(Range_W, [col1,col2,col3,col4]);
% carrier_T = [1198.8, 816.7, 446.8, 328.7];
% n = 3;
%     for i = 1:numel(e_p)
%     E_p = e_p(i);
%     F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
%     126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
%     A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
%     (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
%     126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
%     x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(...
%     (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
%     (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
%     (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
%     (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
%     (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
%     (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
%     (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
%     A2*(4*pi*Eb^3/2)*1/g*(...
%     (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
%     (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
%     (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
%     (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
%     (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
%     (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
%     (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
%     end
% F = F(:);
% end 
% 
% function F = EM_TA_wR4(x, e_p, p, data)
% % load Steady_State_Parameter_Values.mat
% h = 4.1356677*10^-15; % units: eV/ Hz
% c = 3*10^8; % SI units
% kB = 8.617333268*10^-5; % units: eV/ K 
% Wavelength = data(:, 1);% units: nm
% E = (h*c)./(Wavelength*10^-9);
% A1 = p(1,1);
% A2 = p(1,2);
% Eg = p(1,3);
% Eb = p(1,4);
% R = p(1,5);
% g = p(1,6); 
% col1 = 56;
% col2 = 63;
% col3 = 74;
% col4 = 87;
% Range_E = E >= 1.5 & E <= 2.0;
% Range_W = Wavelength >= (h*c)/(2.0*10^-9) & Wavelength <= (h*c)/(1.5*10^-9);
% E_p = E(Range_E); % selected probe energies
% data_new2 = data(Range_W, [col1,col2,col3,col4]);
% carrier_T = [1198.8, 816.7, 446.8, 328.7];
% n = 4;
%     for i = 1:numel(e_p)
%     E_p = e_p(i);
%     F(i) = x(5)*A1*((2*pi*sqrt(x(2)))/E_p)*1/g*(integral(@(E)sech(((E_p - E)./g))*(1 + 10*R*(E - x(1)) + ...
%     126*(R)^2*(E - x(1))^2)/(1 - exp(-2*pi*sqrt(x(2)/(E - x(1))))), x(1), Inf, 'ArrayValued', 1))*(1 - 1/(1 + exp((E_p - x(4))/(kB*carrier_T(1, n)))))^2 - ...
%     A1*(2*pi*sqrt(Eb)/E_p)*(1/g)*...
%     (integral(@(e)sech(((E_p - e)./g))*(1 + 10*R*(e - Eg) + ...
%     126*(R)^2*(e - Eg)^2)/(1 - exp(-2*pi*sqrt(Eb/(e - Eg)))), Eg, Inf, 'ArrayValued', 1)) + ...
%     x(6)*A2*(4*pi*(x(2))^3/2)*1/x(3)*(...
%     (1/1^3)*sech((E_p - x(1) + x(2)/1^2)./x(3)) + ...
%     (1/2^3)*sech((E_p - x(1) + x(2)/2^2)./x(3)) + ...
%     (1/3^3)*sech((E_p - x(1) + x(2)/3^2)./x(3)) + ...
%     (1/4^3)*sech((E_p - x(1) + x(2)/4^2)./x(3)) + ...
%     (1/5^3)*sech((E_p - x(1) + x(2)/5^2)./x(3)) + ...
%     (1/6^3)*sech((E_p - x(1) + x(2)/6^2)./x(3)) + ...
%     (1/7^3)*sech((E_p - x(1) + x(2)/7^2)./x(3))) - ...
%     A2*(4*pi*Eb^3/2)*1/g*(...
%     (1/1^3)*sech((E_p - Eg + Eb/1^2)./g) + ...
%     (1/2^3)*sech((E_p - Eg + Eb/2^2)./g) + ...
%     (1/3^3)*sech((E_p - Eg + Eb/3^2)./g) + ...
%     (1/4^3)*sech((E_p - Eg + Eb/4^2)./g) + ...
%     (1/5^3)*sech((E_p - Eg + Eb/5^2)./g) + ...
%     (1/6^3)*sech((E_p - Eg + Eb/6^2)./g) + ...
%     (1/7^3)*sech((E_p - Eg + Eb/7^2)./g)); 
%     end
% F = F(:);
% end 
% 
% % Curve fitting
% lb = [Eg, Eb, g, 0.3, 0.5, 0.2]; ub = [55, 0.05, 0.05, 20, 2, 1];
% x0 = [1.65, 0.03, 0.03, 1.3, 1, 0.3]; 
% optim_lsq = optimoptions('lsqcurvefit', 'Algorithm', 'levenberg-marquardt', 'MaxFunctionEvaluations',10^5, 'MaxIterations', 10^5, 'FunctionTolerance',10^-20, 'StepTolerance', 10^-20);
% delta_Abs(:,1)= data_new2(:,1);
% x1 = lsqcurvefit(@(x, e_p)EM_TA_wR1(x, e_p, p), x0, E_p, delta_Abs1, lb, ub, optim_lsq);
% x2 = lsqcurvefit(@(x, e_p)EM_TA_wR2(x, e_p, p, data), x1, E_p, delta_Abs2, lb, ub, optim_lsq);
% x3 = lsqcurvefit(@(x, e_p)EM_TA_wR3(x, e_p, p, data), x2, E_p, delta_Abs3, lb, ub, optim_lsq);
% x4 = lsqcurvefit(@(x, e_p)EM_TA_wR4(x, e_p, p, data), x3, E_p, delta_Abs4, lb, ub, optim_lsq);
% 
% 
% plot(E_p, delta_Abs1, 'o','color','blue')
% hold on 
% plot(E_p, delta_Abs2, 'o', 'color', 'yellow')
% plot(E_p, delta_Abs3, 'o', 'color', 'green')
% plot(E_p, delta_Abs4, 'o', 'color', 'magenta')
% plot(E_p, EM_TA_wR1(x1, E_p, p), 'color', 'blue', 'LineWidth', 4.0)
% plot(E_p, EM_TA_wR2(x2, E_p, p, data), 'color', 'yellow', 'LineWidth', 4.0)
% plot(E_p, EM_TA_wR3(x3, E_p, p, data), 'color', 'green', 'LineWidth', 4.0)
% plot(E_p, EM_TA_wR4(x4, E_p, p, data), 'color', 'magenta', 'LineWidth', 4.0)
% xlabel('Probe Photon Energy (eV)')
% ylabel('\Delta A (O.D.)')
% legend('0.5 ps', '1.0 ps', '2.0 ps', '4.0 ps', 'Location', 'southeast')
% hold off

% disp(x1)
% disp(x2)
% disp(x3)
% disp(x4)
% 
% toc

These may vary slightly between runs —

The first version (reading both files in each call to it) requires 153.813613 seconds.

The second version (passing the file outputs as additional parameters) requires 96.784151 seconds. (This version is currently commented-out.)

I had to comment-out the version not being run because MATLAB picked up on the additional definittions and threw an error.

.

[Jack]

Hi Star, thanks for the help! May I ask if there's an alternative method that I may use such that I could reudce the computation time to within 10 seconds?

[Star Strider]

My pleasure!

I looked through your code and I can’t see anything that could be improved through vectorisation or anything else. It just takes time to do those integrations.

(My apologies for the delay. The idiots at Micro$oft are stopping my desktop computer every three minutes and freezing if for a minute afterwards. They then just crash it a few times a day, apparenttly because this is fun for them. Micro$oft should simply be abolished, so that some group can create a good, reliable operating system without all the bloatware and otther <expleitives>.)

.

[Jack]

No worries Star (and XD on the comment on MS)! Understood on the time taken to do the integrations. May I ask, why is it that both the codes does integration but their computational time differs by so much? Any takes on that?

Appreciate your reply even facing technical difficulties!

[Torsten]

For the second case, you have 345 data in E_p instead of 74 for the first case, and you call lsqcurvefit 4 times instead of only once. Maybe the fitting is also more complicated than in the first case - check the number of calls to your 4 functions and compare with the first case which ends up within 3 seconds.

[Star Strider]

Thank you!

As something of a disclaimer, I am not certain what you’re doing. It’s not necessary that I understand, however my lack of knowledge limits my ability to suggest different approaches.

The codes are significantly different, and appear to do different calculations and produce differentt results, as seen in the resulting figures. The first one is much shorter.

In the second section, in the original code, the files are read each time. In my revision of the second code, the files are read once and the information in the files are passed as additional arguments to the functions. That results in about a minute of reduced execution time.

Not all functions have analytic integrals, so it may not be possible to do a symbolic integration of the function you’re integrating, create a function from it, and then plug the variables into that function, avoiding the numerical integration. That’s the only speed improvement I can think of.

[Jack]

I see. Thanks for all the input guys! :)

[Star Strider]

My (our) pleasure!

My apologies for the delay in responding. I had complete ISP failure from 19:00 Z 9 Jul 2024 to 02:00 Z 12 Jul 2024.