% running distillation model 'multivessel' for a ternary mixture
load Xinit                                    % loads state vector from the start_up period
% Number of stages in the top section: 
NT1=25;
% Number of stages in the bottom section (without the reboiler): 
NT2=25;
% Total number of 'stages'  
NT=NT1+NT2+3;                                 % We include decanter, middle vessel and reboiler 

% Initial conditions
zFA0=Xinit(1:NT);                             % Initial composition of light component
zFB0=Xinit(NT+1:2*NT);                        % Initial composition of intermediate component 
M0=Xinit(2*NT+1:3*NT);                        % Initial molar holdup in each stage (kmol) 
T0=Xinit(3*NT+1:4*NT);                        % Initial temperature in each stage (K)
IE10=0;

% Defining options and mass matrix for ODE solver 
M=eye(4*NT+1,4*NT+1);
M(3*NT+1:4*NT,3*NT+1:4*NT)=0;
options.Mass=M; 

% Running column
[t,x]=ode15s(@input_multivessel,[0 5.5],[zFA0;zFB0;M0;T0;IE10],options);

% Save results
save x; 
lengthx=size(x); Xss=x(lengthx(1),:)';
save Xss;
Results=[Xss(1+6+2*NT1+2*NT2), Xss(2+6+3*NT1+2*NT2), Xss(9+3*NT1+3*NT2); Xss(1), Xss(2+NT1), Xss(3+NT1+NT2);...
         Xss(1+3+NT1+NT2), Xss(2+3+2*NT1+NT2), Xss(6+2*NT1+2*NT2);...
         1-Xss(1)-Xss(3+NT1+NT2+1), 1-Xss(2+NT1)-Xss(3+NT1+NT2+2+NT1), 1-Xss(3+NT1+NT2)-Xss(6+2*NT1+2*NT2)];