      SUBROUTINE INITAL
C...
C...  TWO DIMENSIONAL PARTIAL DIFFERENTIAL EQUATION WITH A MIXED PARTIAL
C...  DERIVATIVE
C...
C...  CONSIDER THE FOLLOWING TWO DIMENSIONAL, PDE WITH A MIXED SECOND
C...  ORDER, PARTIAL DERIVATIVE:
C...
C...     U  = U   - 3U   + U                                    (1)
C...      T    XX     XY    YY
C...
C...  A SOLUTION TO EQUATION (1) IS
C...
C...     U(X,Y,T) = EXP(X)*EXP(Y)*EXP(-T)                       (2)
C...
C...  THE FOLLOWING CODING IS FOR INITIAL AND BOUNDARY CONDITIONS
C...  WHICH ARE IDENTICALLY SATISFIED BY EQUATION (2).  THE INITIAL
C...  CONDITION IS
C...
C...     U(X,Y,0) = EXP(X)*EXP(Y)                               (3)
C...
C...  AND THE BOUNDARY CONDITIONS ARE
C...
C...     U(0,Y,T) = EXP(Y)*EXP(-T)                              (4)
C...
C...     U(1,Y,T) = EXP(1)*EXP(Y)*EXP(-T)                       (5)
C...
C...     U(X,0,T) = EXP(X)*EXP(-T)                              (6)
C...
C...     U(X,1,T) = EXP(1)*EXP(X)*EXP(-T)                       (7)
C...
C...  NOTE THAT THESE ARE DIRICHLET BOUNDARY CONDITIONS.  IN THE
C...  FOLLOWING CODING, THESE BOUNDARY CONDITIONS ARE IMPLEMENTED
C...  IN SUBROUTINE DERV BY PROGRAMMING THE CORRESPONDING TIME
C...  DERIVATIVES, E.G., BOUNDARY CONDITION (4) IS PROGRAMMED AS
C...
C...     U (O,Y,T) = -EXP(Y)*EXP(-T)                            (8)
C...      T
C...
C...  THE SOLUTION IS BY THE NUMERICAL METHOD OF LINES ON AN 11 X 11
C...  POINT GRID IN X AND Y.  THE MIXED PARTIAL DERIVATIVE, U  , IN
C...                                                         XY
C...  EQUATION (1) IS COMPUTED FROM THE PARTIAL DERIVATIVES U   AND
C...                                                         XX
C...  U  .  AN ALTERNATIVE APPROACH IS TO COMPUTE THIS MIXED PARTIAL
C...   YY
C...  FROM A FINITE DIFFERENCE APPROXIMATION, USING THE BOUNDARY CON-
C...  DITIONS AS REQUIRED.
C...
C...  COMMON AREA
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=11,NY=11)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(NX,NY)
     2      /F/ UT(NX,NY)
     3      /S/ UX(NX,NY), UY(NX,NY),UXX(NX,NY),UYY(NX,NY),UXY(NX,NY),
     4              X(NX),     Y(NY),        DX,        DY
     5      /I/        IP
C...
C...  GRID SPACING IN X AND Y
      DX=1.0D0/DFLOAT(NX-1)
      DY=1.0D0/DFLOAT(NY-1)
C...
C...  SPATIAL GRID AND INITIAL CONDITION
      DO 1 I=1,NX
      DO 1 J=1,NY
C...
C...     GRID IN X
         X(I)=DX*DFLOAT(I-1)
C...
C...     GRID IN Y
         Y(J)=DY*DFLOAT(J-1)
C...
C...     INITIAL CONDITION IS THE EXACT SOLUTION AT T = 0
         U(I,J)=EXACT(X(I),Y(J),T)
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=11,NY=11)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(NX,NY)
     2      /F/ UT(NX,NY)
     3      /S/ UX(NX,NY), UY(NX,NY),UXX(NX,NY),UYY(NX,NY),UXY(NX,NY),
     4              X(NX),     Y(NY),        DX,        DY
     5      /I/        IP
C...
C...  ARRAYS FOR THE INTERMEDIATE ONE DIMENSIONAL SOLUTION AND DERI-
C...  VATIVES
      DIMENSION US(NX), USX(NX), USY(NY), USXX(NX), USYY(NY)
C...
C...  BC AT X = 0
      DO 1 J=1,NY
         UT(1,J)=-EXACT(X(1),Y(J),T)
1     CONTINUE
C...
C...  BC AT Y = 0
      DO 2 I=1,NX
         UT(I,1)=-EXACT(X(I),Y(1),T)
2     CONTINUE
C...
C...  BC AT X = 1
      DO 3 J=1,NY
         UT(NX,J)=-EXACT(X(NX),Y(J),T)
3     CONTINUE
C...
C...  BC AT Y = 1
      DO 4 I=1,NX
         UT(I,NY)=-EXACT(X(I),Y(NY),T)
4     CONTINUE
C...
C...  UXX
C...
C...     0 LE Y LE 1
         DO 5 J=1,NY
C...
C...     0 LE X LE 1
         DO 6 I=1,NX
C...
C...     TRANSFER U(I,J) TO US(I)
6        US(I)=U(I,J)
C...
C...     USXX
         NL=1
         NU=1
         CALL DSS044(0.0D0,1.0D0,NX,US,USX,USXX,NL,NU)
C...
C...     TRANSFER USXX TO UXX
         DO 7 I=1,NX
7        UXX(I,J)=USXX(I)
5     CONTINUE
C...
C...  UYY
C...
C...     0 LE X LE 1
         DO 8 I=1,NX
C...
C...     0 LE Y LE 1
         DO 9 J=1,NY
C...
C...     TRANSFER U(I,J) TO US(J)
9        US(J)=U(I,J)
C...
C...     USYY
         CALL DSS044(0.0D0,1.0D0,NY,US,USY,USYY,NL,NU)
C...
C...     TRANSFER USYY TO UYY
         DO 10 J=1,NY
10       UYY(I,J)=USYY(J)
8        CONTINUE
C...
C...  UXY
      DO 11 I=2,NX-1
      DO 11 J=2,NY-1
         UXY(I,J)=
     1      (U(I+1,J+1)-2.0D0*U(I,J)+U(I-1,J-1))/(2.0D0*DX*DY)
     2      -UXX(I,J)*(DX/DY)/2.0D0-UYY(I,J)*(DY/DX)/2.0D0
11    CONTINUE
C...
C...  PDE
      DO 12 I=2,NX-1
      DO 12 J=2,NY-1
         UT(I,J)=UXX(I,J)-3.0D0*UXY(I,J)+UYY(I,J)
12    CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=11,NY=11)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/  U(NX,NY)
     2      /F/ UT(NX,NY)
     3      /S/ UX(NX,NY), UY(NX,NY),UXX(NX,NY),UYY(NX,NY),UXY(NX,NY),
     4              X(NX),     Y(NY),        DX,        DY
     5      /I/        IP
C...
C...  DIMENSION THE ARRAYS FOR PRINTING AND PLOTTING
      DIMENSION TP(51), UP(2,51), DIFF(51), UE(NX,NY)
C...
C...  MONITOR THE SOLUTION
      WRITE(*,*)IP,T
C...
C...  PRINT THE SOLUTION EVERY TENTH CALL TO SUBROUTINE PRINT
      IP=IP+1
      IF(((IP-1)/10*10).EQ.(IP-1))THEN
C...
C...  PRINT THE NUMERICAL U(X,Y,T)
      WRITE(NO,1)T,((U(I,J),I=1,NX,2),J=1,NY,2)
1     FORMAT(' T = ',F6.3,//,'   U(X,Y,T)',/,
     1 15X,'X = O',4X,'X = 0.2',3X,'X = 0.4',3X,'X = 0.6',
     2             3X,'X = 0.8',3X,'X = 1.0',/,
     33X,' Y = 0 ',6F10.3,/,3X,'Y = 0.2',6F10.3,/,3X,'Y = 0.4',6F10.3,/,
     43X,'Y = 0.6',6F10.3,/,3X,'Y = 0.8',6F10.3,/,3X,'Y = 1.0',6F10.3,/)
C...
C...  COMPUTE THE EXACT U(X,Y,T)
      DO 2 J=1,NY
      DO 2 I=1,NX
         UE(I,J)=EXACT(X(I),Y(J),T)
2     CONTINUE
C...
C...  PRINT THE EXACT U(X,Y,T)
      WRITE(NO,3)((UE(I,J),I=1,NX,2),J=1,NY,2)
3     FORMAT(//,'  UE(X,Y,T)',/,
     1 15X,'X = 0',4X,'X = 0.2',3X,'X = 0.4',3X,'X = 0.6',
     2             3X,'X = 0.8',3X,'X = 1.0',/,
     33X,' Y = 0 ',6F10.3,/,3X,'Y = 0.2',6F10.3,/,3X,'Y = 0.4',6F10.3,/,
     43X,'Y = 0.6',6F10.3,/,3X,'Y = 0.8',6F10.3,/,3X,'Y = 1.0',6F10.3,/)
C...
C...  PRINT THE NUMERICAL UT(X,Y,T)
      WRITE(NO,4)((UT(I,J),I=1,NX,2),J=1,NY,2)
4     FORMAT(//,'  UT(X,Y,T)',/,
     1 15X,'X = 0',4X,'X = 0.2',3X,'X = 0.4',3X,'X = 0.6',
     2             3X,'X = 0.8',3X,'X = 1.0',/,
     33X,' Y = 0 ',6F10.3,/,3X,'Y = 0.2',6F10.3,/,3X,'Y = 0.4',6F10.3,/,
     43X,'Y = 0.6',6F10.3,/,3X,'Y = 0.8',6F10.3,/,3X,'Y = 1.0',6F10.3,/)
C...
C...  PRINT THE NUMERICAL UXX(X,Y,T)
      WRITE(NO,5)((UXX(I,J),I=1,NX,2),J=1,NY,2)
5     FORMAT(//, ' UXX(X,Y,T)',/,
     1 15X,'X = 0',4X,'X = 0.2',3X,'X = 0.4',3X,'X = 0.6',
     2             3X,'X = 0.8',3X,'X = 1.0',/,
     33X,' Y = 0 ',6F10.3,/,3X,'Y = 0.2',6F10.3,/,3X,'Y = 0.4',6F10.3,/,
     43X,'Y = 0.6',6F10.3,/,3X,'Y = 0.8',6F10.3,/,3X,'Y = 1.0',6F10.3,/)
C...
C...  PRINT THE NUMERICAL UXY IN THE INTERIOR AND THE ANALYTICAL UXY ON
C...  THE BOUNDARIES.  THE LATTER ARE FIRST COMPUTED IN THE FOLLOWING
C...  SERIES OF FOUR DO LOOPS
      DO 6 J=1,NY
6        UXY(1,J) =EXACT(X(1) ,Y(J),T)
      DO 7 J=1,NY
7        UXY(NX,J)=EXACT(X(NX),Y(J),T)
      DO 8 I=1,NX
8        UXY(I,1) =EXACT(X(I) ,Y(1),T)
      DO 9 I=1,NX
9        UXY(I,NY)=EXACT(X(I),Y(NY),T)
      WRITE(NO,10)((UXY(I,J),I=1,NX,2),J=1,NY,2)
10    FORMAT(//, ' UXY(X,Y,T)',/,
     1 15X,'X = 0',4X,'X = 0.2',3X,'X = 0.4',3X,'X = 0.6',
     2             3X,'X = 0.8',3X,'X = 1.0',/,
     33X,' Y = 0 ',6F10.3,/,3X,'Y = 0.2',6F10.3,/,3X,'Y = 0.4',6F10.3,/,
     43X,'Y = 0.6',6F10.3,/,3X,'Y = 0.8',6F10.3,/,3X,'Y = 1.0',6F10.3,/)
C...
C...  PRINT THE NUMERICAL UYY
      WRITE(NO,11)((UYY(I,J),I=1,NX,2),J=1,NY,2)
11    FORMAT(//, ' UYY(X,Y,T)',/,
     1 15X,'X = 0',4X,'X = 0.2',3X,'X = 0.4',3X,'X = 0.6',
     2             3X,'X = 0.8',3X,'X = 1.0',/,
     33X,' Y = 0 ',6F10.3,/,3X,'Y = 0.2',6F10.3,/,3X,'Y = 0.4',6F10.3,/,
     43X,'Y = 0.6',6F10.3,/,3X,'Y = 0.8',6F10.3,/,3X,'Y = 1.0',6F10.3,/)
      END IF
C...
C...  STORE THE NUMERICAL AND ANALYTICAL U(1/2,1/2,T) AND THEIR
C...  DIFFERENCE FOR PRINTING
      TP(IP)=T
      UP(1,IP)=U(6,6)
      UP(2,IP)=EXACT(X(6),Y(6),T)
      DIFF(IP)=(UP(1,IP)-UP(2,IP))/UP(2,IP)*100.0D0
C...
C...  PRINT THE NUMERICAL AND ANALYTICAL U(1/2,1/2,T) AND THEIR
C...  DIFFERENCE VS T
      IF(IP.NE.51)RETURN
      WRITE(NO,12)(TP(I),UP(1,I),UP(2,I),DIFF(I),I=1,51)
12    FORMAT(//,10X,
     1 '   T  ',6X,'U(1/2,1/2,T)','   U(1/2,1/2,T)','    O/O DIFF  ',/,
     2         21X,'     NUM     ','       ANAL    '                ,/,
     3 (10X,F6.2,2F15.4,F12.2))
      RETURN
      END
      DOUBLE PRECISION FUNCTION EXACT(X,Y,T)
C...
C...  FUNCTION EXACT COMPUTES THE EXACT SOLUTION
C...
C...     U(X,Y,T)=EXP(X)*EXP(Y)*EXP(-T)
C...
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      EXACT=DEXP(X)*DEXP(Y)*DEXP(-T)
      RETURN
      END
