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):
*parameters*
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]
% We'll think of detS as a polynomial in w
detS = simplify(det(S));
% Get vector of coefficients of the polynomial in w, and flip them
% so that the first element corresponds to the highest degree of w.
c = fliplr(coeffs(detS,w));
Mbds = [-1500 1500]; % bounds of M range
nM = 100; % number of values of M to get roots at
Mvec = linspace(Mbds(1),Mbds(2),nM)'; % vector of M values
rts = nan(nM,8); % create matrix to hold roots
for j = 1:nM % for each M gridpoint
% substitute in for M, then convert to numeric (double) vector
cj = double(subs(c,M,Mvec(j)));
rtj = roots(cj); % get all 8 roots of polynomial
rtj = rtj(imag(rtj)==0); % drop complex roots
rtj = sort(rtj); % sort remaining roots
nrtj = numel(rtj); % how many real roots there are
% Assign result to roots matrix. Note that if there were complex
% roots, some elements of rts(j,:) won't be assigned and will
% stay as NaN. So if you see an NaN in the end, you'll know that
% entry is not actually a root.
rts(j,1:nrtj) = rtj';
end
figure
hold on
plot(Mvec,rts)
xline(360)
ylim([0 2000])
xlim([-1500 1500])
