G=tf(1,[1,3,3,1]); N=10; Tdly=1; 
[nP,dP]=pade(Tdly,3); G_o=G*tf(nP,dP);
[k,L,T]=getfod(G_o); 
t1=0; t2=0; tend=50; U=1;
[Gc1,Kp,Ti,Td]=ziegler(3,[k,L,T,N]); [Kp,Ti,Td]

G_c=feedback(G_o*Gc1,1); 
step(G_c); hold on
num=G.num{1}; den=G.den{1}; 
[tt1,x,yy1]=sim('c7fsimd',tend,[]); plot(tt1,yy1)

[Kc,p,wc,w]=margin(G_o); Tc=2*pi/wc;
[Gc1,Kp,Ti,Td]=ziegler(3,[Kc,Tc,N]); [Kp,Ti,Td]

G_c=feedback(G_o*Gc1,1); step(G_c); hold on
[tt2,x,yy2]=sim('c7fsimd',tend,[]); plot(tt2,yy2)