function xprime=colamod(t,X,U) 

% Parameters

% Number of stages (including reboiler and total condenser: 
NT = 41; 

% Location of feed stage (stages are counted from the bottom):
NF = 21;

% Relative volatilities: alpha = [alpha1 alpha2...alphaNC]
alpha = [2.0 1.5 1.33];

% Number of components
NC = length(alpha);
zF = [0.4 0.2 0.4];
% Parameters for Franci's Weir Formula L(i) = K*Mow(i)^1.5
Kuf = 21.65032;         % Constant above feed
Kbf = 29.65032;         % Constant below feed
Muw = 0.25;             % Liquid holdup under weir (kmol)

% P-controllers for control of reboiler and condenser hold up.
KcB = 10;  KcD = 10;                     % controller gains
MDs = 0.5; MBs = 0.5;                    % Nominal holdups - these are rather small  
Ds  = 0.5; Bs  = 0.5;                    % Nominal flows

% Splitting the states
x1 = X(1:82);% Liquid compositions from btm to top
x2 = X(124:205);
M1 = X(83:123);              % Liquid hold up from btm to top
M2 = X(206:246);              % Liquid hold up from btm to top

% Rearrange elements of composition vector (x) for later use
Iu = reshape(1:(NC-1)*NT,NT,NC-1);
x1  = x1(Iu); 
x2  = x2(Iu); 
LT1 = U(1) ;
VB1 = U(2);
LT2 = U(3); 
VB2 = U(4); 
D1 = U(5); 
B1 = U(6);
D2 = U(7); 
B2 = U(8);
F = U(9);
qF = U(10);


% The distillation model

% Vapour-liquid equilibria (multicomponent ideal VLE, Stichlmair-Fair, 'Distillation', p. 36, 1998)
y1 = (x1*diag(alpha(1:NC-1)))./((x1*(alpha(1:NC-1) - 1)' + 1)*ones(1,NC-1));

% Vapor flows assuming constant molar flows
V1 = VB1*ones(NT-1,1); V1(NF:NT-1) = V1(NF:NT-1) + (1-qF)*F;

% Liquid flows are given by Franci's Weir Formula L(i)=K*Mow(i)^1.5 
% Liquid flow L(i) dependent only on the holdup over the weir Mow(i)
% M(i)= Mow(i) + Muw(i) (Total holdup = holdup over weir + holdup below weir) 
L1(2:NF,1)      = Kbf*(M1(2:NF) - Muw).^1.5;         % Liquid flows below feed (Kmol/min)
L1(NF+1:NT-1,1) = Kuf*(M1(NF+1:NT-1) - Muw).^1.5;    % Liquid flows above feed (Kmol/min)
L1(NT,1)        = LT1;                               % Condenser's liquid flow (Kmol/min)

% Time derivatives from material balances for: (1) total holdup and (2) component holdup
% Column
j = 2:NT-1; dMdt1(j,1)  = L1(j+1) - L1(j) + V1(j-1) - V1(j);
i = 1:NC-1; dMxdt1(j,i) = diag(L1(j+1))*x1(j+1,i) - diag(L1(j))*x1(j,i) + diag(V1(j-1))*y1(j-1,i) - diag(V1(j))*y1(j,i);

% Correction for feed at the feed stage: The feed is assumed to be mixed into the feed stage
dMdt1(NF)    = dMdt1(NF) + F;
dMxdt1(NF,:) = dMxdt1(NF,:) + F*zF(1,1:NC-1);

% Reboiler (assumed to be an equilibrium stage)
dMdt1(1)     = L1(2) - V1(1) - B1;
dMxdt1(1,i)  = L1(2)*x1(2,i) - V1(1)*y1(1,i) - B1*x1(1,i);

% Total condenser (no equilibrium stage)
dMdt1(NT)    = V1(NT-1) - LT1 - D1;
dMxdt1(NT,i) = V1(NT-1)*y1(NT-1,i) - L1(NT)*x1(NT,i) - D1*x1(NT,i);

% Compute the derivative for the mole fractions from d(Mx) = x dM + M dx
dxdt1 = inv(diag(M1))*(dMxdt1 - diag(dMdt1)*x1);

% Rearrange elements of composition vector derivatives (dxdt)
dxdt1 = dxdt1(:);



% The distillation model

% Vapour-liquid equilibria (multicomponent ideal VLE, Stichlmair-Fair, 'Distillation', p. 36, 1998)
y2 = (x2*diag(alpha(1:NC-1)))./((x2*(alpha(1:NC-1) - 1)' + 1)*ones(1,NC-1));

% Vapor flows assuming constant molar flows
V2 = VB2*ones(NT-1,1); V2(NF:NT-1) = V2(NF:NT-1); %+ (1-qF)*B1;

% Liquid flows are given by Franci's Weir Formula L(i)=K*Mow(i)^1.5 
% Liquid flow L(i) dependent only on the holdup over the weir Mow(i)
% M(i)= Mow(i) + Muw(i) (Total holdup = holdup over weir + holdup below weir) 
L2(2:NF,1)      = Kbf*(M2(2:NF) - Muw).^1.5;         % Liquid flows below feed (Kmol/min)
L2(NF+1:NT-1,1) = Kuf*(M2(NF+1:NT-1) - Muw).^1.5;    % Liquid flows above feed (Kmol/min)
L2(NT,1)        = LT2;                               % Condenser's liquid flow (Kmol/min)

% Time derivatives from material balances for: (1) total holdup and (2) component holdup
% Column
j = 2:NT-1; dMdt2(j,1)  = L2(j+1) - L2(j) + V2(j-1) - V2(j);
i = 1:NC-1; dMxdt2(j,i) = diag(L2(j+1))*x2(j+1,i) - diag(L2(j))*x2(j,i) + diag(V2(j-1))*y2(j-1,i) - diag(V2(j))*y2(j,i);

% Correction for feed at the feed stage: The feed is assumed to be mixed into the feed stage
dMdt2(NF)    = dMdt2(NF) + B1;
s_zF = [x1(1,1) x1(1,2)];
dMxdt2(NF,:) = dMxdt2(NF,:) + B1*s_zF(1,1:NC-1);

% Reboiler (assumed to be an equilibrium stage)
dMdt2(1)     = L2(2) - V2(1) - B2;
dMxdt2(1,i)  = L2(2)*x2(2,i) - V2(1)*y2(1,i) - B2*x2(1,i);

% Total condenser (no equilibrium stage)
dMdt2(NT)    = V2(NT-1) - LT2 - D2;
dMxdt2(NT,i) = V2(NT-1)*y2(NT-1,i) - L2(NT)*x2(NT,i) - D2*x2(NT,i);

% Compute the derivative for the mole fractions from d(Mx) = x dM + M dx
dxdt2 = inv(diag(M2))*(dMxdt2 - diag(dMdt2)*x2);

% Rearrange elements of composition vector derivatives (dxdt)
dxdt2 = dxdt2(:);

xprime=[dxdt1; dMdt1; dxdt2; dMdt2];
