

function DXDT = f(t,x)

% This is the file dist.m  (Version 23Feb2000)
% Distillation column with reboiler (stage 1), a feed stage  (stage 2),
%  ... a stage above this ( stage 3) and a total condenser (stage 4)
% Assumptions; Binary mixture with constant alfa and constant molar flows
%              Molar holdup on all stages is 1 (M=1)
% States x : vector of liquid mole fractions of light component on the stages
%
% Usage:
%  x0 = [0.5, 0.5, 0.5, 0.5]; % initial states (not steady-state)
%  [T,X] = ode15s('dist',[0 1000],x0) % First simulate to t=1000 (steady-state)
%  xss = X(length(X),:) % Save the steady-state mole fractions
%  [T,X] = ode15s('dist',[0 20],xss) % Run new simulation (e.g change F=1.2)

% Assume constant relative volatility
alfa = 4.78;

% Feed rate and feed composition  (may change this)
F=1.0; zF=0.5;

% Flows in the column (feed liquid; constant molar flows)
V=3.55; V1=V; V2=V; V3=V;
L=3.05; L4=L; L3=L; L2=L+F;
D=V3-L4; B=L2-V1;

% Vapor-liquid equlibrium (constant relative voloatility)
y(1) = alfa*x(1)/(1+(alfa-1)*x(1));
y(2) = alfa*x(2)/(1+(alfa-1)*x(2));
y(3) = alfa*x(3)/(1+(alfa-1)*x(3));
y(4)=x(4);  % total condenser

% Component balances
DXDT(1) = L2*x(2)-V1*y(1)-B*x(1);
DXDT(2) = L3*x(3)+V1*y(1)-L2*x(2)-V2*y(2)+F*zF;
DXDT(3) = L4*x(4)+V2*y(2)-L3*x(3)-V3*y(3);
DXDT(4) = V3*y(3)-L4*x(4)-D*x(4);

% Change vector DXDT to a column vector (MATLAB requires this...)
DXDT=DXDT';

