Function 'solve' is working in Matlab 2014a but not in Matlab 2018a, what's possibly wrong?

Hello, everyone!

I've been trying to run a code that I have to solve a 22 ode system in order to get the equillibrium points for the system. I'm using the function solve and I have tested my code in Matlab 2014a and it works, but when I want to run it in Matlab 2018a it gives me this error:

Error using sym/solve

Inconsistent output with 21 variables for input argument with 1 variables.

I'm not sure what's wrong with the code since it works in Matlab 2014a, and I think I understand what the error means, but I haven't figure out what's wrong with it. I checked previous questions about working with solve in Matlab 2018 or further version, and they are more about syntaxis than the error I'm having. I also just checked the code in Matlab Online, but it's the same as in Matlab 2018. If you could please help me out to figure this out, I'd be grateful. Here's the code:

syms x1 x2 x3 x4 x5 x6 x7 x9 x10 x11 x12 x13 x14 x15 x16 x17 x18 x19 x20 x21 x22

x8 = 180;%DVx8;

% Constants:
rBGU = 70;
rRBCU = 10;
rGGU = 20;
rPGUB = 35; 
rHGPB = 155;
TI = 25;
TGm = 65; 
rHGUB = 20;
FKIC = 0.30;
FPIC = 0.15;
FLIC = 0.40;
alpha = 0.0482;
beta = 0.931;
K = 0.00794;
M1 = 0.00747;
M2 = 0.0958;
gamma = 0.575;
Q0 = 6.33;
rMGmC = 9.1;
rPGmRB = 9.1;

% Volumenes:
VBVG = 3.5;
VBIG = 4.5;
VHG = 13.8;
VLG = 25.1;
VGG = 11.2;
VKG = 6.6;
VPVG = 10.4;
VPIG = 67.4;
VBI = 0.26;
VHI = 0.99;
VGI = 0.94;
VLI = 1.14;
VKI = 0.51;
VPVI = 0.74;
VPII = 6.74;
VGm = 11310; 


QBG = 5.9;
QHG = 43.7;
QAG = 2.5;
QLG = 12.6;
QGG = 10.1; 
QKG = 10.1;
QPG = 15.1;
QBI = 0.45;
QHI = 3.12;
QAI = 0.18;
QKI = 0.72;
QPI = 1.05;
QGI = 0.72;
QLI = 0.90;

TB = 2.1;
TPG = 5;
TPI = 20;


IPIB = 4.5; 
GPIB = 117;

GmB = 1;

ILB = 21.43;

GLB = 101;
GHB=91.89;
IHB = 15.15;


MPGUI = 7.03 + 6.52*tanh(0.338*((x15/IPIB) - 5.82));
MPGUG = x8/GPIB;
rPGU = MPGUI*MPGUG*rPGUB;

MHGPGm = 2.7*tanh(0.39*(x19/GmB)) - x18;
MHGPG = 1.42 - 1.41*tanh(0.62*((x5/GLB) - 0.497));
rHGP = x16*MHGPGm+MHGPG*rHGPB;

MHGPIi = 1.21 - 1.14*tanh(1.66*((x12/ILB) - 0.89));

MHGUG = 5.66 + 5.66*tanh(2.44*((x5/GLB) - 1.48));
MHGUIi = 2*tanh(0.55*(x12/ILB));
rHGU = x17*MHGUG*rHGUB;

% if 0<=x(6)<460
    rKGE = 71 + 71*tanh(0.11*(x6 - 460));
% else if x(6)>460
%         rKGE = -330 + 0.872*x(6);
%     end
% end

rPIRB=4.5;
XB=0.1928; 
X = (x3^3.27)/((132^3.27) + 5.93*(x3^3.02));
Y = X^1.11;
Pinf = Y;

% if  X > x(22)
    S = ((M1*Y) + M2*(X-x22))*x20;
% else if X <= x(22)
%         S = M1*Y*x20;
%     end
% end

YB=XB^1.11;  
SB=M1*YB*x20;
rPIR = (S/SB)*rPIRB;

rLIC = FLIC*(QAI*x10 + QGI*x11 + rPIR);
rKIC = FKIC*(QKI*x13);
rPIC = x15/(((1 - FPIC)/FPIC)*(1/QPI) - (TPI/VPII));

rPGmC = rMGmC*x19;

MPGmRG = 2.93 - 2.10*tanh(4.18*((x3/GHB) - 0.61));
MPGmRI = 1.31 - 0.61*tanh(1.06*((x10/IHB) - 0.47));
rPGmR = MPGmRG*MPGmRI*rPGmRB;

Gempt = 0;
I = 0;

dxdt = [(QBG*(x3 - x1) - (VBIG/TB)*(x1 - x2))/VBVG;     % x(1)  dGBV   
        ((VBIG/TB)*(x1 - x2) - rBGU)/VBIG;                   % x(2)  dGBI
        (QBG*x1 + QLG*x5 + QKG*x6 + QPG*x7 - QHG*x3 - rRBCU)/VHG;  % x(3) dGH
        (QGG*(x3 - x4) - rGGU)/VGG;                        % x(4) dGG
        (QAG*x3 + QGG*x4 - QLG*x5 + rHGP - rHGU)/VLG;    % x(5) dGL
        (QKG*(x3 - x6) - rKGE)/VKG;                        % x(6) dGK
        (QPG*(x3 - x7) - (VPIG/TPG)*(x7 - x8))/VPVG;    % x(7) dGPV
        ((VPIG/TPG)*(x7 - x8) - rPGU + Gempt)/VPIG;                 % x(8) dGPI
        (QBI*(x10 - x9))/VBI;                              % x(9) dIB      
        (QBI*x9 + QLI*x12 + QKI*x13 + QPI*x14 - QHI*x10)/VHI;   % x(10) dIH
        (QGI*(x10 - x11))/VGI;                             % x(11) dIG
        (QAI*x10 + QGI*x11 - QLI*x12 + rPIR - rLIC)/VLI; % x(12) dIL
        (QKI*(x10 - x13) - rKIC)/VKI;                      % x(13) dIK
        (QPI*(x10 - x14) - (VPII/TPI)*(x14 - x15))/VPVI; % x(14) dIPV
        ((VPII/TPI)*(x14 - x15) - rPIC + I)/VPII;               % x(15) dIPI
        (1/TI)*(MHGPIi - x16);                               % x(16) dM_HGPI   
        (1/TI)*(MHGUIi - x17);                               % x(17) dM_HGUI 
        (1/TGm)*(((MHGPGm - 1)/2) - x18);                     % x(18) df2
        (rPGmR - rPGmC)/VGm;                                   % x(19) dGm     
        K*(Q0 - x20) + gamma*x21 - S;                      % x(20) dQ      
        alpha*(Pinf - x21);                                  % x(21) dP
        beta*(X - x22)                                       % x(22) dI
        ];
    
%[x1,x2,x3,x4,x5,x6,x7,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22]=solve(dxdt(1)==0,x1,dxdt(2)==0,x2,dxdt(3)==0,x3,dxdt(4)==0,x4,dxdt(5)==0,x5,dxdt(6)==0,x6,dxdt(7)==0,x7,dxdt(9)==0,x9,dxdt(10)==0,x10,dxdt(11)==0,x11,dxdt(12)==0,x12,dxdt(13)==0,x13,dxdt(14)==0,x14,dxdt(15)==0,x15,dxdt(16)==0,x16,dxdt(17)==0,x17,dxdt(18)==0,x18,dxdt(19)==0,x19,dxdt(20)==0,x20,dxdt(21)==0,x21,dxdt(22)==0,x22)
[x1,x2,x3,x4,x5,x6,x7,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22]=solve([dxdt(1)==0,dxdt(2)==0,dxdt(3)==0,dxdt(4)==0,dxdt(5)==0,dxdt(6)==0,dxdt(7)==0,dxdt(9)==0,dxdt(10)==0,dxdt(11)==0,dxdt(12)==0,dxdt(13)==0,dxdt(14)==0,dxdt(15)==0,dxdt(16)==0,dxdt(17)==0,dxdt(18)==0,dxdt(19)==0,dxdt(20)==0,dxdt(21)==0,dxdt(22)==0],[x1,x2,x3,x4,x5,x6,x7,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22])

Warning: Unable to solve symbolically. Returning a numeric solution using vpasolve.

x1 = 

152.51746839000882498782078251139

x2 = 

119.85080172334215832115411584472

x3 = 

164.38187516966984193697332488427

x4 = 

162.40167714986786173895352290407

x5 = 

161.90306860896015758279915028603

x6 = 

164.38187516966984193697332488347

x7 = 

171.74829653820904874906568249659

x9 = 

10.52856143188951116968952089685

x10 = 

10.52856143188951116968952089685

x11 = 

10.52856143188951116968952089685

x12 = 

14.314794100665162301866329373217

x13 = 

8.0988934091457778228380929975767

x14 = 

8.9492772171060844942360927623222

x15 = 

4.0286587437452595410429190494029

x16 = 

1.6121084800737452314752089913977

x17 = 

0.70341080109077170697095000731877

x18 = 

x19 = 

0.96906627607153986571161022353814

x20 = 

26.638894794210821344927810172986

x21 = 

0.42885485950148121688388992118223

x22 = 

0.46638873285385637550428617200014

X = [x1,x2,x3,x4,x5,x6,x7,x8,x9,x10,x11,x12,x13,x14,x15,x16,x17,x18,x19,x20,x21,x22];
% assignin('base','X',X); 

X = vpa(X,4);

[Torsten]

Try the "solve" command in the way I changed it in your code.

[Annie]

Thanks so much! It worked, I hadn't found that way of syntax for solve, thanks