Flash calculations using the Soave-Redlich-Kwong equation of state do not agree with Aspen Plus results

Consider a phenomenon in which a four-component solvent mixture splits into gas and liquid phases at a certain temperature and pressure. We would like to find the composition of the components of the gas and liquid phases in this case. Inlet composition, temperature and pressure are as follows: z=[0.015,0.025, 0.17, 0.79], Temp=473.15[K], Pres=0.5[MPa]. When I ran the program that I will put at the end of the question, the composition of the gas and liquid phases came up without error, but the values were different from the results I calculated using Aspen Plus. I am having trouble finding the error even after reviewing the program many times. Could you please tell me what is wrong? For your information, the calculation formula is based on the following paper: Giorgio Soave, "Equilibrium Constants From a Modified Redlich-Kwong Equation of State", Article in Chemical Engineering Science · June 1972

<Program>

% 
load('parameter','wtm','T_c','P_c','omega_p'); % T_c[K] P_c[bar]

ncomp=4; % Number of components
Pres=0.5; % [MPa]
Temp=473.15; % [K]
z(1:ncomp,1)=[0.015,0.025, 0.17, 0.79];

% Objective Function Definition
objective = @(x) dif_fg(x,ncomp,Pres,Temp,z); % x=[x1,x2,x3,x4,q]
A = [];
b = [];
Aeq = [1 1 1 1 0];
beq = 1;
lb = zeros(5,1); % lower value
ub = ones(5,1); % upper value
x0=[0.05,0.05,0.15,0.75,0.5];

% optimization
% options = optimoptions('fmincon', 'Display', 'iter', 'Algorithm', 'sqp');
[x_opt, fval] = fmincon(objective, x0, A, b, Aeq, beq, lb, ub,[]);

function [dif_fg]=dif_fg(x1,ncomp,Pres,Temp,z)

load('parameter','T_c','P_c','omega_p'); % T_c[K] P_c[bar]

Pres=Pres*10; % [MPa]→[bar]
x = zeros(ncomp, 1);
y = zeros(ncomp, 1);
Tr = zeros(ncomp, 1);
m = zeros(ncomp, 1);
alpha = zeros(ncomp, 1);

x(1:ncomp,1)=x1(1,1:ncomp);
q=x1(1,1+ncomp);
y(1:ncomp,1)=(z(1:ncomp,1)-q*x(1:ncomp,1))/(1-q);

Tr(1:ncomp,1)=Temp./T_c(1:ncomp,1); % critical temperature
m(1:ncomp,1)=0.48+1.574*omega_p(1:ncomp,1)-0.176*omega_p(1:ncomp,1).^2;
alpha(1:ncomp,1)=(1+m(1:ncomp,1).*(1-sqrt(Tr(1:ncomp,1)))).^2;

% Liquid
A_L=0.42747*Pres/Temp^2*(sum(x(1:ncomp,1).*T_c(1:ncomp,1).*sqrt(alpha(1:ncomp,1))./sqrt(P_c(1:ncomp,1))))^2;
B_L=0.08664*Pres/Temp*sum(x(1:ncomp,1).*T_c(1:ncomp,1)./P_c(1:ncomp,1));

% Gas
A_V=0.42747*Pres/Temp^2*(sum(y(1:ncomp,1).*T_c(1:ncomp,1).*sqrt(alpha(1:ncomp,1))./sqrt(P_c(1:ncomp,1))))^2;
B_V=0.08664*Pres/Temp*sum(y(1:ncomp,1).*T_c(1:ncomp,1)./P_c(1:ncomp,1));

% Compressibility factor calculation (solving cubic equations)
% Liquid
coef_L=[1 -1 A_L-B_L-B_L^2 -A_L*B_L];
roots_z_L = roots(coef_L);
real_roots_z_L = roots_z_L(imag(roots_z_L) == 0); % Extract only real solutions

% The smallest value of the real solution is the compressibility factor of the liquid phase
if ~isempty(real_roots_z_L)
    z_L = min(real_roots_z_L);
else
    z_L = NaN; % there is no real solution
end

% Gas
coef_V=[1 -1 A_V-B_V-B_V^2 -A_V*B_V];
roots_z_V = roots(coef_V);
real_roots_z_V = roots_z_V(imag(roots_z_V) == 0); % Extract only real solutions

% The largest value of the real solution is the compressibility factor of the gas phase
if ~isempty(real_roots_z_V)
    z_V = max(real_roots_z_V);
else
    z_V = NaN; % there is no real solution
end

% Fugacity Calculations
% Liquid
a_al_L(1:ncomp,1)=sqrt(alpha(1:ncomp,1)).*T_c(1:ncomp,1)./sqrt(P_c(1:ncomp,1))/sum(x(1:ncomp,1).*T_c(1:ncomp,1).*sqrt(alpha(1:ncomp,1))./sqrt(P_c(1:ncomp,1)));
b_b_L(1:ncomp,1)=T_c(1:ncomp,1)./P_c(1:ncomp,1)/sum(x(1:ncomp,1).*T_c(1:ncomp,1)./P_c(1:ncomp,1));

ln_phi_L(1:ncomp,1)=b_b_L(1:ncomp,1)*(z_L-1)-log(z_L-B_L)-A_L/B_L*(2*a_al_L(1:ncomp,1)-b_b_L(1:ncomp,1))*log((z_L+B_L)/z_L);

phi_L(1:ncomp,1)=exp(ln_phi_L(1:ncomp,1)); % fugacity coefficient

fg_L(1:ncomp,1)=Pres*phi_L(1:ncomp,1).*x(1:ncomp,1); % fugacity

% Gas
a_al_V(1:ncomp,1)=sqrt(alpha(1:ncomp,1)).*T_c(1:ncomp,1)./sqrt(P_c(1:ncomp,1))/sum(y(1:ncomp,1).*T_c(1:ncomp,1).*sqrt(alpha(1:ncomp,1))./sqrt(P_c(1:ncomp,1)));
b_b_V(1:ncomp,1)=T_c(1:ncomp,1)./P_c(1:ncomp,1)/sum(y(1:ncomp,1).*T_c(1:ncomp,1)./P_c(1:ncomp,1));

ln_phi_V(1:ncomp,1)=b_b_V(1:ncomp,1)*(z_V-1)-log(z_V-B_V)-A_V/B_V*(2*a_al_V(1:ncomp,1)-b_b_V(1:ncomp,1))*log((z_V+B_V)/z_V);

phi_V(1:ncomp,1)=exp(ln_phi_V(1:ncomp,1)); % fugacity coefficient

fg_V(1:ncomp,1)=Pres*phi_V(1:ncomp,1).*y(1:ncomp,1); % fugacity

dif_fg=sum((fg_L(1:ncomp,1)-fg_V(1:ncomp,1)).^2);

end