% FIGURE 5.9.1
%                             NTRI = NUMBER OF TRIANGLES
%                             M = NUMBER OF UNKNOWNS
%                             NSTEP = NUMBER OF TIME STEPS
      NTRI = 16;
      M = 5;
      NSTEP = 10;
%                             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)
               X = TRIX(KX,KTRI);
               Y = TRIY(KX,KTRI); 
%                             ASSIGN INITIAL VALUES, AS FUNCTION OF X AND Y
               A(K) = 1.0+X+Y;
               for IX=1:3
                  I = NODES(IX,KTRI);
                  if (I > 0)
                     L = max(L,abs(I-K));
                  end
               end
            end
         end
      end
%                             TAKE NSTEP STEPS OF SIZE DT
      DT = 1.0/NSTEP;
      for KT=1:NSTEP
      TKP1 = KT*DT;
%                             ZERO BAND MATRICES AMAT,BMAT AND VECTOR B.
      for K=1:M
         B(K) = 0.0;
         for ISUB=-L:L
            AMAT(K,L+1+ISUB) = 0.0;
            BMAT(K,L+1+ISUB) = 0.0;
         end
      end
%                             BEGIN ELEMENT-BY-ELEMENT ASSEMBLY OF
%                             AMAT,BMAT AND B.
      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,OMEGAT,OMEGAX,OMEGAY REPRESENT THE
%                             VALUE, T-DERIVATIVE, X-DERIVATIVE,
%                             AND Y-DERIVATIVE OF OMEGA AT THE TRIANGLE
%                             MIDPOINT.
         OMEGA  = 0.0;
         OMEGAT = 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);
               T = TKP1; 
%                             DEFINE BOUNDARY CONDITION COEFFICIENT R HERE
%                             AS FUNCTION OF X,Y AND T.  ALSO DEFINE 
%                             RT = dR/dt. 
               R = (1.0-T)*(1.0+X+Y);
               RT = -(1.0+X+Y);
               OMEGA  = OMEGA  + R*PHI (J);
               OMEGAT = OMEGAT + RT*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;
         T = TKP1;
%                             DEFINE PDE COEFFICIENT FUNCTIONS HERE
%                             AS FUNCTIONS OF X,Y AND T.  PDE IS:
%                             C(X,Y,T)*Ut = DIV(D(X,Y,T)*GRAD(U)) 
%                                         - A(X,Y,T)*U + F(X,Y,T)
%                             CMID,DMID,AMID,FMID ARE THE VALUES OF THE
%                             COEFFICIENT FUNCTIONS AT THE TRIANGLE MIDPOINT.
         CMID = 1.0;
         DMID = 1.0;
         AMID = 0.0;
         FMID = -(1.0+X+Y);
         AREA = abs(DET);
%                             CALCULATE CONTRIBUTIONS OF TRIANGLE KTRI
%                             TO MATRICES A AND B AND VECTOR B,
%                             USING THE MIDPOINT RULE TO APPROXIMATE
%                             ALL INTEGRALS.
         for KX=1:3
            K = NODES(KX,KTRI);
            if (K > 0)
            B(K) = B(K) + ...
                   (-CMID*OMEGAT*PHI(KX) ...
                    -DMID*(OMEGAX*PHIX(KX) + OMEGAY*PHIY(KX)) ...
                    -AMID*OMEGA*PHI(KX)  ...
                    +FMID*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;
                     BMAT(K,L+1+I-K) = BMAT(K,L+1+I-K) + ...
                           CMID*PHI(KX)*PHI(IX)*AREA;
                  end
               end
            end
         end
      end
%                             BACKWARD DIFFERENCE METHOD REQUIRES
%                             SOLUTION OF
%                              (BMAT + DT*AMAT)*AKP1 = DT*B + BMAT*AK
%                             EACH STEP.
      for K=1:M
         B(K) = DT*B(K);
         for I=max(K-L,1):min(K+L,M)
            AMAT(K,L+1+I-K) = BMAT(K,L+1+I-K) + DT*AMAT(K,L+1+I-K);
            B(K) = B(K) + BMAT(K,L+1+I-K)*A(I);
         end
      end
%                             CALL LBAND TO SOLVE BANDED LINEAR SYSTEM
      TKP1
      A = LBAND(AMAT,B,M,L)
      end