Large system of equations, vpaSolve returning empty sym

I have to do a project for my chemical engineering thermodynamics course where I'm attempting to find the equilibrium vapor and liquid mole fractions of a binary mixture. It involves a large system of equations, and when I go to solve them I get an empty struct and an empty sym for all the variables. I'm not sure where my mistake is, any help is much appreciated.

%Project 5
clc; clearvars; format short
R=83.14; %bar cm^3/ mol K
P=40/.986923; %bar
w=[0.224 0.094]; %accentric factors
Tc=[304.2 373.5]; %Critical temperatures, K
Pc=[73.83 89.63]; %Critical pressures, bar
 
%Peng-Robinson Parameters
sigma=1+sqrt(2);
eps=1-sqrt(2);
omega=0.07780;
psi=0.45724;
Zpr=0.30740;
Pr=P./Pc;
T=280;
Tr1=T./(Tc(1));
Tr2=T./(Tc(2));
%Mixture information

%Reproducing T-xy diagram of carbon dioxide(1)/hydrogen sulfide(2)
k12=0;
alpha1=(1+(0.37464 +1.54226*w(1)-0.26992*w(1)^2)*(1-Tr1^.5))^2;
alpha2=(1+(0.37464 +1.54226*w(2)-0.26992*w(2)^2)*(1-Tr2^.5))^2;
a1=psi*((alpha1*R^2*Tc(1)^2)/Pc(1));
a2=psi*((alpha2*R^2*Tc(2)^2)/Pc(2));
a12=(1-k12)*sqrt(a1*a2);
b=omega.*((R.*Tc)./Pc);

syms x1 y1 aliq bliq avap bvap Betaliq qliq aBar1Liq aBar2Liq qBar1Liq qBar2Liq Betavap qvap aBar1vap aBar2vap qBar1vap qBar2vap Zvap phivap1 phivap2 Ivap Zliq Iliq philiq1 philiq2
%Mixing and combining rules, liquid phase
eq(1)=aliq==(x1^2*a1)+(2*x1*(1-x1)*a12)+(((1-x1)^2)*a2);
eq(2)=bliq==x1*b(1)+(1-x1)*b(2);
%Mixing and combining rules, vapor phase
eq(3)=avap==(y1^2*a1)+(2*y1*(1-y1)*a12)+((1-y1)^2*a2);
eq(4)=bvap==y1*b(1)+(1-y1)*b(2);
%Dimensionless parameters, liquid phase
eq(5)=Betaliq==(bliq*P)/(R*T);
eq(6)=qliq==aliq/(bliq*R*T);
eq(7)=aBar1Liq==2*x1*a1+2*(1-x1)*a12-aliq;
eq(8)=aBar2Liq==2*(1-x1)*a2+2*x1*a12-aliq;
eq(9)=qBar1Liq==qliq*(1+(aBar1Liq/aliq)-(b(1)/bliq));
eq(10)=qBar2Liq==qliq*(1+(aBar2Liq/aliq)-(b(2)/bliq));
%Dimensionless parameters, vapor phase
eq(11)=Betavap==(bvap*P)/(R*T);
eq(12)=qvap==avap/(bvap*R*T);
eq(13)=aBar1vap==2*y1*a1+2*(1-y1)*a12-avap;
eq(14)=aBar2vap==2*(1-y1)*a2+2*y1*a12-avap;
eq(15)=qBar1vap==qvap*(1+(aBar1vap/avap)-(b(1)/bvap));
eq(16)=qBar2vap==qvap*(1+(aBar2vap/avap)-(b(2)/bvap));
%Equilibrium and related equations
%Liquid phase
eq(17)=Zliq==Betaliq+(Zliq+eps*Betaliq)*(Zliq+sigma*Betaliq)*((1+Betaliq-Zliq)/(qliq*Betaliq));
eq(18)=Iliq==1/(sigma-eps)*log((Zliq+eps*sigma)/(1+eps*Betaliq));
eq(19)=philiq1==exp(((b(1)/bliq)*(Zliq-1))-log(Zliq-Betaliq)-qBar1Liq*Iliq);
eq(20)=philiq2==exp(((b(2)/bliq)*(Zliq-1))-log(Zliq-Betaliq)-qBar2Liq*Iliq);
%Vapor phase
eq(21)=Zvap==1+Betavap-(qvap*Betavap)*((Zvap-Betavap)/((Zvap+eps*Betavap)*(Zvap+sigma*Betavap)));
eq(22)=Ivap==1/(sigma-eps)*log((Zvap+eps*sigma)/(1+eps*Betavap));
eq(23)=phivap1==exp(((b(1)/bvap)*(Zvap-1))-log(Zvap-Betavap)-qBar1vap*Ivap);
eq(24)=phivap2==exp(((b(2)/bvap)*(Zvap-1))-log(Zvap-Betavap)-qBar2vap*Ivap);
%Equifugacity Equations
eq(25)=y1*phivap1==x1*philiq1;
eq(26)=(1-y1)*phivap2==(1-x1)*philiq2;
% [Aliq, Bliq, Avap, Bvap, BetaLiq, Qliq, ABar1Liq, ABar2Liq, QBar1Liq, QBar2Liq, BetaVap, Qvap, ABar1vap, ABar2vap, QBar1vap, QBar2vap, ZVap, Phivap1, Phivap2, IVap, ZLiq, ILiq, Philiq1, Philiq2, X1, Y1]
S=vpasolve(eq,[aliq bliq avap bvap Betaliq qliq aBar1Liq aBar2Liq qBar1Liq qBar2Liq Betavap qvap aBar1vap aBar2vap qBar1vap qBar2vap Zvap phivap1 phivap2 Ivap Zliq Iliq philiq1 philiq2 x1 y1],[0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf;0 inf]);
Xx1=eval(S.x1)
Yy1=eval(S.y1)

[John D'Errico]

Why are you so sure you made a mistake? Not all large systems of nonlinear equations have an analytical solution, and even if you search for a numerical solution, there is no assurance that a numerical solver (vpasolve) will find a solution.

Note that you can improve your chances of success by providing intelligent starting values for the solver to begin its search. vpasolve can use that, though there is no assurance of a solution.

Not every problem has a simple solution. Some are not solvable at all. An art of mathematics and modeling in general is to know how to reduce problem complexity so a difficult one becomes solvable, to find ways to make an intractable problem tractable. But that does not mean it has to be easy, or that computers can solve all problems. Computers only do that in the movies and on tv.