% Matt Stone under the guidance of Mladen Victor Wickerhauser

% This program plots the zeros of the
% normalized reversed bessel polynomial theta_n(nz)
% note: zeros of these polynomial are poles of the bessel (low-pass) filter

maxN = 40;

w = 1;

coeffBessel = zeros(maxN+1);

for n = 1:maxN
    t = zeros(maxN+1, 1);
    % normalization puts n^n
    t(1) = n^n;
    for i = 1:n
        t(i+1) = t(i) * (((n+i)*(n-i+1))/(n*2*i*w));
    end
     coeffBessel(:, n) = t;
end

% With this, we plot the roots

% uncommenting the hold on and hold off command
%   plots all roots of consecutive normalized reversed Bessel polynomials 
%   on one graph

% hold on;
for n = 3:maxN
    rBessel = roots(coeffBessel(1:n, n));
    plot(rBessel, 'x');
    pause(0.5);
end
% hold off;