% 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];
%                             INLINE FUNCTIONS DEFINING THE
%                             COEFFICIENTS
      FUNC = @(X,Y,T) 1.0;
      FUND = @(X,Y,T) 1.0;
      FUNA = @(X,Y,T) 0.0;
      FUNF = @(X,Y,T) -(1.0+X+Y);
      FUNR = @(X,Y,T) (1.0-T)*(1.0+X+Y);
      FUNRT= @(X,Y,T) -(1.0+X+Y);
      FUNH = @(X,Y) 1.0+X+Y;
%                             CALCULATE HALF-BANDWIDTH, L
      L = 0;
      for KTRI=1:NTRI
         for KX=1:3
            K = NODES(KX,KTRI);
            if (K > 0)
%                             ASSIGN INITIAL VALUES
               A(K) = FUNH(TRIX(KX,KTRI) , TRIY(KX,KTRI));
               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)
               R = FUNR(TRIX(J,KTRI),TRIY(J,KTRI),TKP1);
               RT=FUNRT(TRIX(J,KTRI),TRIY(J,KTRI),TKP1);
               OMEGA  = OMEGA  + R*PHI (J);
               OMEGAT = OMEGAT + RT*PHI(J);
               OMEGAX = OMEGAX + R*PHIX(J);
               OMEGAY = OMEGAY + R*PHIY(J);
            end
         end
%                             CMID,DMID,AMID,FMID ARE THE VALUES OF THE
%                             COEFFICIENT FUNCTIONS C(X,Y,T),D(X,Y,T),
%                             A(X,Y,T) AND F(X,Y,T) AT THE TRIANGLE
%                             MIDPOINT.
         XMID = (X1+X2+X3)/3.0;
         YMID = (Y1+Y2+Y3)/3.0;
         CMID = FUNC(XMID,YMID,TKP1);
         DMID = FUND(XMID,YMID,TKP1);
         AMID = FUNA(XMID,YMID,TKP1);
         FMID = FUNF(XMID,YMID,TKP1);
         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)*a(k+1) 
%                             = DT*b + BMAT*a(k) EACH TIME 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
