<< Graphics`Graphics`
dB[z_] = 20 Log[10, Abs[z]] ; 
    dg[z_] = Im[Log[z]]/Degree ;

Q 6.4    p.143

default

go = 2/(s + 1)/(0.5  s + 1)/(0.2  s + 
            1)
goω = go /. s ->  I ω;

eigenvalues

chr = Numerator[Factor[1 + go]]

pole = Solve[chr == 0, s]

λ = s /. pole[[3]];
Print["ωn=", 
      Abs[λ], "  ζ=", -Re[λ]/
        Abs[λ]]

Bode plot of Go(s)

zu1 = 
      LogLinearPlot[{dB[goω], dg[goω] 40/180}, {ω, 0.1, 
          10}, PlotRange ->  {{0.1, 10}, {-40, 40}}, 
        GridLines ->  Automatic]

Nyquist plot of Go(s)

zu2 = 
      ParametricPlot[{Re[goω], Im[goω]}, {ω, 0.1, 10}, 
        PlotRange ->  {{-2, 2}, {-3, 1}}, 
        AspectRatio ->  1]

phase margin=45 Degree

gok = go k /. k -> 1.9
gokω = gok /. s ->  I ω;

pole = 
      Solve[Numerator[1 + gok] == 0, s]
λ = s /. pole[[3]];
Print["ωn=", 
      Abs[λ], "  ζ=", -Re[λ]/
        Abs[λ]]

zu145 = 
      LogLinearPlot[{dB[ gokω], dg[ gokω] 40/180}, {ω, 
          0.1, 10}, PlotRange ->  {{0.1, 10}, {-40, 40}}, 
        GridLines ->  Automatic]
zu245 = 
      ParametricPlot[{Re[gokω], Im[gokω]}, {ω, 0.1, 10}, 
        PlotRange ->  {{-2, 2}, {-3, 1}}, 
        AspectRatio ->  1]

stablility limit with phase margin=0

gok = go k /. k ->  6.3
gokω = gok  /. s ->  I ω;

pole = 
      Solve[Numerator[1 + gok] == 0, s]
λ = s /. pole[[3]];
Print["ωn=", 
      Abs[λ], "  ζ=", -Re[λ]/
        Abs[λ]]

zu10 = 
      LogLinearPlot[{dB[ gokω], dg[ gokω] 40/180}, {ω, 
          0.1, 10}, PlotRange ->  {{0.1, 10}, {-40, 40}}, 
        GridLines ->  Automatic]
zu20 = 
      ParametricPlot[{Re[gokω], Im[gokω]}, {ω, 0.1, 10}, 
        PlotRange ->  {{-2, 2}, {-3, 1}}, 
        AspectRatio ->  1]

Show[zu1, zu145, zu10]

Show[zu2, zu245, zu20]

Routh-Hurwitz criterion

chrk = Numerator[Factor[1 +  go k]]

h2 = {{0.4, k + 0.5}, {0.05, 0.85}};
MatrixForm[h2]

Factor[
      Det[h2]]

Time history simulation

gc = go k/(1 + go k)

ht1 = Re[InverseLaplaceTransform[gc /. k ->  1, s, t]]

Plot[ht1, {t, 0, 5}]

ht2 = 
        Re[InverseLaplaceTransform[gc /. k ->  1.9, s, 
            t]];
ht3 = 
        Re[InverseLaplaceTransform[gc /. k ->  6.3, s, 
            t]];

Plot[{ht1, ht2, ht3}, {t, 0, 5}, PlotRange ->  {-4, 4}, 
      GridLines ->  Automatic]