%filename tempdyn.m

function Tdot=tempdyn(t,T)
global qsetpt taud taui Kc Tsetpt onoff

%  Use logical block to model the step change at 10 min. 
if t<10
    Tinlet = 60;
else
   Tinlet = 40;
end

qin= qsetpt+Kc*(Tsetpt-T(3))+onoff*Kc/taui*T(4);% total heat sent in 
% use the following statement for part (e)
%qin=max(0,min(2.6*qsetpt,qin));

row(1)= (500*(Tinlet-T(1))+qin)/(4000);	% energy balance
row(2) = (T(1)-T(2)- 0.5*taud*row(1))*2/taud;	% Pade approximation for delay
row(3) = (T(2)-T(3))/5;			% Thermocouple dynamics
row(4) = Tsetpt - T(3);			% the integrated error
% row(4) not needed for part (e), but is calculated anyway
Tdot = row';