function fit_nonl
clear all;clc
format long

data=[    8.6667    2.0000
   10.0000    5.0000
   10.6667   10.0000
   12.6667   15.0000
   13.3333   30.0000
   15.0000   60.0000
   15.3333  120.0000
   16.6667  200.0000
   17.3333  300.0000
   18.3333  550.0000
   20.0000  700.0000
];
xspan=data(:,2);    %xµÄÊý¾Ý£¬ÔÚ´ËÊäÈë
texp=data(:,1);    %tµÄÊý¾Ý£¬ÔÚ´ËÊäÈë

k0=[0.1 0.1 1 100];
lb=-[1 1 1 1]*1e9;
ub=[1 1 1 1]*1e9;


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

% ʹÓú¯Êýlsqnonlin()½øÐвÎÊý¹À¼Æ

options=optimset('maxfunevals',1000);
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@objfunc,k0,lb,ub,options,xspan,texp);

ci = nlparci(k,residual,jacobian);
%residual;
fprintf('\n\nÄâºÏ½á¹û:\n')
fprintf('\n\t²ÎÊý e1 = %.16f',k(1))
fprintf('\n\t²ÎÊý e2 = %.16f',k(2))
fprintf('\n\t²ÎÊý e3 = %.16f',k(3))
fprintf('\n\t²ÎÊý t = %.16f',k(4))
y=kineticseqs(xspan,k);
r2=1-sum((texp-y).^2)./sum((texp-mean(y)).^2);
fprintf('\n\tr^2 = %.16f',r2);
figure
plot(xspan,kineticseqs(xspan,k),'b',xspan,texp,'or'),legend('¼ÆËãÖµ','ʵÑéÖµ','location','best')


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

function f = objfunc(k,xspan,texp)
f=kineticseqs(xspan,k)-texp;

%------------------------------------------------------------------------
function xt = kineticseqs(x,k)
a=2;
e1=k(1);e2=k(2);e3=k(3);t=k(4);
xt=a/e1 +(a/e2)*(1-exp(-x/t))+a*x/e3;