Why lsqnonlin results worst than main simulation

I am optimizing a simulation of experimental data. The problem is that I have more than one series of data. I did the optimization using LSQNONLIN and it fitted the parameters quite well for one series of data. Now, I am trying to do the optimization for all series of data. The problem is that when I do so, I get worse results compared to the first simulation (before optimization) and the results are very far from optimized. My code is here. It doesn’t have any errors or so. I have no idea how to solve the problem. Can anybody help me?

 %%MAIN FUNCTION
 % Import constant data to the main function
 [Input,options] = input_data(mC,mH,mO);     % This function read all excel sheets and calculate some parameters which are constant for each iteration

 % Lower bound
 lb = 1e-30*ones(12,1);

 % Initial guess
 k0 = [6.2489,6.5724,2.1093,1.6735,6.8555,2.1517,3.4014,2.1483,5.7985,9.5707,4.427,1];

 % Options for "lsqnonlin"
 oldoptions = optimset('Display',          'iter'  ,...
     'Diagnostics',       'on'     ,...
     'Algorithm',         'trust-region-reflective'  ,...
     'FinDiffType',       'central'  ,...
     'FinDiffRelStep',    1e-03  ,...
     'MaxFunEvals',       5000   ,...
     'MaxIter',           200    ,...
     'TolX',              1e-16  ,...
     'TolFun',            1e-14);

 % Parameter estimation
 [k,resnorm,residual,exitflag,output] = lsqnonlin(@(k)PFR(k,Input,options,std_dev,MW_E),k0,lb,[],opts);

 ----------------------------------------------------------------------------------------------------------------------------------
 %%OBJECTIVE FUNCTION: This function calculates the difference between simulation concentrations and experimental concentrations for all excel sheets
function final = PFR(k,Input,options,std_dev,MW_E)

 final = [];
 for w = 1:length(Input.ndata)

    % Integration
    [z,f] = ode23tb(@(z,g)PFR_DGL(z,g,Input.n,Input.URH{w,1},k,Input.kTu{w,1},Input.n_u),Input.zspan{w,1},Input.f0{w,1},options);

    % Concentration of species which their experimental concentrations are exist
    final_f = f(:,[2:6]);

    % Experimental concentrations in mol/m³
    conc_E = Input.ndata{w,1}(:,9:12).*Input.rhoRH{w,1}*1000;

    % Difference between simulation concentrations and experimental concentrations
    conc_diff = final_f - final_E;

    % Transfer from matrix to a column vector
    conc_diff_to_vector = conc_diff(:);

    % Adding new conc_diff_to_vector (resulted from the new excel sheet) in the continue of last ones
    final = [final; conc_diff_to_vector];         % final would be a column vector

 end

 ----------------------------------------------------------------------------------------------------------------------------------
 %%THIS FUNCTION IS USED BEY ODE

 function f = PFR_DGL(z,g,n,URH,k,kTu,n_u)

 % Reaction rate function
 ru = rate_Po2_u(g(1:n),k,kTu);   % this function included the parameters which should be estimated

 % Output
 f(1:n) = ((n_u'*ru)/URH);

 -----------------------------------------------------------------------------------------------------------------------------------
 %%THIS FUNCTION INCLUDED THE PARAMETERS WHICH SHOULD BE ESTIMATED

 function ru = rate_Po2_u(c,k,kTu)

 % "k" shows estimated parameters; "kTu" shows the parameters which are calculated in "input_data" function
 ru = [ k(1)*c(1)*c(2)
     kTu(2)*c(3)*c(2)
     k(2)*c(1)*c(5)
     .
     .
     .
     kTu(20)*c(1)*c(7)
     k(11)*c(12)];

Thanks in advance,