load o426.dat
t=o426(:,1);
y=o426(:,2);

plot(t,y,'k')
ylabel('x_1')
xlabel('t')
set(1, 'PaperUnits', 'centimeters')
set(1, 'PaperPosition', [3 10 11 8.25])
print -deps 'id21exm1a.eps'


du = 10;
yn = y/du;
i1=find(t>0.8);
i=i1(1);
t1=t(i); y1=yn(i);
t1 = fix(t1*100+0.5)/100;
y1 = fix(y1*10000+0.5)/10000;
i1=find(t>2);
j=i1(1);
t2=t(j); y2=yn(j);
t2 = fix(t2*100+0.5)/100;
y2 = fix(y2*10000+0.5)/10000;
K = yn(end);
T= (t2-t1)/(log((K-y1)/(K-y2)));
x=(log((K-y1)/K))/(log((K-y2)/K));
D= (t2*x - t1)/(x-1);
T = fix(T*100+0.5)/100;
D = fix(D*100+0.5)/100;

g1=tf(K, [T 1]);
[yo,to]=step(g1,t(end));
to = to + D;
i1=find(to<=t(end));
yod=yo(i1);
tod=to(i1);
plot(t,yn,'k',tod,yod,'k--')
text(t1,y1,'+', 'HorizontalAlignment', 'center')
text(t2,y2,'+', 'HorizontalAlignment', 'center')
ylabel('x_1')
xlabel('t')
legend('x_1', 'x_1^n', 4)
set(1, 'PaperUnits', 'centimeters')
set(1, 'PaperPosition', [3 10 11 8.25])
print -deps 'id21exm1b.eps'
