function c=polrange(Q,xrange,yrange,hx,hy)

% POLRANGE(Q,xrange,yrange,hx,hy) estimates the spectrum and  
% the boundary of the numerical range of a square monic matrix 
% polynomial of the form
%
%               m         m-1
%   P(z)  =  I z  +  A   z   + ... + A z  + A .
%                     m-1             1      0   
%
% (Note that the command tic-toc prints the execution time.)
%
%  INPUT:  Q = {A0,A1,...,Am} array of the matrix-coefficients,
%          xrange = [left bound on x-axis, right bound on x-axis],   
%          yrange = [lower bound on y-axis, upper bound on y-axis],    
%          hx = length of the grid for the x-axis (smaller than 0.5), 
%          hy = length of the grid for the y-axis (smaller than 0.5).
%
%  OUTPUT: c = contour plot points  
%          (boundary of the numerical range and eigenvalues of P(z)).
%
% Authors: Stavros Fatouros and Panayiotis Psarrakos,
%          Department of Mathematics,
%          National Technical University of Athens.
%
% Research supported by a grant of the EPEAEK project PYTHAGORAS II. 
% The project is co-funded by the European Social Fund (75%) and 
% Greek National Resources (25%).

tic

% Setting sizes and the degree of the matrix polynomial P(z)
[t,n]=size(Q);     % here n=m+1=degree+1  
degree=n-1;        % the degree of the matrix polynomial
[t,d]=size(Q{1});  % here d is the order of matrix-coefficients

% Constructing the grid
n1=round((xrange(2)-xrange(1))/hx);  
n2=round((yrange(2)-yrange(1))/hy);  

% Setting initial values on the grid
Z=zeros(n1+1,n2+1);  
Iterations=0;
expoints=0;

% Computing factorials
par=1;

for j=1:degree
  par=par*j;
  p(j)=par;
end  % for

% Main part of the program
for x=xrange(1):hx:xrange(2)
  I=round(((x-xrange(1))/hx)+1);

  for y=yrange(1):hy:yrange(2)
    J=round(((y-yrange(1))/hy)+1);
    mu=x+(y*i);    

    if Z(I,J)==0
      B0=Q{1};
    
      for j=1:degree
        B0=B0+(mu^j)*Q{j+1};
      end  % for
      
      r0=inrad(B0);
      Iterations=Iterations+1;
      
      if r0<0
        expoints=expoints+1;  
        f=maxderiv(Q,mu,p,degree);
        r=-r0/(-r0+f);  % usually smaller than 0.5 
        r=r/1.4;
        q1=round(r/hx);
        q2=round(r/hy);
      
        for k1=0:q1
        
          for k2=0:q2;
            Z(I+k1,J+k2)=r0;
            Y(J+k2)=yrange(1)+(J+k2-1)*hy;
          end  % for
       
          X(I+k1)=xrange(1)+(I+k1-1)*hx;
        end  % for
        
      else
        Z(I,J)=r0;
      end  % if

      Y(J)=yrange(1)+(J-1)*hy;
    end  % if 
 
    X(I)=xrange(1)+(I-1)*hx;
  end  % for
  
end  % for

toc

% Giving the number of iterations and plotting the figure
number_of_total_grid_points=(n1+1)*(n2+1)
number_of_iterations=Iterations
number_of_iterations_at_exterior_points=expoints
contour(X,Y,Z',[0 0],'b')
hold on
axis([xrange(1),xrange(2),yrange(1),yrange(2)])
xlabel('Real  Axis')
ylabel('Imaginary  Axis')

% Computing the companion matrix and the eigenvalues
C0=[];

for j=1:degree
  C0=[C0 -Q{j}];
end  % for

[l1,l2]=size(C0);
C=[zeros(l2-l1,l1) eye(l2-l1); inv(Q{n})*C0];
plot(real(eig(C)),imag(eig(C)),'r+')

hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f=maxderiv(Q,mu,p,degree);

% MAXDERIV computes for the 2-norm ||.||, the maximum of
% the quantities
%
%                (j)  
%      (1/j!) || P (mu) || ,   j = 1,2,...,m.
% 
% It requires also the degree of the matrix polynomial 
% P(z) and the vector p of the necessary factorials.

a=zeros(1,degree);
B={};

for j=1:degree
  B{j}=p(j)*Q{j+1};

  for k=1:degree-j
    B{j}=B{j}+(p(k+j)/p(k))*(mu^k)*Q{k+j+1};
  end  % for
    
  a(j)=norm(B{j},2)/p(j);
end  % for

f=max(a);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function r0=inrad(A)

% INRAD estimates the minimum, r0, among the maximum eigenvalues
% of all matrices  (1/2)*{e^(i*theta)*A + e^(-i*theta)*A}.

E=[];

for k=1:72
  theta=k*2*pi/72;
  A1=(cos(theta)+i*sin(theta))*A;
  H=0.5*(A1+A1');
  E=[E,max(eig(H))];
end  % for
    
r0=min(E);