      SUBROUTINE INITAL
C...
C...  ONE-DIMENSIONAL, NONLINEAR, PARABOLIC PDE BY SECOND AND FOURTH-
C...  ORDER FORMULAS FOR SECOND-ORDER SPATIAL DERIVATIVES
C...
C...  CONSIDER THE FOLLOWING NONLINEAR PARABOLIC PARTIAL DIFFERENTIAL
C...  EQUATION (PDE) (1)
C...
C...     U*U  = (1/X**2)*((X**2)*U*U )  + 5*(U**2) + 4*X*U*U    (1)
C...        T                       X X                     X
C...
C...     U(X,0) = EXP(1 - X**2)                                 (2)
C...
C...     U (0,T) = 0, U(1,T) = EXP(-T)                       (3)(4)
C...      X
C...
C...  THE ANALYTIC SOLUTION TO EQUATIONS (1) TO (4) IS
C...
C...     U(X,T) = EXP(1 - X**2 - T)                             (5)
C...
C...  (1)  THIS PROBLEM WAS PROVIDED BY ROBERT D. SKEEL, DEPARTMENT OF
C...       COMPUTER SCIENCE, UNIVERSITY OF ILLINOIS, AND MARTIN BERZINS,
C...       DEPARTMENT OF COMPUTER STUDIES, UNIVERSITY OF LEEDS
C...
C...  IN THE FOLLOWING CODING, EQUATIONS (1) TO (4) ARE SOLVED BY THE
C...  NUMERICAL METHOD OF LINES USING THE SECOND-DERIVATIVE FORMULAS
C...  REPORTED BY SCHIESSER AND STENGLE (2).  THE T INTEGRATION IS
C...  PERFORMED BY DASSL (3).  BOTH SECOND AND FOURTH-ORDER FORMULAS
C...  ARE USED SO THAT THE EFFECT OF THE ORDER OF THE APPROXIMATION IN
C...  COMPUTING SECOND-ORDER DERIVATIVES CAN BE OBSERVED.
C...
C...     (2)  SCHIESSER, W. E., AND G. A. STENGLE, SECOND DERIVATIVE
C...          FORMULAS FOR USE IN THE NUMERICAL METHOD OF LINES,
C...          PROCEEDINGS OF THE IMACS SIXTH INTERNATIONAL SYMPOSIUM
C...          ON COMPUTER METHODS FOR PARTIAL DIFFERENTIAL EQUATIONS,
C...          LEHIGH UNIVERSITY, BETHLEHEM, PA, JUNE 23-26, 1987.
C...
C...     (3)  DASSL IS A TRANSPORTABLE FORTRAN CODE FOR THE NUMERICAL
C...          SOLUTION OF DIFFERENTIAL-ALGEBRAIC SYSTEMS DEVELOPED BY
C...          LINDA R. PETZOLD OF THE LAWRENCE LIVERMORE LABORATORY
C...
C...  BEFORE PROCEEDING WITH THE CODING FOR EQUATIONS (1) TO (4), A
C...  DETAIL SHOULD BE CONSIDERED.  THE FIRST TERM ON THE RIGHT HAND
C...  SIDE OF EQUUATON (1) CAN BE EXPANDED BY THE CHAIN RULE FOR
C...  DIFFERENTIATION OF A PRODUCT
C...
C...     (1/X**2)*((X**2)*U*U )  = U*U   + (U )**2 + (2/X)*U*U
C...                         X X      XX     X                X
C...
C...  ADDITIONALLY, THE SECOND TERM IS INDETERMINATE FOR X = 0 AS A
C...  CONSEQUENCE OF BOUNDARY CONDITION (3), I.E.,
C...
C...       LIM   (2/X)*U*U  = 2*(U*U   + (U )**2)               (6)
C...      X--+0           X         XX     X
C...
C...  THIS RESULT FOR X = 0 IS USED IN THE SUBSEQUENT METHOD OF LINES
C...  CODING.
C...
C...  THE CODING FOR THE SECOND-ORDER FORMULAS IS CONTAINED IN SUBROU-
C...  TINES INTAL1, DERV1 AND PRINT1 WHICH ARE CALLED BY SETTING NORD=2
C...  IN THE FOLLOWING CODE.  THE CODING FOR THE FOURTH-ORDER FORMULAS
C...  IS IN SUBROUTINES INTAL2, DERV2 AND PRINT2 WHICH ARE CALLED BY
C...  SETTING NORD=4.
C...
C...  COMMON AREA FOR INTERFACE WITH SUBROUTINES INITAL, DERV AND
C...  PRINT
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=161)
      COMMON/T/        T,     NSTOP,     NORUN
     1      /Y/    U(NX)
     2      /F/   UT(NX)
     3      /S/   UX(NX),   UXX(NX)
     4      /C/    X(NX),         N,      NORD,         IP
C...
C...  SELECT THE ORDER OF THE SECOND-DERIVATIVE FORMULAS
      NORD=2
      NORD=4
C...
C...  INITIAL CONDITION
      IF(NORD.EQ.2)CALL INTAL1
      IF(NORD.EQ.4)CALL INTAL2
      RETURN
      END
      SUBROUTINE DERV
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=161)
      COMMON/T/        T,     NSTOP,     NORUN
     1      /Y/    U(NX)
     2      /F/   UT(NX)
     3      /S/   UX(NX),   UXX(NX)
     4      /C/    X(NX),         N,      NORD,         IP
C...
C...  SELECT THE SECOND-DERIVATIVE FORMULAS
      IF(NORD.EQ.2)CALL DERV1
      IF(NORD.EQ.4)CALL DERV2
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=161)
      COMMON/T/        T,     NSTOP,     NORUN
     1      /Y/    U(NX)
     2      /F/   UT(NX)
     3      /S/   UX(NX),   UXX(NX)
     4      /C/    X(NX),         N,      NORD,         IP
C...
C...  SELECT THE OUTPUT ROUTINE
      IF(NORD.EQ.2)CALL PRINT1(NI,NO)
      IF(NORD.EQ.4)CALL PRINT2(NI,NO)
      RETURN
      END
      SUBROUTINE INTAL1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=161)
      COMMON/T/        T,     NSTOP,     NORUN
     1      /Y/    U(NX)
     2      /F/   UT(NX)
     3      /S/   UX(NX),   UXX(NX)
     4      /C/    X(NX),         N,      NORD,         IP
C...
C...  SET THE NUMBER OF GRID POINTS
      IF(NORUN.EQ.1)N=11
      IF(NORUN.EQ.2)N=21
      IF(NORUN.EQ.3)N=41
      IF(NORUN.EQ.4)N=81
      IF(NORUN.EQ.5)N=161
C...
C...  INITIAL CONDITION
      DO 1 I=1,N
         X(I)=DFLOAT(I-1)/DFLOAT(N-1)
         U(I)=DEXP(1.0D0-X(I)**2)
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV1
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=161)
      COMMON/T/        T,     NSTOP,     NORUN
     1      /Y/    U(NX)
     2      /F/   UT(NX)
     3      /S/   UX(NX),   UXX(NX)
     4      /C/    X(NX),         N,      NORD,         IP
C...
C...  FIRST DERIVATIVE UX INCLUDING THE DIRCHLET BOUNDARY CONDITION
C...  AT X = 1
      U (N)= DEXP(-T)
      UT(N)=-DEXP(-T)
      CALL DSS002(0.D0,1.D0,N,U,UX)
C...
C...  SECOND DERIVATIVE INCLUDING THE SYMMETRY (NEUMANN) BOUNDARY
C...  CONDITION AT X =0
      NL=2
      NU=1
      UX(1)=0.D0
      CALL DSS042(0.D0,1.D0,N,U,UX,UXX,NL,NU)
C...
C...  PDE
      DO 1 I=1,N-1
C...     X = 0
         IF(I.EQ.1)THEN
            UT(1)=(1.D0/U(1))*3.D0*(U(1)*UXX(1)+UX(1)**2)+5.D0*U(1)
C...
C...     0 LT X LT 1
         ELSE
            UT(I)=(1.D0/U(I))*(U(I)*UXX(I)+UX(I)**2+2.D0/X(I)*U(I)*UX(I)
     1           +5.D0*U(I)**2+4.D0*X(I)*U(I)*UX(I))
         END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT1(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=161)
      COMMON/T/        T,     NSTOP,     NORUN
     1      /Y/    U(NX)
     2      /F/   UT(NX)
     3      /S/   UX(NX),   UXX(NX)
     4      /C/    X(NX),         N,      NORD,         IP
C...
C...  ARRAYS FOR PLOTTING
      DIMENSION SLOGN(11), SLOGE(11)
C...
C...  MONITOR THE CALCULATION
      WRITE(*,*)NORD,NORUN,N
C...
C...  COMPUTE THE EXACT SOLUTION, THE DIFFERENCE BETWEEN THE NUMERICAL
C...  AND ANALYTICAL SOLUTIONS, AND SAVE THE LARGEST NUMERICAL ERROR
      ERRMAX=0.D0
      DO 1 I=1,N-1
         UA=DEXP(1.0D0-X(I)**2-T)
         ERR=DABS(UA-U(I))
         IF(ERRMAX.LT.ERR)ERRMAX=ERR
1     CONTINUE
C...
C...  LABEL AND PRINT THE LARGEST ERROR (NOTE - AT T = 0, THE ERROR
C...  IS ZERO)
      IP=IP+1
      IF(IP.EQ.1)THEN
C...
C...     PRINT A HEADING AT T = 0
         WRITE(NO,3)
3        FORMAT(//,8X,' SOLUTION BY SECOND ORDER APPROXIMATIONS',//,
     1   19X,'T',4X,'N',9X,'ERRMAX')
         ERRPLT=ERRMAX
      ELSE IF(IP.NE.1)THEN
C...
C...     PRINT ERRMAX VS T
         WRITE(NO,2)T,N,ERRMAX
2        FORMAT(10X,F10.3,I5,D15.4)
         IF(ERRMAX.GT.ERRPLT)ERRPLT=ERRMAX
      END IF
C...
C...  STORE LOG10(ERRPLT) VS LOG10(N) FOR SUBSEQUENT PLOTTING
      IF(IP.EQ.11)THEN
         SLOGN(NORUN)=DLOG10(DFLOAT(N))
         SLOGE(NORUN)=DLOG10(ERRPLT)
C...
C...  AT THE END OF THE FIFTH RUN (WITH N = 161), PLOT LOG10(ERR)
C...  VS LOG10(N)
      IF(NORUN.EQ.5)THEN
         WRITE(NO,4)(SLOGN(I),SLOGE(I),I=1,NORUN)
4        FORMAT(//,19X,' DATA FOR PLOT',//,
     1          15X,'LOG10(N)',6X,'LOG10(ERR)',/,
     2          (7X,2F15.4))
         CALL SPLOTS(1,NORUN,SLOGN,SLOGE)
C...
C...     COMPUTE AND PRINT ESTIMATED ORDER
         SLOPE=(SLOGE(1)-SLOGE(NORUN-1))/(SLOGN(1)-SLOGN(NORUN-1))
         WRITE(NO,5)SLOPE
5        FORMAT(//,' LOG10(MAXERR) VS LOG10(N), APPARENT ORDER = ',
     1             F6.2)
      END IF
      END IF
      RETURN
      END
      SUBROUTINE INTAL2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=161)
      COMMON/T/         T,     NSTOP,     NORUN
     1      /Y/    U(NX)
     2      /F/   UT(NX)
     3      /S/   UX(NX),  UXX(NX)
     4      /C/    X(NX),         N,      NORD,         IP
C...
C...  SET THE NUMBER OF GRID POINTS
      IF(NORUN.EQ.1)N=11
      IF(NORUN.EQ.2)N=21
      IF(NORUN.EQ.3)N=41
      IF(NORUN.EQ.4)N=81
      IF(NORUN.EQ.5)N=161
C...
C...  INITIAL CONDITION
      DO 1 I=1,N
      X(I)=DFLOAT(I-1)/DFLOAT(N-1)
      U(I)=DEXP(1.0D0-X(I)**2)
1     CONTINUE
C...
C...  INITIAL DERIVATIVES
      CALL DERV
      IP=0
      RETURN
      END
      SUBROUTINE DERV2
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=161)
      COMMON/T/        T,     NSTOP,     NORUN
     1      /Y/    U(NX)
     2      /F/   UT(NX)
     3      /S/   UX(NX),   UXX(NX)
     4      /C/    X(NX),         N,      NORD,         IP
C...
C...  FIRST DERIVATIVE UX INCLUDING THE DIRCHLET BOUNDARY CONDITION
C...  AT X = 1
      U (N)= DEXP(-T)
      UT(N)=-DEXP(-T)
      CALL DSS004(0.D0,1.D0,N,U,UX)
C...
C...  SECOND DERIVATIVE INCLUDING THE SYMMETRY (NEUMANN) BOUNDARY
C...  CONDITION AT X = 0
      NL=2
      NU=1
      UX(1)=0.D0
      CALL DSS044(0.D0,1.D0,N,U,UX,UXX,NL,NU)
C...
C...  PDE
      DO 1 I=1,N-1
C...     X = 0
         IF(I.EQ.1)THEN
            UT(1)=(1.D0/U(1))*3.D0*(U(1)*UXX(1)+UX(1)**2)+5.D0*U(1)
C...
C...     0 LT X LT 1
         ELSE
            UT(I)=(1.D0/U(I))*(U(I)*UXX(I)+UX(I)**2+2.D0/X(I)*U(I)*UX(I)
     1           +5.D0*U(I)**2+4.D0*X(I)*U(I)*UX(I))
         END IF
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT2(NI,NO)
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
      PARAMETER (NX=161)
      COMMON/T/        T,     NSTOP,     NORUN
     1      /Y/    U(NX)
     2      /F/   UT(NX)
     3      /S/   UX(NX),   UXX(NX)
     4      /C/    X(NX),         N,      NORD,         IP
C...
C...  ARRAYS FOR PLOTTING
      DIMENSION SLOGN(11), SLOGE(11)
C...
C...  MONITOR THE CALCULATION
      WRITE(*,*)NORD,NORUN,N
C...
C...  COMPUTE THE EXACT SOLUTION, THE DIFFERENCE BETWEEN THE NUMERICAL
C...  AND ANALYTICAL SOLUTIONS, AND SAVE THE LARGEST NUMERICAL ERROR
      ERRMAX=0.D0
      DO 1 I=1,N-1
         UA=DEXP(1.0D0-X(I)**2-T)
         ERR=DABS(UA-U(I))
         IF(ERRMAX.LT.ERR)ERRMAX=ERR
1     CONTINUE
C...
C...  LABEL AND PRINT THE LARGEST ERROR (NOTE - AT T = 0, THE ERROR
C...  IS ZERO)
      IP=IP+1
      IF(IP.EQ.1)THEN
C...
C...     PRINT A HEADING AT T = 0
         WRITE(NO,3)
3        FORMAT(//,8X,' SOLUTION BY FOURTH ORDER APPROXIMATIONS',//,
     1   19X,'T',4X,'N',9X,'ERRMAX')
         ERRPLT=ERRMAX
      ELSE IF(IP.NE.1)THEN
C...
C...     PRINT ERRMAX VS T
         WRITE(NO,2)T,N,ERRMAX
2        FORMAT(10X,F10.3,I5,D15.4)
         IF(ERRMAX.GT.ERRPLT)ERRPLT=ERRMAX
      END IF
C...
C...  STORE LOG10(ERRPLT) VS LOG10(N) FOR SUBSEQUENT PLOTTING
      IF(IP.EQ.11)THEN
         SLOGN(NORUN)=DLOG10(DFLOAT(N))
         SLOGE(NORUN)=DLOG10(ERRPLT)
C...
C...  AT THE END OF THE FIFTH RUN (WITH N = 161), PLOT LOG10(ERR)
C...  VS LOG10(N)
      IF(NORUN.EQ.5)THEN
         WRITE(NO,4)(SLOGN(I),SLOGE(I),I=1,NORUN)
4        FORMAT(//,19X,' DATA FOR PLOT',//,
     1          15X,'LOG10(N)',6X,'LOG10(ERR)',/,
     2          (7X,2F15.4))
         CALL SPLOTS(1,NORUN,SLOGN,SLOGE)
         WRITE(NO,5)
5        FORMAT(//,' LOG10(MAXERR) VS LOG10(N)')
      END IF
      END IF
      RETURN
      END
