WT = 0.4*s/(s + 1);

w = logspace(-2, 2, 200); % omega

WT_jw = frd(WT, w); % frequency response
P0_jw = frd(P0, w);

figure(3);
bodemag(WT_jw, 'r');
hold on;
for c=0.8:0.05:1.2
  Ptilde = 1/((s-1)*(c*s + 1));
  Ptilde_jw = frd(Ptilde, w);
  bodemag(Ptilde_jw/P0_jw - 1, 'b');
end
legh = legend('$W_T$', '$\Delta$', 'Location', 'SouthEast');
set(legh, 'interpreter', 'latex');

for phi=0:30:360
  Delta_jw = exp(1i*phi*pi/180);
  P_jw = P0_jw * (1 + WT_jw * Delta_jw);
  figure(1);
  bode(P_jw, 'g');
end
