// 3ŽŸŒn (ƒÖ^3)/(s^3+ƒżƒÖs+ƒŔƒÖ^2s+ƒÖ^3) ‚Ěƒ{[ƒhü}
//       = (wn/m^2)*(m*wn)^2/( (s+wn/m^2)*(s^2+2Ā(mwn)s+(mwn)^2) )

s=poly(0,'s');

zeta = 0.5;
m = (2*zeta)^(1/3);

aa = 2*zeta*m+1/m^2;
bb = m^2+2*zeta/m;

//aa = 3.15;
//bb = 3.18;

fn = 1;               // [Hz]
wn = 2*%pi*fn;         // [rad/sec]

h=syslin('c',wn^3/(s^3+aa*wn*s^2+bb*wn^2*s+wn^3));

bode(h);