I may be misunderstanding things, but it looks to me like the equation you're trying (det(S)=0) to solve is an 8th-order polynomial in w. Algebraic solutions to 8th-order polynomials don't exist except in very special cases. When you apply solve to a case like this, the "solution" isn't typically actually an explicit algebraic solution. Instead, it's an instruction on how to get an explicit numeric solution if necessary. Specifically, the solutions found by solve will look something like:
speeds(1) = root(<polynomial in z>,z,1)
speeds(2) = root(<polynomial in z>,z,2)
speeds(8) = root(<polynomial in z>,z,8)
where root(<polynomial in z>,z,j) means the j-th root of the specified polynomial in z. If and when MATLAB needs to find an actual explicit solution numerically, it finds these roots numerically.
In your case, I'm guessing the cofficients of the above polynomial in z also contain one (and only one) generic symbolic variable, with everything else numeric. For the sake of the explanation, let's suppose M is the generic symbolic. When you call fplot(speeds(j)), MATLAB first creates a numeric grid for M, and then for each point in that grid, calls a numeric root-finding algorithm to get the desired root.
This is quite inefficient, for (at least) two reasons. First, MATLAB uses variable-precision arithmetic to get numeric solutions to symbolic equations, and that can be quite slow. Unless you need a very high degree of accuracy (e.g., solutions determined to 32 significant digits), you can likely achieve speed improvements by solving using regular double-precision arithmetic.
Second, there's some duplication of effort here. For example, consider the first gridpoint of M, call it m1. When you do fplot(speeds(1)), MATLAB first sets M=m1 in the polynomial and then finds the 1st root. When you do fplot(speeds(2)), MATLAB does the same thing except finds the 2nd root, fplot(speeds(3)) the 3rd root, etc. This means you are finding each of the roots of the polynomial with M=m1 separately. Since nearly all of the work involved in finding all the roots of the polynomial with M=m1 will already have been done when you call fplot(speeds(1)), there are wasted calculations when you then do fplot(speeds(2)), and then fplot(speeds(3)), etc.
Instead of how you've done it, I would do something like the following (again, for the sake of argument assuming M is the only non-numeric variable in S other than w):
S = [k(1,1)-M*w^2 k(1,2) 0 0;
k(1,2) k(2,2)-Idgear*w^2 0 i*Ipgear*w^2*Omega;
0 0 k(1,1)-M*w^2 k(1,2);
0 -i*Ipgear*w*Omega k(1,2) k(2,2)-Idgear*w^2]
detS = simplify(det(S));
c = fliplr(coeffs(detS,w));
Mbds = [-1500 1500];
nM = 100;
Mvec = linspace(Mbds(1),Mbds(2),nM)';
rts = nan(nM,8);
for j = 1:nM
cj = double(subs(c,M,Mvec(j)));
rtj = roots(cj);
rtj = rtj(imag(rtj)==0);
rtj = sort(rtj);
nrtj = numel(rtj);
rts(j,1:nrtj) = rtj';