Parameters b1, b2[0,1], b3[1,2];
Variable a, x, y, z;
Constant R=1.073533812;
Function R*a=(1+exp(b1*(b2-b3)*z*(x-R*y)))^(-1)*b2*x+(1-(1+exp(b1*(b2-b3)*z*(x-R*y)))^(-1))*b3*x;
data;

function fit_nonl
clear all;clc
format long

data=[
];
R=1.073533812;
Xt_1=data(:,2);
Xt_2=data(:,3);
Xt_3=data(:,4);
Xt=data(:,1);
X=[Xt_1, Xt_2, Xt_3];
Y=R*Xt;
b0=[-10; 0.5; 1.5];
lb=[-100; 0; 1]*1e9;
ub=[0; 1; 2]*1e9;


%-------------------------------------------------------------------------

% ʹÓú¯Êýlsqnonlin()½øÐвÎÊý¹À¼Æ
OPTIONS=optimset('MaxFunEvals',1000);
[b,resnorm,residual,jacobian] = ...
    lsqnonlin(@ObjFunc,b0,lb,ub,OPTIONS,X,Y);
ci = nlparci(b,residual,jacobian);
%residual;

%-------------------------------------------------------------------------

function f = ObjFunc(b,X,Y)
f=KineticsEqs(X,b)-Y;

%------------------------------------------------------------------------
function xt = KineticsEqs(b,X)
b1=b(2);
b2=b(3);
x1=X(:,1);
xt=(myfunc)^(-1)*b1*x1+(1-(myfunc)^(-1))*b2*x1;

function n=myfunc(b,X)
R=1.073533812;
B=b(1);
b1=b(2);
b2=b(3);
n=1+exp(B*(b1-b2)*X(:,3)*(X(:,2)-R*X(:,1)));