s=tf('s')
wfal=8;
for n=2:10
e=sqrt(10^(wfal/10)-1);
u=(1/n)*asinh(1/e) ;
omc=1/(cosh(1/n*acosh(sqrt((2*e^2+1)/e^2))));
if mod(n,2)==0
    F=1;
    for k=1:(n/2)
    alp=((2*k-1)/(2*n))*pi;
    B=1/(cosh(u)^2-cos(alp)^2);
    A=2*B*sinh(u)*cos(alp);
    a=A/omc;
    b=B/(omc.^2);
    L=1;
    M=[b a 1]
    h=tf(L,M)
    F=F*h
    end
    H(n)=F
else
    F= 1/(1/sinh(u)/omc*s +1)
    for k=2:((n+1)/2)
        bet=((k-1)/n)*pi;
        B=1/(cosh(u)^2-cos(bet)^2)
        A=2*B*sinh(u)*cos(bet)
        a=A/omc
        b=B/(omc.^2)
        L=1;
        M=[b a 1]
        h=tf(L,M)
        F=F*h
    end
    H(n)=F
    end
end
%% BODE
figure(1)
bode(H(2),H(3),H(4),H(5),H(6),H(7),H(8),H(9),H(10))
legend('n=2','n=3','n=4','n=5','n=6','n=7','n=8','n=9','n=10')
%% JEDNOSTKOWY
figure(2)
step(H(2),H(3),H(4),H(5),H(6),H(7),H(8),H(9),H(10))
legend('n=2','n=3','n=4','n=5','n=6','n=7','n=8','n=9','n=10')
%% BIEGUNY
figure(3)
pzmap(H(2),H(4),H(6),H(8),H(10))
legend('n=2','n=4','n=6','n=8','n=10')
figure(4)
pzplot(H(3),H(5),H(7),H(9))
legend('n=3','n=5','n=7','n=9')
xlim([-1.2 0.1])
ylim([-1.1 1.1])