% FIGURE 3.6.1
      function POLY = DEGENP(A,B,N,LAM0)
%
%  FUNCTION DEGENP CALLS HESSQZ TO REDUCE THE GENERALIZED EIGENVALUE PROBLEM
%
%                A*X = LAMBDA*B*X
%
%  TO A SIMILAR PROBLEM WITH THE SAME EIGENVALUES, WHERE A IS UPPER
%  HESSENBERG AND B IS UPPER TRIANGULAR, THEN CALCULATES THE COEFFICIENTS 
%  OF THE POLYNOMIAL DET(A-LAMBDA*B).  THE ROOTS OF THIS POLYNOMIAL WILL 
%  BE THE EIGENVALUES OF THE GENERALIZED PROBLEM. 
%
%  ARGUMENTS
%
%             ON INPUT                          ON OUTPUT
%             --------                          ---------
%
%    A      - THE N BY N A MATRIX.               
%
%    B      - THE N BY N B MATRIX.               
%
%    N      - THE SIZE OF MATRICES A AND B.
%
%    LAM0   - A SCALAR, SEE POLY. (SET 
%             LAM0=0 TO GET THE USUAL
%             POLYNOMIAL COEFFICIENTS)
%
%    POLY   -                                   VECTOR OF LENGTH N+1,
%                                               CONTAINING THE POLYNOMIAL
%                                               COEFFICIENTS. DET(A-LAMBDA*B)=
%                                                  SUM FROM I=0 TO N OF
%                                                  POLY(I+1)*(LAMBDA-LAM0)^I
%
%-----------------------------------------------------------------------
%
%                        CALL HESSQZ TO REDUCE A TO UPPER HESSENBERG AND
%                        B TO UPPER TRIANGULAR FORM
      [A B] = HESSQZ(A,B,N);
%                        DEGENP USES A RECURRENCE RELATION TO CALCULATE
%                        DET(A-LAMBDA*B) AND ALL N DERIVATIVES AT 
%                        LAMBDA = LAM0, FROM WHICH THE POLYNOMIAL
%                        COEFFICIENTS CAN BE FOUND. 
      DET(1,N+1) = 1.0;
      for I=2:N+1
         DET(I,N+1) = 0.0;
      end
      DET(1,N) = A(N,N)-LAM0*B(N,N);
      DET(2,N) = -B(N,N);
      for I=3:N+1
        DET(I,N) = 0.0;
      end
      for K=N-1:-1:1
         DET(1,K) = (A(K,K)-LAM0*B(K,K))*DET(1,K+1);
         for I=1:N
            DET(I+1,K) = (A(K,K)-LAM0*B(K,K))*DET(I+1,K+1)-B(K,K)*DET(I,K+1);
         end
         FACT = 1.0;
         for J=K+1:N
            FACT = -FACT*A(J,J-1);
            if (A(K,J)==0.0 & B(K,J)==0.0) 
               continue
            end
            DET(1,K) = DET(1,K) + FACT*(A(K,J)-LAM0*B(K,J))*DET(1,J+1);
            for I=1:N
               DET(I+1,K) = DET(I+1,K) + ...
               FACT*((A(K,J)-LAM0*B(K,J))*DET(I+1,J+1)-B(K,J)*DET(I,J+1));
            end
         end
      end
      for I=1:N+1
         POLY(I) = DET(I,1);
      end
