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