% Problem 9.18
%clear
%  System model
num=[0 0 0 40];
den=conv([1 4 0],[1 10]);
G=[];
P=[];
Wg=[];
Wp=[];
k=[];
for K=10:10:100,
    [Gm,Pm,Wcg,Wcp] = margin(K*num,den);
    G=[G 20*log10(Gm)];
    P=[P Pm];
    Wg=[Wg Wcg];
    Wp=[Wp Wcp];
    k=[k K];
end
i=1:10; i=i';
subplot(221)
plot(i,P)
grid
xlabel('(a) K*10')
ylabel('phase margin')
r=rlocus(num,den,k);
subplot(222)
plot(i,G)
grid
xlabel('(b) K*10')
ylabel('gain [dB]')
subplot(223)
plot(i,real(r(:,1)))
grid
xlabel('(c) K*10')
ylabel('Re(p1)')
subplot(224)
plot(i,real(r(:,2)))
grid
xlabel('(d) K*10')
ylabel('Re(p2), Re(p3)')
