% FIGURE 0.4.4
      function X = LBAND(A,B,N,L)
%
%     ARGUMENT DESCRIPTIONS
%
%     A - (INPUT) A IS AN N BY 2*L+1 ARRAY CONTAINING THE BAND MATRIX.
%         A(I,L+1+J-I) CONTAINS THE MATRIX ELEMENT IN ROW I, COLUMN J.
%     X - (OUTPUT) X IS THE SOLUTION VECTOR OF LENGTH N.
%     B - (INPUT) B IS THE RIGHT HAND SIDE VECTOR OF LENGTH N.
%     N - (INPUT) N IS THE NUMBER OF EQUATIONS AND NUMBER OF UNKNOWNS
%         IN THE LINEAR SYSTEM.
%     L - (INPUT) L IS THE HALF-BANDWIDTH, DEFINED AS THE MAXIMUM VALUE
%         OF ABS(I-J) SUCH THAT AIJ IS NONZERO.
%
%                              COPY B TO X
      MD = L+1;
      for I=1:N
         X(I) = B(I);
         for J=-L:2*L
%                              COPY A TO A2, WHICH HAS L EXTRA
%                              DIAGONALS ABOVE, TO HOLD FILL-IN
            if (J <= L)
               A2(I,MD+J) = A(I,MD+J);
            else
               A2(I,MD+J) = 0.0;
            end
         end
      end
%                              BEGIN FORWARD ELIMINATION
      for K=1:N-1
         IBIG = K;
         BIG = abs(A2(K,MD));
%                              FIND THE LARGEST POTENTIAL PIVOT
         for I=K:min(K+L,N)
            if (abs(A2(I,MD+K-I)) > BIG)
               BIG = abs(A2(I,MD+K-I));
               IBIG = I;
            end
         end
         if (BIG == 0.0)
            error('The matrix is singular')
         end
%                              SWITCH ROW K WITH THE ROW (IBIG) 
%                              CONTAINING THE LARGEST POTENTIAL PIVOT
         for J=K:min(K+2*L,N)
            TEMP = A2(IBIG,MD+J-IBIG);
            A2(IBIG,MD+J-IBIG) = A2(K,MD+J-K);
            A2(K,MD+J-K) = TEMP;
         end
         TEMP = X(IBIG);
         X(IBIG) = X(K);
         X(K) = TEMP;
         for I=K+1:min(K+L,N)
            AMUL = -A2(I,MD+K-I)/A2(K,MD);
            if (AMUL ~= 0.0)
%                              ADD AMUL TIMES ROW K TO ROW I
               for J=K:min(K+2*L,N)
                  A2(I,MD+J-I) = A2(I,MD+J-I) + AMUL*A2(K,MD+J-K);
               end
               X(I) = X(I) + AMUL*X(K);
            end
         end
      end
      if (A2(N,MD) == 0.0)
         error('The matrix is singular')
      end
%                              BEGIN BACK SUBSTITUTION
      X(N) = X(N)/A2(N,MD);
      for K=N-1:-1:1
         for J=K+1:min(K+2*L,N)
            X(K) = X(K) - A2(K,MD+J-K)*X(J);
         end
         X(K) = X(K)/A2(K,MD);
      end