clc
n = 25;
z1 = exp(2*pi*i*[1:n+1]/(n+1));
z2 = -W(z1/exp(1));
p = 1./factorial(n:-1:0).*((n+3).^[n:-1:0]);
C = compan(p);
z3 = eig(C);

plot(...
       real(z2), imag(z2), 'o', ...
       real(z3), imag(z3), '+');

axis(0.6*[-1 +1 -1 +1])