this came up again in a suggested answer on stackoverflowit seems that matlab returns the results for multi-varialbe inputs to vpasolve as a struct with a field for each variable containing an ordered array of answers:

>> syms x1 x2
eqn1= -x1+2.*x1.^3+x2==0;
eqn2= -x1-x2==0;
>> sol=vpasolve([eqn1,eqn2],[x1,x2],[0,0])

sol = 

  struct with fields:

    x1: [3×1 sym]
    x2: [3×1 sym]

>> sol.x1
 
ans =
 
   0
 1.0
-1.0
 
>> sol.x2
 
ans =
 
   0
-1.0
 1.0

noted that there's already an open report about vpasolve not finding all roots that is unlikely to be fixed soon but output format matching may be worth doing for compatibility even if the data was incomplete