% programa 5.6, W.L.Luyben, pag. 126
% "nonisothermal cstr" 
% modif. E.Quintero.M 18nov'98
%--------------------------------------------------
% disturbance is step change in feed composition at
% time zero from 0.50 to 0.55
clear
format short e
% initial conditions
ca=0.245; t=600; tj=594.59; xx=[];
v=48; time=0; vc=v*ca; vt=v*t; 
% parameter values
tjo=530; fo=40; to=530;
kc=4; delta=0.01; tprint=0;
% disturbance 0.55
cao=0.55; %s.s. 0.50
while(time<6)
%  feedback controllers   
   fj=49.9-kc*(600-t);
   f=40-10*(48-v);
%  reaction rate
   k=7.08e10*exp(-30000/(1.99*t));
   q=150*250*(t-tj);
%  evaluate all derivatives
   vdot=fo-f;
   vcdot=fo*cao-f*ca-v*k*ca;
   vtdot=fo*to-f*t+(30000*v*k*ca-q)/(0.75*50);
   tjdot=fj*(tjo-tj)/3.85+q/240;
   if(time>=tprint)
      xx=[xx;time,ca,t,v,f,tj,fj,q];
      tprint=tprint+0.2;
   end
   v=v+vdot*delta;
   vc=vc+vcdot*delta;
   vt=vt+vtdot*delta;
   tj=tj+tjdot*delta;
   time=time+delta;
   ca=vc/v;
   t=vt/v;
end
xx
subplot(4,2,1),plot(xx(:,1),xx(:,7)),ylabel('Fj'),grid
subplot(4,2,3),plot(xx(:,1),xx(:,8)),ylabel('Qc'),grid
subplot(4,2,5),plot(xx(:,1),xx(:,6)),ylabel('Tj'),grid
subplot(4,2,7),plot(xx(:,1),xx(:,5)),ylabel('F'),grid
xlabel('Time, h')
subplot(4,2,2),plot(xx(:,1),xx(:,4)),ylabel('V'),grid
subplot(4,2,4),plot(xx(:,1),xx(:,3)),ylabel('T'),grid
subplot(4,2,6),plot(xx(:,1),xx(:,2)),ylabel('Ca'),grid
xlabel('Time, h')