CODE:

function parafit

clear all;

t=[0 2 4 8 24 48];

y=[0.69 0.645 0.635 0.62 0.61 0.61];

y0=0.69;

% Nonlinear least square estimate using lsqnonlin()

beta0=0.5;

lb=[0];ub=[inf];

[beta,resnorm,residual,exitflag,output,lambda,jacobian] = ...

lsqnonlin(@Func,beta0,lb,ub,[],t,y,y0);

ci = nlparci(beta,residual,jacobian);

beta;

% result

fprintf('\n Estimated Parameters by Lsqnonlin():\n')

fprintf('\t k = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))

fprintf('  The sum of the residual squares is: %.1e\n\n',sum(residual.^2))

% plot of fit results

tspan = [0  max(t)];

[tt yc] = ode45(@ModelEqs,tspan,y0,[],beta);

tc=linspace(0,max(t),200);

yca = spline(tt,yc,tc);

plot(t,y,'ro',tc,yca,'r-');

hold on

xlabel('Time');

ylabel('Concentration');

hold off

% =======================================

function f1 = Func(beta,t,y,y0)        % Define objective function

tspan =t;

[tt yy] = ode45(@ModelEqs,tspan,y0,[],beta);

yc= spline(tt,yy,t);

f1=y-yc;

% ==================================

function dydt = ModelEqs(t,y,beta)          % Model equations

dydt = -beta*(0.7-y).^(1/3)*(y-0.61);