Where exactly do you get the output 0.016 from the first snippet of code?
Do you mean the sigma2 value? If yes, then you are misunderstanding the solution. The solutions are the elements of the 3x1 vector, where the terms have been grouped. And the information below the solutions show, which terms have been grouped into the respective variables.
Also, there's a mistake in the 2nd piece of code, you have written a*x^3 + b*x^2 + c*x + b == 0, it should be d instead of b. Correcting the mistake leads to the correct output, see below -
clc,clear
syms x a
u = symunit;
a = 1.01325*u.bar;
b = -0.049576382*u.bar*u.m^3/u.mol;
c = .000003657553277*u.bar*u.m^6/(u.mol)^2;
d = -.00000000015676624*u.bar*u.m^9/u.mol^3;
eqn = a*x^3+b*x^2+c*x+b == 0
vpa(solve(eqn,x,"MaxDegree",3), 10)
clc,clear
syms x
u = symunit;
a = 1.01325;
b = -0.049576382;
c = 3.657553277*10^-6;
d = -1.5676624*10^-10;
% v
eqn = a*x^3+b*x^2+c*x+d == 0;
solve(eqn,x,"MaxDegree",3);
vpa(solve(eqn,x,'MaxDegree',3),10)
clear
a = 1.01325;
b = -0.049576382;
c = 3.657553277*10^-6;
d = -1.5676624*10^-10;
x = [a b c d];
vpa(roots(x),10)
u = symunit;
a = 1.01325*u.bar;
b = -0.049576382*u.bar*u.m^3/u.mol;
c = .000003657553277*u.bar*u.m^6/(u.mol)^2;
d = -.00000000015676624*u.bar*u.m^9/u.mol^3;
x = [a b c d];
vpa(roots(x),10)
