function c=sqprange(A,B,C)

% SQPRANGE(A,B,C) approximates the spectrum and the boundary 
% of the numerical range of the selfadjoint matrix polynomial  
% P(z) = A z^2 + B z + C  (i.e., A,B,C: Hermitian), where the  
% leading coefficient A is positive definite.
%
% (Note that the command tic-toc prints the execution time.)
%
% Authors: Peter Lancaster,
%          Department of Mathematics and Statistics,
%          University of Calgary. 
%
%          Panayiotis Psarrakos,
%          Department of Mathematics,
%          National Technical University of Athens.

tic

S=inv(sqrt(A));
B1=S'*B*S;  
C1=S'*C*S;
i=sqrt(-1);
K=C1+i*B1;
[n,m]=size(K);
r=2000;   
d=(2*pi)/r;   % step-size for partition   
z=[];   
d1=0.2;       % tolerance for linear interpolation
co=0;         % counter of number of linear segments
I=eye(n);
Z=zeros(n); 

for k=1:r
  AA=(cos(k*d)+i*sin(k*d))*K;
  HH=0.5*(AA+AA');   % Hermitian part of AA
  [X,L]=eig(HH);
  g=diag(real(L));   % vector of the eigenvalues of HH
  v=max(g)==g;
  x=X*v;
  metro=sqrt(x'*x);
  x=x/metro;
  z(k)=x'*C1*x+i*x'*B1*x;
end  % for

z(r+1)=z(1);

for  k=2:r+1
  zd=z(k)-z(k-1);
    
  if sqrt(zd'*zd)>d1,
          
    for l=1:100,
      zs(co*100+l)=z(k-1)+l*0.01*zd;
    end  % for 
          
    co = co +1 ;
  end  % if
  
end  % for

hold on

if co>0,
  plot(real(zs),imag(zs),'r.')
end  % if

hold  off

for  k=1:r,
  polnm=[1 imag(z(k)) real(z(k))] ;
  rt=roots(polnm);
  u(k)=rt(1) ;   
  w(k)=rt(2); 
end  % for

if co>0,

  for  k=1:co*100,
    polnm=[1 imag(zs(k)) real(zs(k))];  
    rt=roots(polnm);
    us(k)=rt(1);   
    ws(k)=rt(2); 
  end  % for

end  % if

C2=[Z I ; -C1 -B1];               
ee=eig(C2);  % the eigenvalues of the matrix polynomial

toc

plot(real(ee),imag(ee),'+')
hold on
plot(real(u),imag(u),'b.')
plot(real(w),imag(w),'c.')

if  co > 0,
  plot(real(us),imag(us),'r.')
  plot(real(ws),imag(ws),'r.')
end  % if

xlabel('Real  Axis')
ylabel('Imaginary  Axis')
hold off