I would use the bode function with three output arguments, then use findpeaks with the ‘mag’ and ‘wout’ arguments:
[mag,phase,wout] = bode(sys);
[pks,locs] = findpeaks(mag, wout);
Add the name-value pair arguments you need to get the result you want. The findpeaks function has considerable flexibility, but how much depends on your version of MATLAB, so be sure to read the relevant documentation for your version.