function ydot = reactor(t,y)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% This function defines the catalytic reactor dynamics.
%
% Input variables: t - time
% y - state column vector
% y(1) = P
% y(2) = T
% y(3) = Pp
% y(4) = Tp
%
% Output variables: ydot - state derivatives column vector
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Reaction rate constant
K = 6e-4*exp(20.7-15000/y(2));
ydot(1) = 0.1+320*y(3)-321*y(1);
ydot(2) = 1752+266.667*y(4)-269.267*y(2);
ydot(3) = 1866.8*(y(1)-(1+K)*y(3));
ydot(4) = 10369*K*y(3)+1.2964*(y(2)-y(4));
ydot = ydot(:);
end;