function fig8_16

% fit logistic equation using RK4 + cubic spline + nead-medler

global td yd tmax


%  generate data for regression
a=1; alpha0=2; beta0=3;
n=80;
tmax=2;
%td=linspace(0,tmax,n);
tdd=tmax*rand(1,n);
td=sort(tdd);
qd=0.2*(rand(1,n)-0.5);
yd=zeros(1,n);
td(1)=0; yd(1)=a;
for in=2:n
    yd(in)=(1+qd(in))*a*beta0/(a+(beta0-a)*exp(-alpha0*beta0*td(in)));
end

[sol,fv,ee,o] = fminsearch(@F, [1, 2],optimset('MaxFunEvals',1000000,'MaxIter',1000000,'TolX',1e-6))
alpha=sol(1); beta=sol(2);
fprintf('\n  Computed Solution =  %e  %e\n\n',alpha, beta)

nt=100;
y0=a;
tt=linspace(0,tmax,nt);
k=tt(2)-tt(1);
yy=rk4(tt,y0,k,nt,alpha,beta);

clf
% get(gcf)
set(gcf,'Position', [1 1115 658 230])
subaxis(1,1,1,1,'MT',0.01,'MB',0.16,'MR',-0.02,'ML',0.03,'P',0.04)
co = [0 0 1;
    0 0.5 0;
    1 0 0;
    0 0.75 0.75;
    0.75 0 0.75;
    0.75 0.75 0;
    0.25 0.25 0.25];
set(groot,'defaultAxesColorOrder',co)
plot(tt,yy,'r','LineWidth',1.6)
hold on
grid on
axis([0 tmax 1 3.5])
plot(td,yd,'b.','MarkerSize',20)
xlabel('t-axis')
ylabel('y-axis')

set(gca,'FontSize',16,'FontWeight','bold')

%exportgraphics(gca,'/Users/mark/Desktop/fitle.eps')

%function s=F(alpha,beta,td,yd)
function s=F(p)
global td yd tmax
alpha=p(1); beta=p(2);
y0=1;
m=6;
t=linspace(0,tmax,m+1);
k=t(2)-t(1);
y_rk4=rk4(t,y0,k,m+1,alpha,beta);
dy0=f(0,y0,alpha,beta);
dy1=f(tmax,y_rk4(m+1),alpha,beta);
yf = splineA(t,y_rk4,td,[1 1],[dy0 dy1]);
s=norm(yf-yd,2)^2;

function dy=f(t,y,alpha,beta)
dy(1)=alpha*y*(beta-y);

% RK4 method
function ypoints=rk4(t,y0,h,n,alpha,beta)
y=y0;
ypoints=y0;
for i=2:n
    k1=h*f(t(i-1),y,alpha,beta);
    k2=h*f(t(i-1)+0.5*h,y+0.5*k1,alpha,beta);
    k3=h*f(t(i-1)+0.5*h,y+0.5*k2,alpha,beta);
    k4=h*f(t(i),y+k3,alpha,beta);
    yy=y+(k1+2*k2+2*k3+k4)/6;
    ypoints=[ypoints, yy];
    y=yy;
end


















