function cdot = ex56(t, c)
% 
% this file defines the equation set in Figure 5.4,

global coeff T xb xt xz JT right_hand

 xb(1)=c(1); xt(1)=c(2);
 xb(2)=c(3); xt(2)=c(4);
 xb(3)=c(5); xt(3)=c(6);
 xb(4)=c(7); xt(4)=c(8);
 xb(5)=c(9); xt(5)=c(10);
 for I = 1:5,
   xz(I) = 1-xb(I)-xt(I);
 end

holdup(1) = 0.05;
holdup(2) = 0.0125;
holdup(3) = 0.0125;
holdup(4) = 0.0125;
holdup(5) = 0.6;
 MWb = 78.1; MWt = 92.1; MWz = 106.2;
 lambda = 155;
 Q1 = 85000;
 F = 2.5;
 reflux = 2.0;
 xtf = .3;
 TF = 250;
 P = 1;
 xbf = (P -feval('Pz',TF)-xtf*(feval('Pt',TF)-feval('Pz',TF)))...
       /(feval('Pb',TF)-feval('Pz',TF));
 QF = F*0.5*TF*(MWb*xbf +MWt*xtf +MWz*(1-xbf-xtf));
%
% solve nonlinear equations to get y an T
%
for J = 1:5,
   JT = J;
   T(J) = fzero('impT',T(J));
   yb(J) = feval('Pb',T(J))/P*xb(J);
   yt(J) = feval('Pt',T(J))/P*xt(J);
   yz(J) = feval('Pz',T(J))/P*xz(J);
end
%
% calc HL and HV
%
for J = 1:5,
   HL(J) = (MWb*xb(J)+MWt*xt(J)+MWz*xz(J))*(.5*T(J));
   HV(J) = (MWb*yb(J)+MWt*yt(J)+MWz*yz(J))*(.5*T(J)+lambda);
end
%
% update matrix
%
 coeff(3,1) = HL(1);
 coeff(3,3) = -HL(2);
 coeff(3,4) = -HV(2);
 coeff(3,6) = HV(3);
 coeff(5,3) = HL(2);
 coeff(5,5) = -HL(3);
 coeff(5,6) = -HV(3);
 coeff(5,8) = HV(4);
 coeff(7,5) = HL(3);
 coeff(7,7) = -HL(4);
 coeff(7,8) = -HV(4);
 coeff(7,10) = HV(5);
 coeff(9,7) = HL(4);
 coeff(9,9) = -HL(5);
 coeff(9,10) = -HV(5);

LV = coeff\right_hand;
L(1)=LV(1); L(2)=LV(3); L(3)=LV(5); L(4)=LV(7); L(5)=LV(9);
V(1)=LV(2); V(2)=LV(4); V(3)=LV(6); V(4)=LV(8); V(5)=LV(10);
%
% derivatives
%
cdot(1) = (V(2)*yb(2) - L(1)*c(1) - V(1)*yb(1))/holdup(1);
cdot(2) = (V(2)*yt(2) - L(1)*c(2) - V(1)*yt(1))/holdup(1);
cdot(3) = (V(3)*yb(3) + L(1)*c(1) - V(2)*yb(2) - L(2)*c(3))...
          /holdup(2);
cdot(4) = (V(3)*yt(3) + L(1)*c(2) - V(2)*yt(2) - L(2)*c(4))...
          /holdup(2);
cdot(5) = (L(2)*c(3) + V(4)*yb(4) + F*xbf - V(3)*yb(3)...
           -L(3)*c(5))/holdup(3);
cdot(6) = (L(2)*c(4) + V(4)*yt(4) + F*xtf - V(3)*yt(3)...
          -L(3)*c(6))/holdup(3);
cdot(7) = (L(3)*c(5) + V(5)*yb(5) - L(4)*c(7) - V(4)*yb(4))...
          /holdup(4);
cdot(8) = (L(3)*c(6) + V(5)*yt(5) - L(4)*c(8) - V(4)*yt(4))...
          /holdup(4);
cdot(9) = (L(4)*c(7) - V(5)*yb(5) - L(5)*c(9))/holdup(5);
cdot(10) = (L(4)*c(8) - V(5)*yt(5) - L(5)*c(10))/holdup(5);
cdot = cdot(:);
