clear; clc;
to=0; tf=150; yo=[5;50;50;5];
% vector y = [F1; T; Ts; F1spec,est]
soln=ode15s(@f, [to tf], yo);
fig1=figure(1);
x=linspace(0,tf,150);
y1=deval(soln,x,1);
y2=deval(soln,x,2);
y3=deval(soln,x,3);
y4=deval(soln,x,4);
subplot(2,1,1);
plot(x,y3,'k-',x,y2,'r-','linewidth',1);
legend('Ts','T','Location','northeast');
xlabel('time (s)');ylabel('o C');
subplot(2,1,2);
plot(x,y1,'r-',x,y4,'k-','linewidth',1);
legend('F1','F1spec','Location','northeast');
xlabel('time (s)');ylabel('kg/s');

function dydt=f(t,y)
F2=5.0;M=100.0;T1=25.0;T2=75.0;tauv=2.0;tauts=6.0;F10=5.0;
% line below sets the F1 perturbation "h" 
h=0.1;T0=50;
F1spec=F10;
% next four lines implement the F1spec variation: for time 0 to 10 F1spec is
% kept at the initial value F10, then from 10 to 20 F1spec is decreased by 
% h to cause T (i.e. y(2)) to deviate from its steady state value T0 (this starts the ATV process), 
% and after that F1spec is cycled by h above/below F10 each time the measured signal Ts (i.e. y(3)) crosses T0
if T0-y(3)>=0; F1spec=F10-h; end
if T0-y(3)<0; F1spec=F10+h; end
if t<20; F1spec=F10-h; end
if t<10; F1spec=F10; end
dydt(1)=(1/tauv)*(F1spec-y(1));
dydt(2)=(1/M)*(y(1)*T1+F2*T2-(y(1)+F2)*y(2));
dydt(3)=(1/tauts)*(y(2)-y(3));
dydt(4)=(1/.01)*(F1spec-y(4));
dydt=dydt';
end