estimating parameter for an ode system using lsqcurve

Hi,

i am trying to estimate two parameters (kr and delta) of an ode system by fitting it to experiemtal data (yield_ref) however i do not get a good fit when using lsqcurfit. I have tried to change the values of my initial guesses of the parameters to be estimated, changed the algorithm in the lsqcurvefit function and increased the tolerance and yet the obtained results are not quite the same as the experimental one. So i don't know if i have a problem in my code or it is just the case of the problem which gives such incoherent results. If you can suggest any other method of solving this or find an error in my way of coding i will be extremely grateful.

yield_ref=[5.32 5.11 4.42 3.46];    % represent the experimental yield for the highest velocity couple
C0=[211.8 116.89 85.98 42.47];      % vector of initial concentrations (mol/m^3)
kr=1;                               % initial guess of the rate constant value (m^3/mol.s)
delta=1e-9;                         % initial guess of the diffusion film thickness
P0=[kr delta];
lb=[0 1e-9];
ub=[100 1e-5];
options=optimoptions(@lsqcurvefit,'Algorithm','levenberg-marquardt','StepTolerance',1e-14,'OptimalityTolerance',1e-10);
[P,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@find_yield,P0,C0,yield_ref,lb,ub,options);
yield=find_yield(P,C0);
plot(C0,yield_ref,"Marker","+",'Color','r');
hold on 
plot(C0,yield,"Marker","+",'Color','b');

the function used in the code are:

function yield=find_yield(P,C0)

% This function code the ode system consisting of the following equations
%       dA/dz=-((kr*B0)/(u_aq*K_DB))*A
%       dCorgi/dz=K_DC*((kr*B0)/(u_aq*K_DB))*A
%       dC/dz=(Dcorg/(L-delta)*delta*u_org)*(Corgi-C)
%
% All parameters are fixed constants except for kr and delta which should
% be estimated
% P is a parameter vector where P(1) is kr and P(2) is delta 
% C0 is a vector of initial concentrations

vaq=1.264;                          % the mean value of the highest velocity in the aqueous phase
vorg=0.797;                         % the mean value of the highest velocity in the organic phase
KD_B=1900;
Dc_org=3e-10;
B0=1.125e3;
L=0.0125;
KD_C=[9.5 15.4 21.1 21.6]; 
n=length(C0);
yield=zeros(1,n);
zspan=linspace(0,L,20);

for i=1:n
    q0=[C0(i);0;0];
    [z1,q1]=ode45(@(z,q)internal_f(z,q,P,vaq,vorg,L,KD_B,Dc_org,B0,KD_C(i)),zspan,q0);
    Corg_out=q1(length(zspan),3);
    yield(i)=(Corg_out/C0(i))*100;
   
end
    
    function dydt=internal_f(z,q,P,vaq,vorg,L,KD_B,Dc_org,B0,KD_C)
                
        C1=(P(1)*B0)/(vaq*KD_B);
        C2=Dc_org/((L-P(2))*P(2)*vorg);
        
        dydt=[-C1*q(1);KD_C*C1*q(1);C2*(q(2)-q(3))];
        
    end

end

Thank you in advance for any help !