e9p3d.m

% Example 9.3 
phim = 50; w1=1.7; a0 = 1; Gp = tf([4],[1 3 2 0]);
Gpjw1 = evalfr(Gp,j*w1);
Gpjw1mag = abs(Gpjw1);
theta = -pi + phim/57.296 - angle(Gpjw1);
   a1 = (1 - a0*Gpjw1mag*cos(theta))/...
            (w1*Gpjw1mag*sin(theta))
   b1 = (cos(theta) - a0*Gpjw1mag)/(w1*sin(theta))
Gc = tf([a1/b1 a0/b1],[1 1/b1]);
T = minreal(Gc*Gp/(1+Gc*Gp));
pole(T), pause, margin(Gc*Gp)