% FIGURE 5.6.1
%                             NTRI = NUMBER OF TRIANGLES
%                             M = NUMBER OF UNKNOWNS
      NTRI = 16;
      M = 5;
%                             ARRAYS TRIX,TRIY HOLD THE VERTICES OF THE
%                             TRIANGLES.  THE THREE VERTICES OF TRIANGLE
%                             KTRI ARE
%                                (TRIX(J,KTRI) , TRIY(J,KTRI))  (J=1,2,3)
%                             AND THE CORRESPONDING NODE NUMBERS ARE
%                                NODES(J,KTRI).
%                             THE NUMBER OF A VERTEX ON THE BOUNDARY IS 0
      TRIX = [0 0 1 1 0.5 0 0.5 1 1 0 0 0 0 0 -1 -1;
              0 1 1 0.5 0 0.5 1 1 0 0 0 0 0 -1 -1 0;
      0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 -0.5 -0.5 -0.5 -0.5 -0.5];
      TRIY = [0 -1 -1 0 0 0 0 0 1 1 0.5 0 0.5 1 1 0;
              -1 -1 0 0 0 0 0 1 1 0.5 0 0.5 1 1 0 0;
      -0.5 -0.5 -0.5 -0.5 -0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5 0.5];
      NODES = [0 0 0 0 2 0 2 0 0 0 4 0 4 0 0 0;
               0 0 0 2 0 2 0 0 0 4 0 4 0 0 0 0;
               1 1 1 1 1 3 3 3 3 3 3 5 5 5 5 5];
%                             CALCULATE HALF-BANDWIDTH, L
      L = 0;
      for KTRI=1:NTRI
         for KX=1:3
            K = NODES(KX,KTRI);
            if (K > 0)
               for IX=1:3
                  I = NODES(IX,KTRI);
                  if (I > 0)
                     L = max(L,abs(I-K));
                  end
               end
            end
         end
      end
%                             ZERO BANDED COEFFICIENT MATRIX AMAT AND
%                             RIGHT-HAND SIDE VECTOR B.
      for K=1:M
         B(K) = 0.0;
         for ISUB=-L:L
            AMAT(K,L+1+ISUB) = 0.0;
         end
      end
%                             BEGIN ELEMENT-BY-ELEMENT ASSEMBLY OF
%                             MATRIX AND RIGHT-HAND SIDE VECTOR.
      for KTRI=1:NTRI
         X1 = TRIX(1,KTRI);
         Y1 = TRIY(1,KTRI);
         X2 = TRIX(2,KTRI);
         Y2 = TRIY(2,KTRI);
         X3 = TRIX(3,KTRI);
         Y3 = TRIY(3,KTRI);
         DET = (X2-X1)*(Y3-Y1) - (X3-X1)*(Y2-Y1);
%                             PHI(J),PHIX(J),PHIY(J) REPRESENT THE
%                             VALUE, X-DERIVATIVE AND Y-DERIVATIVE
%                             OF PHI-J AT THE TRIANGLE MIDPOINT.
         PHI (1) = 1.0/3.0;
         PHIX(1) = (Y2-Y3)/DET;
         PHIY(1) = (X3-X2)/DET;
         PHI (2) = 1.0/3.0;
         PHIX(2) = (Y3-Y1)/DET;
         PHIY(2) = (X1-X3)/DET;
         PHI (3) = 1.0/3.0;
         PHIX(3) = (Y1-Y2)/DET;
         PHIY(3) = (X2-X1)/DET;
%                             OMEGA,OMEGAX,OMEGAY REPRESENT THE
%                             VALUE, X-DERIVATIVE AND Y-DERIVATIVE
%                             OF OMEGA AT THE TRIANGLE MIDPOINT.
         OMEGA  = 0.0;
         OMEGAX = 0.0;
         OMEGAY = 0.0;
         for J=1:3
            if (NODES(J,KTRI) <= 0)
               X = TRIX(J,KTRI);
               Y = TRIY(J,KTRI);
%                             DEFINE BOUNDARY CONDITION COEFFICIENT 
%                             R(X,Y) HERE.
               R = 0;
               OMEGA  = OMEGA  + R*PHI (J);
               OMEGAX = OMEGAX + R*PHIX(J);
               OMEGAY = OMEGAY + R*PHIY(J);
            end
         end
         X = (X1+X2+X3)/3.0;
         Y = (Y1+Y2+Y3)/3.0;
%                             DEFINE PDE COEFFICIENT FUNCTIONS HERE
%                             AS FUNCTIONS OF X AND Y.  PDE IS:
%                             DIV(D(X,Y)*GRAD(U)) - A(X,Y)*U + F(X,Y) = 0
%                             DMID,AMID,FMID ARE THE VALUES OF THE
%                             COEFFICIENT FUNCTIONS AT THE TRIANGLE MIDPOINT. 
         DMID = 1;
         AMID = 0;
         FMID = 1;
         AREA = abs(DET);
%                             CALCULATE CONTRIBUTIONS OF TRIANGLE KTRI
%                             TO MATRIX AND RIGHT-HAND SIDE INTEGRALS,
%                             USING THE MIDPOINT RULE TO APPROXIMATE
%                             ALL INTEGRALS.
         for KX=1:3
            K = NODES(KX,KTRI);
            if (K > 0)
               B(K) = B(K) +  ...
                   (FMID*PHI(KX)  ...
                   -DMID*(OMEGAX*PHIX(KX) + OMEGAY*PHIY(KX))  ...
                   -AMID*OMEGA*PHI(KX)) *AREA;
               for IX=1:3
                  I = NODES(IX,KTRI);
                  if (I > 0)
                     AMAT(K,L+1+I-K) = AMAT(K,L+1+I-K) +  ...
                           (DMID*(PHIX(KX)*PHIX(IX) + PHIY(KX)*PHIY(IX)) ...
                           +AMID*PHI(KX)*PHI(IX)) *AREA;
                  end
               end
            end
         end
      end
%                             CALL LBAND TO SOLVE BANDED LINEAR SYSTEM
      A = LBAND(AMAT,B,M,L)