C FIGURE 3.4.2
      IMPLICIT DOUBLE PRECISION(A-H,O-Z)                                
C                              N = NUMBER OF POINTS IN X,Y DIRECTIONS   
C                              M = NUMBER OF TIME STEPS                 
      PARAMETER (N=40, M=200)                                           
      DIMENSION UKM1(0:N,0:N),UK(0:N,0:N),UKP1(0:N,0:N)                 
      DX = 1.D0/N                                                       
      DY = DX                                                           
      DT = 2.D0/M                                                       
C                              BOUNDARY AND INTERIOR VALUES SET TO ZERO 
C                              (BOUNDARY VALUES WILL NEVER CHANGE)      
      DO 5 I=0,N                                                        
      DO 5 J=0,N                                                        
         UKM1(I,J) = 0.0                                                
         UK  (I,J) = 0.0                                                
         UKP1(I,J) = 0.0                                                
    5 CONTINUE                                                          
C                              RESET VALUES OF UKM1,UK AT INTERIOR      
C                              POINTS TO INITIAL VALUES AT T0,T1        
      DO 10 I=1,N-2                                                     
      DO 10 J=1,N-I-1                                                   
         XI = I*DX                                                      
         YJ = J*DY                                                      
         UKM1(I,J) = XI*YJ*(XI+YJ-1)                                    
         UK  (I,J) = UKM1(I,J) + DT*DT*(XI+YJ)                          
   10 CONTINUE                                                          
C                              BEGIN MARCHING FORWARD IN TIME           
      DO 30 K=1,M                                                       
         DO 15 I=1,N-2                                                  
         DO 15 J=1,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/DY**2 )/ (1+DT/2)     
   15    CONTINUE                                                       
C                              UPDATE UKM1,UK                           
         DO 20 I=0,N                                                    
         DO 20 J=0,N                                                    
            UKM1(I,J) = UK(I,J)                                         
            UK  (I,J) = UKP1(I,J)                                       
   20    CONTINUE                                                       
         TKP1 = (K+1)*DT                                                
C                             OUTPUT SOLUTION AT X=0.25,Y=0.25          
         PRINT 25, TKP1,UKP1(N/4,N/4)                                   
   25    FORMAT (2E15.5)                                                
   30 CONTINUE                                                          
      STOP                                                              
      END