function DYDT=f(t,y)
% This is file cstrTeq.m
% INPUT: Time t and state vector y
% OUTPUT: derivatives DYDT
% Usage with ode15s:   options.Mass=diag([1 1 1 0],0);
%    [T,Y]= ode15s(@cstrTeq,[0 50],[1454 8546 453.4 854.6],options)
% Plot: 	plot(T,Y(:,3))

% All in SI units except time which is in minutes.
% I. Data (parameters and independent variables)
cAf=10000;             % feed concentration of A [mol/m3]
cBf=0;                 % feed concentration of B [mol/m3]
Tf=300;                % feed temperature [K]
V = 0.9;               % reactor volume [m3]
q = 0.12;              % volumetric flowrate [m3/min] (KAN ENDRES)
dhr=  -80e3;           % Heat of reaction [J/mol]
cp= 4.e3;              % heat capacity [J/kg]
rho = 1000;            % mean density [kg/m3]
U = 1000*60;           % heat transfer coefficient [J/min,m2,K]
A = 5;                 % heat transfer area [m2]
Tc = 430;              % cooling temperature [K]
k1_400 = 0.1;          % rx1: rate constant at 400K [1/min]
k2_400 = 0.001;        % rx2: rate constant at 400K [1/min]
E1 = 60e3;             % rx1: activation energy [J/mol]
E2 = E1 - dhr;         % rx2: activation energy
R = 8.31;              % [J/K mol]

% II. Extract present value of states
cA=y(1);
cB=y(2);
T=y(3);
xi=y(4);

% IIIa. Intermediate calculations
k1 = k1_400 * exp(-(E1/R) * (1/T - 1/400));
k2 = k2_400 * exp(-(E2/R) * (1/T - 1/400));
K = k1/k2  % NEW! Eq. constant
Q = U*A*(Tc-T);

% IIIb. Evaluate derivatives of states
Vdcadt = q*cAf - q*cA - xi;                    % [mol A/min]
Vdcbdt = q*cBf - q*cB + xi;                    % [mol B/min]
mcpdTdt = rho*q*cp*(Tf-T) + xi*(-dhr) + Q;     % [J/min]

f1 = Vdcadt/V;
f2 = Vdcbdt/V;
f3 = mcpdTdt / (rho*cp*V);
f4 = K - cB/cA;  % NEW 

DYDT=[f1; f2; f3; f4];  

