```function zs=verroots(a)
%    VERROOTS       Verified roots of a complex (or real) polynomial.
%
%    This is an INTLAB file. It requires to have INTLAB installed under
%    MATLAB to function properly.
%
%    For a complex (or real) vector a of length n+1 with a(1)=1,
%        zs=verroots(a)
%    returns a vector zs of verified ALL roots of the polynomial
%        a(1)*z^n+...+a(n)*z+a(n+1),   (1)
%    or yields no verified result.
%
%    COMMENT. It must be a(1)=1 (otherwise the computation breaks down), so
%    that the polynomial is monic.
%
%    Possible outcomes:
%
%    ~isnan(zs.inf(1,1)) :  zs is  a (generally complex) interval vector
%                           verified to contain the vector of all roots of
%                           the polynomial (1); it is sorted in ascending
%                           order of the real parts of the midpoints,
%     isnan(zs.inf(1,1)) :  no verified result (the interval vector zs
%                           consists of NaN's).
%

%
%    Polynomial roots computed as eigenvalues of the companion matrix.
%
%    This work was done during author's employment at the Anglo-American
%    University in Prague, Czech Republic.
%
%    WARRANTY
%
%    Because the program is licensed free of charge, there is
%    no warranty for the program, to the extent permitted by applicable
%    law. Except when otherwise stated in writing the copyright holder
%    and/or other parties provide the program "as is" without warranty
%    of any kind, either expressed or implied, including, but not
%    limited to, the implied warranties of merchantability and fitness
%    for a particular purpose. The entire risk as to the quality and
%    performance of the program is with you. Should the program prove
%    defective, you assume the cost of all necessary servicing, repair
%    or correction.
%
%    History
%
%    2008-05-15   first version
%    2008-05-28   version for posting
%
gr=getround;
setround(0);
% a: coefficients of the monic polynomial a(1)*z^n+...+a(n)*z+a(n+1), a(1)=1
a=a(:);
n=length(a)-1;  % degree of the polynomial
if n<1||a(1)~=1 % must be monic
zs=repmat(infsup(NaN,NaN),n,1);
setround(gr); return
end
% companion matrix
C=diag(ones(1,n-1),-1);
C(1,:)=-a(2:n+1); % first row
% verified eigenvalues of C
L=ol(C);    % diagonal matrix
zs=diag(L); % vector of verified zeros
% sorting in ascending order of real parts of the midpoints
[y,I]=sort(real(mid(zs)));
zs=zs(I);
setround(gr);
```