% FIGURE 2.5.1
%                             N = NUMBER OF POINTS IN X,Y DIRECTIONS
%                             M = NUMBER OF TIME STEPS
      N = 20;
      M = 3200;
      DX = 1.0/N;
      DY = DX;
      DT = 0.5/M;
%                             BOUNDARY AND INTERIOR VALUES SET TO 1.
%                             (BOUNDARY VALUES WILL NOT CHANGE)
      for I=1:N+1
      for J=1:N+2
         UK(I,J) = 1.0;
         VK(I,J) = 1.0;
         UKP1(I,J) = 1.0;
         VKP1(I,J) = 1.0;
      end
      end
%                             BEGIN MARCHING FORWARD IN TIME
      for K=1:M
         for I=2:N
         for J=2:N+1
            UKP1(I,J) = UK(I,J) + DT*( VK(I,J)^2-3*UK(I,J)   ...
             + (UK(I+1,J)-2*UK(I,J)+UK(I-1,J))/DX^2          ...
             + (UK(I,J+1)-2*UK(I,J)+UK(I,J-1))/DY^2);
            VKP1(I,J) = VK(I,J) + DT*( -2*VK(I,J)^2 + 6*UK(I,J) ...
             + 4*(VK(I+1,J)-2*VK(I,J)+VK(I-1,J))/DX^2           ...
             + 4*(VK(I,J+1)-2*VK(I,J)+VK(I,J-1))/DY^2);
         end
         end
%                             HANDLE NORMAL DERIVATIVE BOUNDARY
%                             CONDITION AT Y=1.
         for I=2:N
            UKP1(I,N+2) = UKP1(I,N);
            VKP1(I,N+2) = VKP1(I,N);
         end
%                             UPDATE UK,VK
         for I=1:N+1
         for J=1:N+2
            UK(I,J) = UKP1(I,J);
            VK(I,J) = VKP1(I,J);
         end
         end
         TKP1 = K*DT
%                             PRINT SOLUTION AT (0.5,0.5)
         U = UKP1(N/2+1,N/2+1)
         V = VKP1(N/2+1,N/2+1)
      end