MATLAB: Matrix dimensions must agree. Error Message appears.

errorinputmatrixmatrix array

Hi
I'm running a function to calculate fugacity and other ouput. The function works some of the time. Like if i run it with the input below:
Which is a gas mixture of 3 components-.
``Pc=[4.59E6 3.36E6 2.47E6] Tc=[190.6 469.7 568.7]w=[0.011 0.251 0.396]kappa=zeros(3)eta=zeros(3)P=0.1E6T=300z=[0.85 0.09 0.06]Ko=[332 0.75 0.03]``
But when i try to run a gas mixture of 9 components, with this input:
``Pc=[4.599E6 4.872E6 4.248E6 3.796E6 3.37E6 3.025E6 2.49E6 1.817E6 1.401E6]; Tc=[190.6 305.3 369.8 425.1 469.7 507.6 568.7 658.1 722];T=320;Omega=[0.012 0.1 0.152 0.2 0.252 0.301 0.346 0.577 0.721];P=0.2E6;kappa=zeros(9);eta=zeros(9);ki=(Pc/P.*exp(5.37.*(1+Omega).*(1-Tc/T)));Ko=ki;z= [0.7 0.08 0.07 0.03 0.01 0.02 0.04 0.02 0.03];``
this error message Appears :
``Matrix dimensions must agree.Error in PTFLASH_VLE_PRVDW (line 28)    while max(abs((fV-fL)./fV))>0.00001 && iter<1000 && psi_0*psi_1<    end``
The complete Code is below with different subprograms.
The function:
``function [x,y,alphaV,fL,fV] = PTFLASH_VLE_PRVDW(Pc, Tc, Omega, kappa, eta, P, T, z, Ko)%The funtion PTFLASH_VLE_PRVDW makes a PT-FLASH calculation for a mixture%at given P, T and overall composition (z) using PR-CEoS and VDW MIXING RULES.%This function requires MIXRULES_VDW and FUGACITY_INMIX_PRVDW funtions to run.%All input/output data are expressed in SI.%P[Pa], T[K], w[dimensionless], V[m^3/mol], Z[dimensionless]%giving a firt value for cc=length(z);%firt verification for alphaV between 0 and 1psi_0=sum(z.*(Ko-1));psi_1=sum(z.*(Ko-1)./Ko);if psi_0*psi_1>=0    x=0;    y=0;    alphaV=0;    fL=0;    fV=0;    else    iter=0;    fL=zeros(1,c);    fV=[1,1,1];    K=Ko;        while max(abs((fV-fL)./fV))>0.00001 && iter<1000 && psi_0*psi_1<0        [ alphaV ] = RACHFORDRICE_BISECT( K,z );        x=z./((1-alphaV)+alphaV*K);        y=K.*x;        [fiL,~,fL,~] = fugacity_INMIX_PRVDW( Pc,Tc,Omega,kappa,eta,P,T,x );        [~,fiV,~,fV] = fugacity_INMIX_PRVDW( Pc,Tc,Omega,kappa,eta,P,T,x );        K=fiL./fiV;        psi_0=sum(z.*(K-1));        psi_1=sum(z.*(K-1)./K);        iter=iter+1;          endendend``
The subprograms:
``function [fiL,fiV,fL,fV] = FUGACITY_INMIX_PRVDW(Pc, Tc, w, kappa, eta, P, T, x)%The funtion FUGACITY_INMIX_PRVDW calculates the fugacity (f) and fugacity %coefficient (fi) for a fluid IN A MIXTURE at given pressure P, temperature T %and composition using PR-CEoS and VDW_MIXING RULES. This function requires %ZMIX_PRVDW funtion to run.%All input/output data are expressed in SI.    %P[Pa], T[K], f[Pa], fi[dimensionless]%universal gas constantR=8.314;%PR attractive constant and covolume at critical pointb=0.07780*R*Tc./Pc;k=0.37464+1.54226*w-0.26992*w.^2;alpha=(1+k.*(1-sqrt(T./Tc))).^2;a=0.45724*alpha.*(R.^2*(Tc.^2))./Pc;%calling MIXRULES_VDW function[am, bm, aa] = MIXRULES_VDW(a, b, kappa, eta, x);%calculation of Am and BmAm=am*P/(R*T)^2;Bm=bm*P/(R*T);%calculation of A and BA=aa*P/(R*T)^2;B=b*P/(R*T);%calculation of ZL and ZV using ZMIX_PRVDW function[Z,V,ZL,ZV,VL,VV] = ZMIX_PRVDW(Pc, Tc, w, kappa, eta, P, T, x);%giving a firt value for cc=length(x);%calculation of fugacity coeffiecientfor i=1:c    c(i)=x*(A(i,:))';    lnfiL(i)=B(i)/Bm*(ZL-1)-log(ZL-Bm)-Am/(2*sqrt(2)*Bm)*(2*c(i)/Am-B(i)/Bm)*log((ZL+(1+sqrt(2))*Bm)/(ZL+(1-sqrt(2))*Bm));    lnfiV(i)=B(i)/Bm*(ZV-1)-log(ZV-Bm)-Am/(2*sqrt(2)*Bm)*(2*c(i)/Am-B(i)/Bm)*log((ZV+(1+sqrt(2))*Bm)/(ZV+(1-sqrt(2))*Bm));    fiL(i)=exp(lnfiL(i));    fiV(i)=exp(lnfiV(i));end%calculation of fugacity fL=P*fiL.*x;fV=P*fiV.*x;end``
Another subprogram:
``function [Z,V,ZL,ZV,VL,VV] = ZMIX_PRVDW(Pc, Tc, w, kappa, eta, P, T, x)%The funtion ZMIX_PRVDW calculates the compressibility factor z and the specific%molar volume V for a mixture at given pressure P, temperature T and concentration%using PR-CEoS and VDW MIXING RULES. This function requires REALROOTS_CARDANO and% MIXRULES_VDW funtions to run.%All input/output data are expressed in SI.%P[Pa], T[K], w[dimensionless], V[m^3/mol], Z[dimensionless]%universal gas constantR=8.314;%PR attractive constant and covolume at critical pointb=0.07780*(R*Tc)./Pc;k=0.37464+1.54226*w-0.26992*w.^2;alpha=(1+k.*(1-sqrt(T./Tc))).^2;a=0.45724*alpha.*(R^2.*(Tc.^2))./Pc;%calling MIXRULES_VDW function[am, bm, aa] = MIXRULES_VDW(a, b, kappa, eta, x);%calculation of Am and BmAm=am*P/(R*T)^2;Bm=bm*P/(R*T);%calculation of A and Ba2=-1+Bm;a1=Am-3*Bm^2-2*Bm;a0=-Am*Bm+Bm^2+Bm^3;%calculation of Z and V using REALROOTS_CARDANO functionZ=REALROOTS_CARDANO(a2,a1,a0);V=Z*R*T/P;%selection of Z and V values for liquid and vaporif min(V)<b    ZL=max(Z);    ZV=max(Z);    VL=max(V);    VV=max(V);else    ZL=min(Z);    ZV=max(Z);    VL=min(V);    VV=max(V);    endend``
The last Subprogram:
``function [RR]=REALROOTS_CARDANO(a2,a1,a0)%This function calculates the real roots of a polynomial of degree 3 using %an algorithm based on the Cardano's formula. %The polynomial must be written as: x^3 + a2*x^2 + a1*x + a0. %Input data: a2, a1, a0.%Output data: array containing the real roots. %The output array can be 1x1 or 1x3. Q=(a2^2-3*a1)/9;R=(2*a2^3-9*a1*a2+27*a0)/54;    if Q==0 && R==0 %Case of single real root with multiplicity 3        RR=-a2/3;    else        if Q^3-R^2>=0             theta=acos(R/(Q^(3/2)));            RR(1)=-2*sqrt(Q)*cos(theta/3)-a2/3;            RR(2)=-2*sqrt(Q)*cos((theta+2*pi)/3)-a2/3;            RR(3)=-2*sqrt(Q)*cos((theta+4*pi)/3)-a2/3;        else            RR=-sign(R)*((sqrt(R^2-Q^3)+abs(R))^(1/3)+Q/(sqrt(R^2-Q^3)+abs(R))^(1/3))-a2/3;        end    endend``

``    fL=zeros(1,c);    fV=[1,1,1];``
``      fV=ones(1,c)``