figure(2); clf; subplot(121);

%--- MY TRANSFER FUNCTION (below) ---
u1=(tf([130 1500],[1 10 100 100]))

%--- FIRST, JUST LET MATLAB CALCULATE GAIN AND PHASE MARGINS ---
margin(u1)
[g,p,gf,pf]=margin(u1)
subplot(222);

%--- NOW, PLOT THE ROOT LOCUS ---
rlocus(u1); title('Root Locus [3 poles, 1 zero]'); grid on;
axis equal
subplot(224);

%--- FIND WHERE ROOT LOCUS CROSSES IMAG AXIS (on margin of stability) ---
rlocus(u1,10.^[-2:.002:2]); title('Close-up of Root locus');
axis([-.1 .1 25 27]);
grid on

%--- USING ROOT LOCUS TO CALCULATE GAIN MARGIN (remainder below) ---
[z,p,n]=zpkdata(u1,'v')
gazc=26.2; % 'guess at the zero crossing' from close-up of root locus!

%--- CALCULATE CENTROID and plot ASYMPTOTES ---
sigma=(sum(p)-sum(z))/(length(p)-length(z))
subplot(222); hold on; 
plot([0 0]+sigma,[-40 40],'m--');
axis([-40 20 -30 30]);
plot([0 0],[gazc -gazc],'ks');

%--- ...THIS USES DISTANCES TO THE ZERO CROSSINGS! ---

dvalsp=sqrt((real(p)).^2 + (imag(p)-gazc).^2)
dvalsz=sqrt((real(z)).^2 + (imag(z)-gazc).^2)
k_crit=prod(dvalsp)/prod(dvalsz)
k_crit_test=n*g
fprintf('Matlab calculates a GAIN MARGIN of %.2f', g);
fprintf(' (%.2f dB)',20*log10(g));
fprintf('\nwhile root locus methods predict a gain margin of %.2f',k_crit/n);
fprintf(' (%.2f dB)\n',20*log10(k_crit/n));