% FIGURE 3.4.2
%                              N = NUMBER OF POINTS IN X,Y DIRECTIONS
%                              M = NUMBER OF TIME STEPS
      N = 40;
      M = 200;
      DX = 1.0/N;
      DY = DX;
      DT = 2.0/M;
%                              BOUNDARY AND INTERIOR VALUES SET TO ZERO
%                              (BOUNDARY VALUES WILL NEVER CHANGE)
      for I=1:N+1
      for J=1:N+1
         UKM1(I,J) = 0.0;
         UK  (I,J) = 0.0;
         UKP1(I,J) = 0.0;
      end
      end
%                              RESET VALUES OF UKM1,UK AT INTERIOR
%                              POINTS TO INITIAL VALUES AT T0,T1
      for I=2:N-1
      for J=2:N-I+1
         XI = (I-1)*DX;
         YJ = (J-1)*DY;
         UKM1(I,J) = XI*YJ*(XI+YJ-1);
         UK  (I,J) = UKM1(I,J) + DT*DT*(XI+YJ);
      end
      end
%                              BEGIN MARCHING FORWARD IN TIME
      for K=1:M
         for I=2:N-1
         for J=2:N-I+1
            UKP1(I,J) = ( 2*UK(I,J) + (-1+DT/2)*UKM1(I,J)    ...
           +(UK(I+1,J)-2*UK(I,J)+UK(I-1,J))*DT^2/DX^2        ...
           +(UK(I,J+1)-2*UK(I,J)+UK(I,J-1))*DT^2/DX^2 )/ (1+DT/2);
         end
         end
%                              UPDATE UKM1,UK
         for I=1:N+1
         for J=1:N+1
            UKM1(I,J) = UK(I,J);
            UK  (I,J) = UKP1(I,J);
         end
         end
         TKP1 = (K+1)*DT
%                             OUTPUT SOLUTION AT X=0.25,Y=0.25
         U = UKP1(N/4+1,N/4+1)
      end