*DECK TEST1D1
      SUBROUTINE INITAL
C...
C...  THE FOLLOWING ONE-DIMENSIONAL PARABOLIC PARTIAL DIFFERENTIAL
C...  EQUATION (PDE) WITH NONLINEAR BOUNDARY CONDITIONS IS INTEGRATED
C...  BY THE NUMERICAL METHOD OF LINES.  THE SPATIAL DERIVATIVES IN
C...  IN X ARE APPROXIMATED BY ORTHOGONAL COLLOCATION ON FINITE ELEMENTS
C...  AS IMPLEMENTED IN SUBROUTINE DSS046 WHICH IS CALLED IN SUBROUTINE
C...  DERV.
C...
C...  DSS046 IS A GENERAL DIFFERENTIATION ROUTINE FOR ONE-DIMENSIONAL
C...  DOMAINS.  THE USER CAN SPECIFY THE ORDER OF THE APPROXIMATING
C...  ORTHOGONAL POLYNOMIAL AND THE NUNBER OF FINITE ELEMENTS.
C...
C...  ADDITIONAL DETAILS CONCERNING THE IMPLEMENTATION OF ORTHOGONAL
C...  COLLOCATION ON FINITE ELEMENTS ARE GIVEN IN THE SUBORDINATE
C...  SUBROUTINES CALLED BY DSS046.
C...
C...     SINGLE PRECISION
C...
C...     U  = U
C...      T    XX
C...
C...     U (O,T) = 0
C...      X
C...
C...                    4         4
C...     U (1,T) = E*(UA  - U(1,T) )
C...      X
C...
C...     U(X,0) = 600
C...
C...     UA = 500,   E = 1.73E-09
C...
      COMMON      /T/      T,       NFIN,      NRUN
     1            /Y/   U(4)
     2            /F/  UT(4)
     3           /UX/  UX(4),     UXX(4)
     4            /C/      N,         UA
     5            /B/     XL,         XU,       X(4)
     6            /O/    NEL,       NDEG
C...
C...  SET THE NUMBER OF ELEMENTS (UNITY FOR ORTHOGONAL COLLOCATION)
      NEL=1
C...
C...  SPECIFY THE DEGREE OF THE ORTHOGONAL POLYNOMIAL (CUBIC LOBATTO
C...  QUADRATURE IS SELECTED)
      NDEG=4
C...
C...  SET THE NUMBER OF NODAL POINTS (= NEL*(NDEG-1)+1)
      N=NEL*(NDEG-1)+1
C...
C...  SET THE BOUNDARY VALUES OF X
      XL=0.E0
      XU=1.E0
C...
C...  SET THE AMBIENT TEMPERATURE, INITIAL TEMPERATURE
      UA=500.E0
      DO 1 I=1,N
1     U(I)=600.E0
      RETURN
      END
      SUBROUTINE DERV
      COMMON      /T/      T,       NFIN,      NRUN
     1            /Y/   U(4)
     2            /F/  UT(4)
     3           /UX/  UX(4),     UXX(4)
     4            /C/      N,         UA
     5            /B/     XL,         XU,       X(4)
     6            /O/    NEL,       NDEG
C...
C...  COMPUTE THE PARTIAL DERIVATIVE UX
      CALL DSS046(NEL,NDEG,N,XL,XU,U,UX)
C...
C...  APPLY THE BOUNDARY CONDITION AT X = XL AS A CONSTRAINT
      UX(1)=0.E0
C...
C...  APPLY THE BOUNDARY CONDITION AT X = XU AS A CONSTRAINT
      UX(N)=1.73E-09*(UA**4-U(N)**4)
C...
C...  COMPUTE THE PARTIAL DERIVATIVE UXX
      CALL DSS046(NEL,NDEG,N,XL,XU,UX,UXX)
C...
C...  ASSEMBLE THE PARTIAL DIFFERENTIAL EQUATION
      DO 1 I=1,N
      UT(I)=UXX(I)
1     CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON      /T/      T,       NFIN,      NRUN
     1            /Y/   U(4)
     2            /F/  UT(4)
     3           /UX/  UX(4),     UXX(4)
     4            /C/      N,         UA
     5            /B/     XL,         XU,       X(4)
     6            /O/    NEL,       NDEG
C...
C...  ARRAY UP CONTAINS THE INTERPOLATED DEPENDENT VARIABLE
      DIMENSION UP(6)
C...
C...  ARRAY XP CONTAINS THE VALUES OF X AT WHICH THE DEPENDENT VARIABLE
C...  IS INTERPOLATED
      DIMENSION XP(6)
C...
C...  ARRAY D CONTAINS THE NEWTON DIVIDED DIFFERENCES FOR INTERPOLATION
      DIMENSION D(4)
C...
C...  AT THE BEGINNING OF THE RUN, SET UP THE PARAMETERS FOR THE
C...  ORTHOGONAL COLLOCATION SOLUTION
      IF(T.EQ.0.E0)THEN
C...
C...  STORE THE PRECOMPUTED COLLOCATION POINTS IN ARRAY X
      CALL DSS46R(NEL,NDEG,N,XL,XU,X)
C...
C...  SPECIFY THE NUMBER OF POINTS AT WHICH THE DEPENDENT VARIABLE IS
C...  INTERPOLATED
      NP=6
C...
C...  COMPUTE THE POINTS IN X FOR THE INTERPOLATION
      DO 1 I=1,NP
1     XP(I)=FLOAT(I-1)/FLOAT(NP-1)
      END IF
C...
C...  COMPUTE THE NEWTON DIVIDED DIFFERENCES
      CALL INTERP(N,X,U,D)
C...
C...  INTERPOLATE THE SOLUTION
      DO 2 J=1,NP
      UP(J)=POLNST(N,XP(J),X,D)
2     CONTINUE
C...
C...  PRINT THE SOLUTION
      WRITE(NO,3)T
3     FORMAT(/,2X,7HTIME = ,E9.2,/,1X,
     111H     X=0   ,11H    X=0.2  ,11H    X=0.4  ,
     211H    X=0.6  ,11H    X=0.8  ,11H    X=1.0  )
      WRITE(NO,4)(UP(I),I=1,NP)
4     FORMAT(6E11.3)
      RETURN
      END

