% FIGURE 0.4.3
      function X = LBAND0(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);
      end
%                              BEGIN FORWARD ELIMINATION
      for K=1:N-1
         if (A(K,MD) == 0.0)
            error ('Zero pivot encountered')
         end
         for I=K+1:min(K+L,N)
            AMUL = -A(I,MD+K-I)/A(K,MD);
            if (AMUL ~= 0.0)
%                              ADD AMUL TIMES ROW K TO ROW I
               for J=K:min(K+L,N)
                  A(I,MD+J-I) = A(I,MD+J-I) + AMUL*A(K,MD+J-K);
               end
               X(I) = X(I) + AMUL*X(K);
            end
         end
      end
      if (A(N,MD) == 0.0)
         error ('Zero pivot encountered')
      end
%                              BEGIN BACK SUBSTITUTION
      X(N) = X(N)/A(N,MD);
      for K=N-1:-1:1
         for J=K+1:min(K+L,N)
            X(K) = X(K) - A(K,MD+J-K)*X(J);
         end
         X(K) = X(K)/A(K,MD);
      end