ap=1
as=30
Os=2.35
epsilon=sqrt(10^(ap/10)-1)
k1=sqrt(10^(ap/10)-1)/sqrt(10^(as/10)-1)
n=ceil(acosh(1/k1)/acosh(Os))
B=(sqrt(1+1/epsilon^2)+1/epsilon)^(1/n)
a=(B-1/B)/2
b=(B+1/B)/2
%format long
mu=(1:2:2*n-1)'
p=(-a*sin(mu*pi/(2*n))+i*b*cos(mu*pi/(2*n)))
Hjme=real(poly(p))           %jmenovatel
Hcit=1/epsilon/2^(n-1)       %èitatel
Phiz=i*cos(mu*pi/(2*n))        %nuly char. funkce
Phi=real(poly(Phiz))*2^(n-1) %char. funkce