% FIGURE 5.11.2
%                              BREAK POINTS ARE XPTS(K), K=1,N+1
%                              UNKNOWNS ARE A(I), I=1,M
      global N XPTS M A OMEGA LOCATE
      OMEGA = @OMEGA5112;
      LOCATE = @LOCATE5112;
%                              N = NUMBER OF SUBINTERVALS
      N = 4;
%                              M = NUMBER OF UNKNOWNS
      M = 2*N;
      ALPHA = 4.0;
      for K=1:N+1
         XPTS(K) = (K-1)/N;
      end
%                              INITIAL VALUES ASSIGNED
      for I=1:M
         [HORS,K] = LOCATE5112(I);
         if (HORS == 'H')
            A(I) = 1.0;
         else
            A(I) = 0.0;
         end
      end
      BETA(1) = 0.5 - 0.5/sqrt(3.0);
      BETA(2) = 0.5 + 0.5/sqrt(3.0);
%                              DO ITMAX INVERSE POWER ITERATIONS
      ITMAX = 15;
      for ITER=1:ITMAX
%                              SAVE PREVIOUS A VECTOR
         for I=1:M
            AOLD(I) = A(I);
         end
%                              CALCULATE MATRIX A-ALPHA*B, AND F=B*AOLD
         L = 2;
         for LL=0:N-1
         for J=1:2
            K = 2*LL+J;
            ZK = XPTS(LL+1) + BETA(J)*(XPTS(LL+2)-XPTS(LL+1));
            F(K) = 0.0;
            for I=max(K-L,1):min(K+L,M)
               AMAT(K,L+1+I-K) = -PHI(I,2,ZK) + 2*PHI(I,1,ZK) - PHI(I,0,ZK);
               BMAT = PHI(I,0,ZK);
               F(K) = F(K) + BMAT*A(I);
               AMAT(K,L+1+I-K) = AMAT(K,L+1+I-K) - ALPHA*BMAT;
            end
         end
         end
%                              SOLVE LINEAR SYSTEM USING BAND SOLVER
         A = LBAND(AMAT,F,M,L);
%                              ESTIMATE EIGENVALUE
         ANUM = 0.0;
         ADEN = 0.0;
         for I=1:M
            ANUM = ANUM + F(I)*AOLD(I);
            ADEN = ADEN + F(I)*A(I);
         end
         EIGEN = ANUM/ADEN + ALPHA;
%                              CALCULATE MAXIMUM ERROR
         ERMAX = 0.0;
         for J=0:100
            X = J/100.0;
            EXACTX = exp(X)*sin(2.0287579*X);
            EXACT1 = exp(1.0)*sin(2.0287579*1.0);
            ERR = abs(USOL(0,X)/USOL(0,1.0) - EXACTX/EXACT1);
            ERMAX = max(ERMAX,ERR);
         end
         ITER
         EIGEN
         ERMAX
      end
