You can use the vpasolve function to solve the polynomial, use the sym2poly function to automate extracting the vector of coefficients to a double precision vector, or use the coeffs function to automate extracting the vector of coefficients to a symbolic vector.
%%
syms P
sympref('FloatingPointOutput', false);
% Matrix definitions
KE_red =[13.0000 0 -8.4853;
0 13.0000 0;
-8.4853 0 8.0000];
KG_red =[ -(3*2^(1/2)*P)/65, 0, P/130;
0, -(3*2^(1/2)*P)/65, 0
P/130, 0, -(2*2^(1/2)*P)/195];
% Determinant
D=det(KE_red + KG_red)
% Solve for the roots of the determinant = 0
p = [-9*2^(1/2)/219700, 167/3380, -1328*2^(1/2)/195, 416];
r = roots(p)
s = vpasolve(D)
p2 = sym2poly(D)
r2 = roots(p2)
p3 = coeffs(D, 'all') % Handle the case where one of the powers is not present
r3 = roots(p3)
To check, let's sort the various vectors. I'll use vpa to approximate the symbolic answers.
vpa([sort(r), sort(s), sort(r2), sort(r3)], 8)
Those look to be in pretty close agreement.
