Issue Changing Script to function in MatLab

I am trying to give a/set 'rho_real' to this code and receive P_Mpa as an/set of outputs, however the initial 'rho_real' does not seem to get changed as I run the function.

Any Help is greatly appreciated.

function P_Mpa = Density_vs_Pressure_Ethane (rho_real)

% calculating Pressure vs Density curve using 

mi = [1.6069,1.6069]; 
sigmai = [3.5206,3.5206]; 1
% epsilon = [2.873268218e-21,2.873268218e-21]; 
epsilon_ki = [191.42,191.42]; 
xi = [0.5,0.5]; % mole fraction of component i
T = 298.15; % temperature of the system
% k_bolt = 1.3806488e-23; % Boltzmann constant = 1.3806488 × 10-23 m2 kg s-2 K-1

% Initial guess for density

rho_real = 10000;

rho = rho_real*6.022e-7;

% Initial Kij

Kij = 0;

% Temperatuure-dependent segment diameter di of component i - matrix save

di = zeros(1,2);

giihs = zeros (1,2);

rho_giihs = zeros (1,2);


m_bar = 0; 

zi0=0;
zi1=0;
zi2=0;
zi3=0;
for i=1:2
m_bar=m_bar+((xi(1,i))*(mi(1,i)));


di(1,i) = (sigmai(1,i))*(1-(0.12*exp(-3*(epsilon_ki(1,i))/T)));

% Calculating Zieeeeeeh for Hard Sphere compressibility (Rechcecked)

 zi0 = zi0+((pi/6)*rho)*(xi(1,i)*(mi(1,i))*((di(1,i))^0));
 zi1 = zi1+((pi/6)*rho)*(xi(1,i)*(mi(1,i))*((di(1,i))^1));
 zi2 = zi2+((pi/6)*rho)*(xi(1,i)*(mi(1,i))*((di(1,i))^2));
 zi3 = zi3+((pi/6)*rho)*(xi(1,i)*(mi(1,i))*((di(1,i))^3));
end

for i=1:2 

giihs (1,i)= (1/(1-zi3))+((((di(1,i))^2)/(2*(di(1,i))))*(3*zi2)/((1-zi3)^2))+((di(1,i)^2)/(2*(di(1,i)^2))^2)*((2*(zi2)^2)/((1-zi3)^3)); 

% Derevative of site to site radial distribution function of hard spheres
% with respect to density

rho_giihs (1,i) = ((zi3/((1-zi3)^2)))+(((((di(1,i))^2)/(2*(di(1,i)))))*...
((((3*zi2)/((1-zi3)^2)))+((6*zi2*zi3)/((1-zi3)^3))))+...
(((((di(1,i))^2)/(2*(di(1,i))))^2)*(((4*(zi2^2))/((1-zi3)^3))+...
(((6*zi2^2)*zi3)/((1-zi3)^4)))); 
end


eta = zi3;

% Calculating the hard sphere compresibility factor ( Rechecked)

zhs = (zi3/(1-zi3)+((3*zi1*zi2)/((zi0*(1-zi3)^2)))+(((3*(zi2^3))-(zi3*(zi2^3)))/(zi0*((1-(zi3)^3)))));

zhc=0;

for i = 1:2
% Calculating the hard chain compressibility factor (Rechecked)

zhc= zhc+((xi(1,i)*(mi(1,i)-1)*((giihs(1,i))^-1)*((rho_giihs(1,i)))));
end
  zhc=m_bar*zhs-zhc;

% a0 constants 

a0 = [0.9105631445,0.6361281449,2.6861347891,-26.547362491,97.759208784,-159.59154087,91.297774084]; % Checked and Correct Sadowski 2001

% a1 constants

a1 = [-0.3084016918,0.1860531159,-2.5030047259,21.419793629,-65.255885330,83.318680481,-33.746922930]; % Checked and Correct Sadowski 2001

% a2 constants

a2 = [-0.0906148351,0.4527842806,0.5962700728,-1.7241829131,-4.1302112531,13.776631870,-8.6728470368]; % Checked and Correct Sadowski 2001

% b0 constants

b0 = [0.7240946941,2.2382791861,-4.0025849485,-21.003576815,26.855641363,206.55133841,-355.60235612]; % Checked and Correct Sadowski 2001

% b1 constants

b1 = [-0.5755498075,0.6995095521,3.8925673390,-17.215471648,192.67226447,-161.82646165,-165.20769346]; % Checked and Correct Sadowski 2001

% b2 constants

b2 = [0.0976883116,-0.2557574982,-9.1558561530,20.642075974,-38.804430052,93.626774077,-29.666905585]; 

aim = zeros(1,7);
bim = zeros(1,7);

for i = 1:7
aim(1,i)=a0(1,i)+(((m_bar-1)/m_bar)*a1(1,i))+(((m_bar-1)/m_bar)*((m_bar-2)/m_bar))*a2(1,i); % Checked and Correct Sadowski 2001 (rechecked) 
bim(1,i)=b0(1,i)+(((m_bar-1)/m_bar)*b1(1,i))+(((m_bar-1)/m_bar)*((m_bar-2)/m_bar))*b2(1,i); 


c1 = ((1 + (m_bar*(((8*eta)-(2*eta^2))/(1-eta)^4)+((1-m_bar)*(((20*eta)-(27*eta^2)+(12*eta^3)-(2*eta^4))/(((1-eta)*(2-eta))^2)))))^-1); 

c2 = ((-c1^2)*(m_bar*(((-4*eta^2)+(20*eta)+8)/((1-eta)^5))+((1-m_bar)*(((2*eta^3)+(12*eta^2)-(48*eta)+40)/(((1-eta)*(2-eta))^3))))); 

% Integrals of perturbation theory without/with respect to eta (rechecked)
dI1 = 0;
dI2 = 0;
I1 = 0;
I2 = 0;
for j = 1:7

I1 = I1+ aim(1,j)*eta^(j-1); 
I2 = I2+ bim(1,j)*eta^(j-1); 

dI1 = dI1 + (aim(1,j)*((j-1)+1)*(eta^(j-1))); 
dI2 = dI2 + (bim(1,j)*((j-1)+1)*(eta^(j-1))); 
end
% Segment abriviations 

m2e = 0;
m2e2 = 0;
sigma_ij = zeros(1,2);
epsilon_ij = zeros (1,2);
for i = 1:2
for j = 1:2
    % Mixing rules for sigma and eta respectively.
    sigma_ij(i,j)=0.5*(sigmai(1,i)+sigmai(1,j)); 
    epsilon_ij(i,j)=sqrt((epsilon_ki(1,i)*epsilon_ki(1,j))*(1-Kij)); 

    m2e=m2e+(xi(1,i)*xi(1,j)*mi(1,i)*mi(1,j)*((epsilon_ij(i,j)/(T)))*(sigma_ij(i,j)^3)); 
    m2e2=m2e2+(xi(1,i)*xi(1,j)*mi(1,i)*mi(1,j)*(((epsilon_ij(i,j)/(T)))^2)*(sigma_ij(i,j)^3)); 
end

end



% The dispersion contribution to the comprehensibility factor 

zdis = (-2*pi*rho*dI1*m2e)-(pi*rho*m_bar*((c1*dI2)+(c2*eta*I2))*m2e2); 

% Compresibility Factor (rechecked)

z = 1 + zhc + zdis;

zdis;
zhc;

P = z*8.314*T*rho/6.022e-7; % (rechecked)
P_Mpa = P*1e-6;