clear; clc;
to=0; tf=200; yo=[5;50;50;0];
% y = [F1; T; Ts; Error Integral]
% Note: estimate Kp=DeltaT/DeltaF1=(75-50)/(0-5) = -5.0 o s/kg; taup = M/(F1+F2)=100/(5+5)=10 sec 
soln=ode15s(@f, [to tf], yo);
x=linspace(0,tf,100);
y1=deval(soln,x,1);
y2=deval(soln,x,2);
y3=deval(soln,x,3);
plot(x,y2,'k-','linewidth',2);
%hold on;
%plot(x,y3,'r-','linewidth',2);
xlabel('t (s)');ylabel('T (C)');

function dydt=f(t,y)
F2=5;M=100;T1=25;T2=75;tauv=2;tauts=6;
% PID settings (Kc, tauI, tauD) can be adjusted in the line below, then
% run the program again to see the closed-loop response with new settings
Kc=-0.82;tauI=14.0;tauD=2.1;Tsp=55;
if t<10; Tsp=50; end
% derivative kick is suppressed in the line below which sets de/dt =
% -dTs/dt
erder=-(1/tauts)*(y(2)-y(3));
er=Tsp-y(3);
F1spec=5+Kc*(er+y(4)/tauI+tauD*erder);
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)=Tsp-y(3);
dydt=dydt';
end