a=[-2 -4];b=[0 0 -1 -8 -6];c=1818; Gs=zpk(a,b,c) T=0.1; Gz=c2d(Gs,T) HGz=c2d(Gs,T,'zoh') [z,p,k]=zpkdata(HGz) HGz1=zpk(z,p,k,T,'variable','z^-1') syms z a0 a1 a2 e0 Gcz=z^-1*(1+3.034*z^-1)*(a0+a1*z^-1+a2*z^-2) f1=subs(Gcz,z,1)-1 f2=subs(diff(Gcz,1),z,1) f3=subs(diff(Gcz,2),z,1) [a0j,a1j,a2j]=solve(f1,f2,f3) Gcz=subs(Gcz,[a0 a1 a2],[a0j a1j a2j]) Gez=(1-z^-1)^3*(1+e0*z^-1) f4=subs(Gez,z,-3.034)-1 e00=solve(f4) Gez=subs(Gez,e0,e00) Gz=(0.24551*z^-1)*(1+3.034*z^-1)*(1-0.8187*z^-1)*(1-0.6703*z^-1)*(1+0.2104*z^-1)/((1-z^-1)^2)/(1-0.9048*z^-1)/(1-0.5488*z^-1)/(1-0.4493*z^-1) Guz=Gcz/Gz Dyz=Gcz/Gez/Gz simplify(Gcz) [N,D]=numden(simplify(Gcz)) numc=sym2poly(N) denc=sym2poly(D) [N,D]=numden(simplify(Guz)) numu=sym2poly(N) denu=sym2poly(D) t=0:0.1:1 u=T*(t.^2)/2 hold on dlsim(numc,denc,u) dlsim(numu,denu,u) hold off [N,D]=numden(simplify(Dyz)) numdy=sym2poly(N) dendy=sym2poly(D)