function r=VapLiqMix(Var,No,Siv,Sil,Sov,Sol,Mode)

GlobalVars;

Neq=2*NC+6;                        % Number of equations
S=sparse(diag([
  1.1030e+000
  2.0252e+000
  4.1250e+000
  7.2497e+000
  4.1250e+000
  6.9591e-005
  1.4406e-002
  8.4323e-002
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e+000
  1.0000e-001
  8.3333e-002
  1.2500e-003
]));
if Unscaled
  S=speye(size(S));
end

I=ones(1,NC);

% Variables
% Vapor feed stream
Ix=IxStr(Siv);
Ixxiv = Ix(1:NC);   xiv = Var(Ixxiv);
IxFiv = Ix(NC+1);   Fiv = Var(IxFiv);
IxTiv = Ix(NC+2);   Tiv = Var(IxTiv);
IxPiv = Ix(NC+3);   Piv = Var(IxPiv);

% Liquid feed stream
Ix=IxStr(Sil);
Ixxil = Ix(1:NC);   xil = Var(Ixxil);
IxFil = Ix(NC+1);   Fil = Var(IxFil);
IxTil = Ix(NC+2);   Til = Var(IxTil);
IxPil = Ix(NC+3);   Pil = Var(IxPil);

% Vapor feed stream
Ix=IxStr(Sov);
Ixxov = Ix(1:NC);   xov = Var(Ixxov);
IxFov = Ix(NC+1);   Fov = Var(IxFov);
IxTov = Ix(NC+2);   Tov = Var(IxTov);
IxPov = Ix(NC+3);   Pov = Var(IxPov);

% Liquid feed stream
Ix=IxStr(Sol);
Ixxol = Ix(1:NC);   xol = Var(Ixxol);
IxFol = Ix(NC+1);   Fol = Var(IxFol);
IxTol = Ix(NC+2);   Tol = Var(IxTol);
IxPol = Ix(NC+3);   Pol = Var(IxPol);

if strcmp(Mode,'Mass')
    %Compute or define the mass matrix for this model
    M=spalloc(Neq,Nx,0);
    r = M;
    
elseif strcmp(Mode,'InitVal')
    % Return initial values for variables reciding in this module
    Var(Ixxov) = [0.94 .01 .02 .02 .01]';            % Vapour stream composition   
    Var(IxFov) = 0.2;                                % Vapour stream component flow     [kmol/s]
    Var(IxTov) = 359.15;                             % Vapour stream temperature        [C] 
    Var(IxPov) = 10;                                 % Vapour stream pressure           [bar a]

    Var(Ixxol) = [0.02 .00 .30 .54 .14]';            % Liquid stream composition
    Var(IxFol) = 0.2;                                % Liquid stream component flow     [kmol/s]
    Var(IxTol) = 359.15;                             % Liquid stream temperature        [C] 
    Var(IxPol) = 10;                                 % Liquid stream pressure           [bar a]
    
    r = Var;
    
elseif strcmp(Mode,'Jacobian')
    % Compute Jacobian matrix
    J=spalloc(Neq,NVar,2*NC^2+8*NC+9);
    II=speye(NC);
    hiv=h(Tiv,1);
    hil=h(Til,0);
    hov=h(Tov,1);
    hol=h(Tol,0);
    dhiv=dh_dT(Tiv,1);
    dhil=dh_dT(Til,0);
    dhov=dh_dT(Tov,1);
    dhol=dh_dT(Tol,0);

    % Component balance
    J(1:NC,Ixxiv)= Fiv*II;        J(1:NC,Ixxil)= Fil*II;        J(1:NC,Ixxov)=-Fov*II;        J(1:NC,Ixxol)=-Fol*II;
    J(1:NC,IxFiv)= xiv;           J(1:NC,IxFil)= xil;           J(1:NC,IxFov)=-xov;           J(1:NC,IxFol)=-xol;
    % Energy balance
    J(NC+1,Ixxiv)= Fiv*hiv';      J(NC+1,Ixxil)= Fil*hil';      J(NC+1,Ixxov)=-Fov*hov';      J(NC+1,Ixxol)=-Fol*hol';
    J(NC+1,IxFiv)= xiv'*hiv;      J(NC+1,IxFil)= xil'*hil;      J(NC+1,IxFov)=-xov'*hov;      J(NC+1,IxFol)=-xol'*hol;
    J(NC+1,IxTiv)= Fiv*xiv'*dhiv; J(NC+1,IxTil)= Fil*xil'*dhil; J(NC+1,IxTov)=-Fov*xov'*dhov; J(NC+1,IxTol)=-Fol*xol'*dhol;
    % Equlibrium 
    J(NC+2:2*NC+1,Ixxov)= II;                                   J(NC+2:2*NC+1,Ixxol) =-K(Tov,Pov); 
    J(NC+2:2*NC+1,IxTov)=-dK_dT(Tov,Pov)*xol;                   J(NC+2:2*NC+1,IxPov) =-dK_dP(Tov,Pov)*xol;
    % Sum of components equal to one
    J(2*NC+2,Ixxov)= I;                  
    J(2*NC+3,Ixxol)= I;
    % Equal stream pressure
    J(2*NC+4,IxPil)= 1;  J(2*NC+4,IxPov)=-1;
    J(2*NC+5,IxPiv)= 1;  J(2*NC+5,IxPov)=-1;
    % Equal outlet temperature
    J(2*NC+6,IxTov)= 1;  J(2*NC+6,IxTol)=-1;

    r = S*J;

else
    % Compute dxdt
    r=spalloc(Neq,1,Neq);
    hiv=h(Tiv,1);
    hil=h(Til,0);
    hov=h(Tov,1);
    hol=h(Tol,0);
    
    r =[Fiv*xiv+Fil*xil-Fov*xov-Fol*xol                      % Component balance          [kmol/h]
        Fiv*xiv'*hiv+Fil*xil'*hil-Fov*xov'*hov-Fol*xol'*hol  % Energy balance             [kJ/h]
        xov-K(Tov,Pov)*xol                                   % Equlibrium 
        I*xov-1                                              % Sum of components equal to one
        I*xol-1                                              % Sum of components equal to one
        Pil-Pov                                              % Equal pressure
        Piv-Pov                                              % Equal pressure
        Tov-Tol];                                            % Outlet stream temperature  [C]
   
    r = S*r;
end

