C FIGURE 2.5.1
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C                             N = NUMBER OF POINTS IN X,Y DIRECTIONS    
C                             M = NUMBER OF TIME STEPS                  
      PARAMETER (N=20, M=3200)                                          
      DIMENSION UK(0:N,0:N+1),UKP1(0:N,0:N+1)                           
      DIMENSION VK(0:N,0:N+1),VKP1(0:N,0:N+1)                           
      DX = 1.D0/N                                                       
      DY = DX                                                           
      DT = 0.5D0/M                                                      
C                             BOUNDARY AND INTERIOR VALUES SET TO 1.    
C                             (BOUNDARY VALUES WILL NOT CHANGE)         
      DO 5 I=0,N                                                        
      DO 5 J=0,N+1                                                      
         UK(I,J) = 1.0                                                  
         VK(I,J) = 1.0                                                  
         UKP1(I,J) = 1.0                                                
         VKP1(I,J) = 1.0                                                
    5 CONTINUE                                                          
C                             BEGIN MARCHING FORWARD IN TIME            
      DO 30 K=1,M                                                       
         DO 10 I=1,N-1                                                  
         DO 10 J=1,N                                                    
            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)                 
   10    CONTINUE                                                       
C                             HANDLE NORMAL DERIVATIVE BOUNDARY         
C                             CONDITION AT Y=1.                         
         DO 15 I=1,N-1                                                  
            UKP1(I,N+1) = UKP1(I,N-1)                                   
            VKP1(I,N+1) = VKP1(I,N-1)                                   
   15    CONTINUE                                                       
C                             UPDATE UK,VK                              
         DO 20 I=0,N                                                    
         DO 20 J=0,N+1                                                  
            UK(I,J) = UKP1(I,J)                                         
            VK(I,J) = VKP1(I,J)                                         
   20    CONTINUE                                                       
         TKP1 = K*DT                                                    
C                             PRINT SOLUTION AT (0.5,0.5)               
         PRINT 25, TKP1,UKP1(N/2,N/2),VKP1(N/2,N/2)                     
   25    FORMAT (3E15.5)                                                
   30 CONTINUE                                                          
      STOP                                                              
      END