diff --git a/Examples/testcase3.m b/Examples/testcase3.m
new file mode 100644
index 0000000..fcacad4
--- /dev/null
+++ b/Examples/testcase3.m
@@ -0,0 +1,26 @@
+% test case
+% liquid-two-phase phase transition
+clc; clear;
+% Define the components and load pure physsical properties
+[component, comp_flag] = addComponents({'dimethyl ether', 'water', 'C10H22'});
+% Define the thermodynamic models
+T0 = 300; % [K]
+p0 = 100e5; % [Pa]
+thermo1 = addThermo();
+thermo1.EOS = @PREOS;
+mixture1 = addMixture(component, T0, p0);
+mixture1.mole_fraction = [0.1, 0.45, 0.45];
+BIP = struct();
+BIP.EOScons = [0 0.01 0.1; 0.01 0 0.1; 0.1 0.1 0];
+BIP.EOStdep = zeros(3);
+mixture1.bip = BIP;
+% Define flash options
+options.accuracy = 1e-7;
+options.iteration = 100;
+% Negative flash
+[vapor_y, liquid_x, vapor_frac] = ...
+        lleflash(mixture1, thermo1, options)
+%mixture1.temperature=T0;
+[liquid_z, vapor_z, fugacity, HR] = PREOS(mixture1, thermo1);    
+
+[s1, sl, sv] = stabilityLLETest(mixture1, thermo1)
diff --git a/Flash/kvalueLLE.m b/Flash/kvalueLLE.m
index 9c0f546..e22be44 100644
--- a/Flash/kvalueLLE.m
+++ b/Flash/kvalueLLE.m
@@ -1,5 +1,4 @@
-function ki = kvalueLLE(component, BIP, liquid_x, vapor_y, pressure, ...
-        temperature, eosf, mixrule, activityfun)
+function ki = kvalueLLE(mixture, thermo, liquid_x, vapor_y)
 % 
 % SYNOPSIS:
 %   
@@ -44,10 +43,15 @@
 %}
 
 
-[zl,zv,liq_fug,hh]=eosf(component, BIP, liquid_x, pressure, temperature, ...
-        mixrule, activityfun, 1, 1);
-[zl,zv,vap_fug,hh]=eosf(component, BIP, vapor_y, pressure, temperature, ...
-        mixrule, activityfun, 1, 1);
+thermo.phase=1;
+thermo.fugacity_switch=1; % switch on
+eosf = thermo.EOS;
+mixture1 = mixture;
+mixture1.mole_fraction = liquid_x;
+[zl,zv,liq_fug,hh]=eosf(mixture1, thermo);
+thermo.phase=1;
+mixture1.mole_fraction = vapor_y;
+[zl,zv,vap_fug,hh]=eosf(mixture1, thermo);
 % [zl,zv,vap_fug,hh]=EOS(critical_temp,critical_pres,vapor_y, ...
 %     acentric_fact,[0 0;0 0],pressure,temperature,eos_type,1,2);
 N = length(liquid_x);
diff --git a/Flash/lleflash.m b/Flash/lleflash.m
new file mode 100644
index 0000000..16d272f
--- /dev/null
+++ b/Flash/lleflash.m
@@ -0,0 +1,103 @@
+function [vapor_y, liquid_x, vapor_frac]=vleflash(mixture, thermo, options)
+%IMPORTANT: every variable should be in the form of row vectors [1 x N]
+% 
+% SYNOPSIS:
+%   
+% 
+% PARAMETERS:
+%   
+% 
+% RETURNS:
+%   
+% 
+% EXAMPLE:
+% 
+% SEE ALSO:
+%     
+
+%{
+Copyright (c) 2012, 2013, Ali Akbar Eftekhari
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or 
+without modification, are permitted provided that the following 
+conditions are met:
+
+    *   Redistributions of source code must retain the above copyright notice, 
+        this list of conditions and the following disclaimer.
+    *   Redistributions in binary form must reproduce the above 
+        copyright notice, this list of conditions and the following 
+        disclaimer in the documentation and/or other materials provided 
+        with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 
+PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR 
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, 
+OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%}
+
+    
+    eps1 = options.accuracy; 
+    max_itr = options.iteration;
+    %Initial estimate for k-values using an empirical equation
+    ki = kval_estimate(mixture);
+    composition = mixture.mole_fraction;
+    %Initialization of variables:
+    vapor_frac=0.5;
+    error1=1;
+    error2=1;
+    error3=1;
+    j=0;
+    while ((error1>eps1) || (error2>eps1) || (error3>eps1))
+    j=j+1;
+    if (j>max_itr)  %check the maximum number of itr to avoid infinite loop
+      break;
+    end
+    [f, dfdv] = massbalfunc(composition, ki, vapor_frac);
+    %New value for variable vapor_frac using Newton's method:
+    vapor_frac=vapor_frac-f/dfdv;
+    %physical constraint on vapor_frac
+%     ss = 0.7;
+%     dvf = f/dfdv;
+%     while ((vapor_frac_new<0) || (vapor_frac_new>1))
+%         dvf = ss*dvf;        
+%         vapor_frac_new=vapor_frac-dvf;
+%     end
+%     vapor_frac = vapor_frac_new;
+    if (vapor_frac<0)
+      vapor_frac=0;
+    elseif (vapor_frac>1)
+      vapor_frac=1;
+    end
+    
+    [liquid_x, vapor_y] = xy_calc(composition, vapor_frac, ki);
+    
+    error1=abs(sum(liquid_x)-1);
+    error2=abs(sum(vapor_y)-1);
+    %Another if statement based on my experience:
+    if (abs(dfdv)<eps1)
+        error3=f;
+    else                  
+        error3=abs(f/dfdv);
+    end
+    % normalize the mole fractions
+    liquid_x = mynormalize(liquid_x);
+    vapor_y = mynormalize(vapor_y);
+    
+    ki = kvalueLLE(mixture, thermo, liquid_x, vapor_y);
+    end
+  
+  % Final physical constraint for 1-phase result
+  if ((vapor_frac==0) || (vapor_frac==1))
+      vapor_y=composition;
+      liquid_x=composition;
+  end
+  vapor_y = mynormalize(vapor_y);
+  liquid_x = mynormalize(liquid_x);
diff --git a/Flash/stabilityLLETest.m b/Flash/stabilityLLETest.m
new file mode 100644
index 0000000..382212c
--- /dev/null
+++ b/Flash/stabilityLLETest.m
@@ -0,0 +1,164 @@
+function [stability_flag, SL, SV] = stabilityTest(mixture, thermo, varargin)
+% Michelsen stability test; I have used the algorithm described here:
+% https://www.e-education.psu.edu/png520/m17_p7.html
+% 
+% SYNOPSIS:
+%   
+% 
+% PARAMETERS:
+%   
+% 
+% RETURNS:
+%   
+% 
+% EXAMPLE:
+% 
+% SEE ALSO:
+%     
+
+%{
+Copyright (c) 2012, 2013, Ali Akbar Eftekhari
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or 
+without modification, are permitted provided that the following 
+conditions are met:
+
+    *   Redistributions of source code must retain the above copyright notice, 
+        this list of conditions and the following disclaimer.
+    *   Redistributions in binary form must reproduce the above 
+        copyright notice, this list of conditions and the following 
+        disclaimer in the documentation and/or other materials provided 
+        with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
+AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, 
+THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR 
+PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR 
+CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, 
+EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, 
+PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, 
+OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 
+LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 
+NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+%}
+
+% extract the data from flash options
+if nargin>2
+    trivial_eps = varargin{1}.trivialSolutionMaxError;
+    convergence_eps = varargin{1}.convergenceMaxError;
+    max_itr = varargin{1}.maxIteration;
+else
+    % default values
+    trivial_eps = 1e-5;
+    convergence_eps = 1e-10;
+    max_itr = 50;
+end
+
+% extract EOS function
+eosf = thermo.EOS;
+
+% switch on the fugacity calculation in thermo structure
+thermo.fugacity_switch=1; % switch on
+
+% initialize pseudo second phases
+gasmixture = mixture;
+liquidmixture = mixture;
+gasthermo = thermo;
+gasthermo.phase = 1;
+liquidthermo = thermo;
+liquidthermo.phase = 1;
+
+
+% extract the total composition and pressure
+composition = mixture.mole_fraction;
+p = mixture.pressure;
+
+% --------------------------- FIRST TEST ----------------------------------
+% calculate the fugacity of the mixture, assuming it is a liquid
+[~,~,fug_coef,~]=eosf(mixture, liquidthermo);
+mixfug = fug_coef.*composition*p; %[Pa]
+
+% Initial estimate for k-values using an empirical equation
+ki = kval_estimate(mixture);
+
+% assign large number to error values to begin the loop
+conv_error = 1+convergence_eps;
+triv_error = 1+trivial_eps;
+j = 0;
+while (conv_error>convergence_eps) && (triv_error>trivial_eps) && (j<max_itr)
+    j = j+1; % loop counter
+    % create a vapor-like second phase
+    Yi = composition.*ki;
+    SV = sum(Yi);
+
+    % normalize the vapor-like mole fractions
+    yi = Yi/SV;
+
+    % calculate the fugacity of the vapor-like phase using the thermo structure
+    gasmixture.mole_fraction = yi;
+    [~,~,fug_coef,~]=eosf(gasmixture, gasthermo);
+    gasfug = fug_coef.*yi*p; %[Pa]
+
+    % correct K-values
+    Ri = mixfug./gasfug*(1/SV);
+    ki = ki.*Ri;
+
+    % calculate the convergence and trivial solution error values
+    conv_error = sum((Ri-1).^2);
+    triv_error = sum(log(ki(ki>0)).^2);
+end
+
+% analyze the first test results
+if (conv_error>convergence_eps)
+    stability_flag(1) = 1; % converged to trivial solution
+elseif (triv_error>trivial_eps)
+    stability_flag(1) = 2; % converged
+elseif (j>=max_itr)
+    stability_flag(1) = 3; % maximum iteration reached
+end
+
+% --------------------------- SECOND TEST ---------------------------------
+% calculate the fugacity of the mixture, assuming it is a gas
+[~,~,fug_coef,~]=eosf(mixture, gasthermo);
+mixfug = fug_coef.*composition*p; %[Pa]
+
+% Initial estimate for k-values using an empirical equation
+ki = kval_estimate(mixture);
+
+% assign large number to error values to begin the loop
+conv_error = 1+convergence_eps;
+triv_error = 1+trivial_eps;
+j = 0;
+while (conv_error>convergence_eps) && (triv_error>trivial_eps) && (j<max_itr)
+    j = j+1; % loop counter
+    % create a liquid-like second phase
+    Xi = composition./ki;
+    SL = sum(Xi);
+
+    % normalize the liquid-like mole fractions
+    xi = Xi/SL;
+
+    % calculate the fugacity of the liquid-like phase using the thermo structure
+    liquidmixture.mole_fraction = xi;
+    [~,~,fug_coef,~]=eosf(liquidmixture, liquidthermo);
+    liquidfug = fug_coef.*xi*p; %[Pa]
+
+    % correct K-values
+    Ri = liquidfug./mixfug*SL;
+    ki = ki.*Ri;
+
+    % calculate the convergence and trivial solution error values
+    conv_error = sum((Ri-1).^2);
+    triv_error = sum(log(ki(ki>0)).^2);
+end
+
+% analyze the first test results
+if (conv_error>convergence_eps)
+    stability_flag(2) = 1; % converged to trivial solution
+elseif (triv_error>trivial_eps)
+    stability_flag(2) = 2; % converged
+elseif (j>=max_itr)
+    stability_flag(2) = 3; % maximum iteration reached
+end