e9p4.m

% Example 9.4 
phim = 40; w1=[1.7 2 3]; a0 = 1; Gp = tf([4],[1 3 2 0]);
for k=1:3
   Gpjw1 = evalfr(Gp,j*w1(k));
   Gpjw1mag = abs(Gpjw1);
   theta = -pi + phim/57.296 - angle(Gpjw1);
      a1 = (1 - a0*Gpjw1mag*cos(theta))/...
            (w1(k)*Gpjw1mag*sin(theta))
      b1 = (cos(theta) - a0*Gpjw1mag)/(w1(k)*sin(theta))
   Gc = tf([a1/b1 a0/b1],[1 1/b1]), pause
   T = minreal(Gc*Gp/(1+Gc*Gp));
   pole(T), pause, [Gm,Pm,Wcg,Wcp] = margin(Gc*Gp), pause
end