Error : Objective function is returning undefined values at initial point. lsqcurvefit cannot continue

I am having this ''Objective function is returning undefined values at initial point. lsqcurvefit cannot continue'' errror message from the code. Not sure what I should change. Please let me know your suggestions. I have attached the experimental data file as well.

clear all, close all,  clc
EQEACRXTN = xlsread('EQEACRXTN','Sheet1','A21:B53');% loading raw data from excel file
CurrentData = EQEACRXTN(:,1);
EQEData = EQEACRXTN(:,2);

%parameter estimation (units per ns)
B0 = [1e-10,1e-15]; %sequence [k_rs,k_ISC,k_RISC,k_NRT,k_SS,k_ST,k_TT] ]
lb = [1e-10,1e-15];% lower bound estimation
ub = [10e-10,5e-15];% upper bound estimation
options = optimset('Algorithm', 'trust-region-reflective');   
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@fitEQE,B0,CurrentData,EQEData,lb,ub,options);%sending data for fitting
F = fitEQE(B,CurrentData);% fitted data saving

figure(1)
plot(CurrentData,EQEData,'bo')%raw data
hold on
plot(CurrentData,F,'r','linewidth',2)%fit data 4.6
% grid on
% grid minor
xlabel('Current Density (mA/cm^2)');
ylabel('EQE (%)');
legend('Experimental EQE','EQE fit');
set(gca,'xscale','log')
set(gca,'yscale','log')

% 
function nEQE = fitEQE(B,CurrentData)
% % % NonLinear Part
syms Ns Nt 
assume(Ns > 0);
assume(Nt > 0);
k_rs = 1.1000e-02*1e9; 
k_ISC = 8.7545e-03*1e9; 
k_RISC = 9.9000e-04*1e9;
k_NRT = 1.9998e-04*1e9; 
d = 15e-7 ; 
e = 1.6e-19*1e3 ; 
J = CurrentData.';

sol_Ns = double((nan(numel(J),1)));
sol_Nt = double((nan(numel(J),1)));
for i = 1:numel(J)
eq1 = -(k_rs+k_ISC)*Ns+k_RISC*Nt-B(1)*Ns*Nt+0.25*B(2)*Nt.^2+(J(i)/(4*d*e))==0 ; 
eq2 = k_ISC*Ns-(k_RISC+k_NRT)*Nt-1.25*B(2)*Nt.^2+((3*J(i))/(4*d*e))==0 ; 
[sol_Ns(i),sol_Nt(i)] = vpasolve([eq1, eq2],[Ns,Nt]); 
end

%linear Part

syms Ns0 Nt0
assume(Ns0 > 0);
assume(Nt0 > 0);
%  nlinear1 =0;
%  nlinear2 = 0;
Singlet = double((nan(numel(J),1)));
Triplet = double((nan(numel(J),1)));
for ii = 1:numel(J)
    nlinear1 =-(k_rs+k_ISC)*Ns0+k_RISC*Nt0+(J(ii)/(4*d*e));
    nlinear2 = k_ISC*Ns0-(k_RISC+k_NRT)*Nt0+((3*J(ii))/(4*d*e));
    [Singlet(ii),Triplet(ii)] = solve([nlinear1,nlinear2],[Ns0,Nt0]);
end 

% % % EQE predition
nEQE0 = 3.8;
nEQE = nEQE0*(sol_Ns./Singlet);
end

Answer 1

[Walter Roberson]

You have a couple of issues:

1. The last entry of your second column of your file is empty, so you are getting NaN values. This is your most immediate problem.
2. vpasolve() is ignoring the implicit restriction of Ns and Nt to real numbers. Some of the results it is producing are complex valued. If you try to use vpasolve() with range restriction from eps(realmin) to infinity then vpasolve() will not be able to find solutions. The equations turn out to be cubic; if you use solve() then even though it is apparent that there should be no imaginary part for one of the three solutions, when you use vpa() or double() then there is a small imaginary residue due to round-off error.

You can manipulate the solution to the form

Nt = root() of a cubic
Ns = B2 * Nt^2/constant1 + Nt*constant2 - J * constant3

and you can solve the cubic and take the first root (which will typically be the real root), and you can take real() of it to eliminate the noise complex part, and then calculate Ns from there.

Or you could store the results from vpasolve() into temporary variables and store real() of them, making the assumption that vpasolve() will end up choosing the one of the three cubic roots that is just noise from being real-valued.