function Example_8
clear, clc, format short g, format compact
global Feval SumSave
load P6_9.txt  % Load experimental data as a text file. 
                                % First column - independent variable. 2nd
                                % and higher columns dependent variables.
                                % All columns must be full. Missing values
                                % are not permitted
Beta=[0.1 4 0.02 0.02]; % Initial estimate for parameter values
Feval=0;
SumSave=realmax;
[BetaOpt,Fval,exitflag,Output]=fminsearch(@SumSqr,Beta,[],P6_9);
disp(' ');
disp('Optimal Results ');
disp(' Parm. No.  Value')
for i=1:size(Beta,2)
    disp([i BetaOpt(i)]);
end
disp([' Final Sum of Squares ' num2str(Fval) '  Function evaluations  ' num2str(Feval)]);
tspan = [P6_9(1,1) P6_9(end,1)];
y0=[P6_9(1,2) P6_9(2,2)];
[t,y]=ode45(@ODEfun,tspan,y0,[],BetaOpt);
for i=1:size(y,2)
    hold off
    plot(t,y(:,i));
    title([' Plot of dependent variable y' int2str(i)]);
    xlabel(' Independent variable (t)');
    ylabel([' Dependent variable y' int2str(i)]);
    hold on
    plot(P6_9(:,1),P6_9(:,i+1),'*')
    pause
 end
function sum=SumSqr(Beta,data)
global Feval SumSave
Feval=Feval+1;
time=data(:,1);
ym=data(:,2:end);
yc(1,:)=ym(1,:);
[nr,nc]=size(ym); % nr - number of data points; 
                  % nc - number of dependent variables
tspan = [time(1) time(2)]; % Range for the independent variable 
y0 = ym(1,:); % Initial values for the dependent variables 
for i=2:nr
    [t,y]=ode45(@ODEfun,tspan,y0,[],Beta);
    yc(i,:)=y(end,:);
    if i<nr
        tspan = [time(i) time(i+1)];
        y0 = yc(i,:);
    end
end
sum=0;
for ic=1:nc
    sum=sum+(ym(:,ic)-yc(:,ic))'*(ym(:,ic)-yc(:,ic));
end
if sum< SumSave*0.9
    disp([' Function evaluation No. ' num2str(Feval) '  Sum of Squres ' num2str(sum)]);
    SumSave=sum;
end
%- - - - - - - - - - - - - - - - - - - - - -

function dYfuncvecdt = ODEfun(t,Yfuncvec,Beta); 
y1 = Yfuncvec(1); 
y2 = Yfuncvec(2); 
b1=Beta(1);
b2=Beta(2);
b3=Beta(3);
b4=Beta(4);
dy1dt = b1 * y1 * (1 - y1 / b2); 
dy2dt = b3 * y1 - b4 * y2; 
dYfuncvecdt = [dy1dt; dy2dt];