*DECK SEP16
      SUBROUTINE INITAL
C...
C...  NONLINEAR DIFFUSION
C...
C...  THE ONE-DIMENSIONAL DIFFUSION EQUATION WITH TWO NONLINEAR COEFFI-
C...  CIENTS
C...
C...          (D(U)U )  = C(U)U                                      (1)
C...                X X        T
C...
C...  IS TO BE INTEGRATED SUBJECT TO THE AUXILIARY CONDITIONS
C...
C...          U (0,T) = 0, U(1,T) = 1                             (2)(3)
C...           X
C...
C...          U(X,0) = 0                                             (4)
C...
C...  THE TWO NONLINEAR COEFFICIENTS IN THIS CASE ARE
C...
C...          D(U) = C(U) = 1 + U                                 (5)(6)
C...
C...  THE FOLLOWING NUMERICAL SOLUTION WAS REPORTED BY ZELLNER FOR
C...  COMPARISON WITH THE METHOD OF LINES SOLUTION
C...
C...          U(0.0,0.2)     0.29738
C...
C...          U(0.2,0.2)     0.33817
C...
C...          U(0.4,0.2)     0.45173
C...
C...          U(0.6,0.2)     0.61651
C...
C...          U(0.8,0.2)     0.80710
C...
C...          U(1.0,0.2)     1.00000
C...
C...  (ZELLNER, M. G., PH.D. DISSERTATION, LEHIGH UNIVERSITY, BETHLEHEM,
C...  PA).  EQUATION (1) CAN ALSO BE LINEARIZED BY A KIRCHOFF TRANS-
C...  FORMATION AND THE RESULTING LINEAR EQUATION SOLVED ANALYTICALLY.
C...  THIS LINEARIZATION THEREFORE PROVIDES A RIGOROUS TEST OF THE
C...  NUMERICAL SOLUTION.  THE DETAILS ARE GIVEN BY ZELLNER.
C...
C...  EQUATION (1) IS INTEGRATED WITH THREE AND FIVE-POINT CENTERED
C...  DIFFERENCES.  ALSO, THE LEFT HAND SIDE OF EQUATION (1) CAN BE
C...  EXPANDED TO GIVE
C...
C...          D(U)U   + D(U) U  = C(U)U                              (7)
C...               XX       X X        T
C...
C...  EQUATION (7) IS ALSO INTEGRATED BY THREE AND FIVE-POINT CENTERED
C...  DIFFERENCES SO THAT A TOTAL OF FOUR RUNS IS PROGRAMMED.  A COMPAR-
C...  ISON OF THE NUMERICAL SOLUTIONS FOR EQUATIONS (1) AND (7) WILL
C...  GIVE SOME INDICATION WHETHER THE PARTICULAR FORM OF THE DIFFER-
C...  ERENTIAL EQUATION AFFECTS THE SOLUTION SIGNIFICANTLY.
C...
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(11)
     2      /F/ UT(11)
     3     /SD/ UX(11),UXX(11)
     4      /P/      N,     IP
C...
C...  SET THE NUMBER OF GRID POINTS, COUNTER FOR THE PLOTTED SOLUTION
      N=11
      IP=0
C...
C...  INITIAL CONDITION (4)
      DO 1 I=1,N
      U(I)=0.
1     CONTINUE
C...
C...  CALL SUBROUTINE DERV TO INITIALIZE ALL OF THE DERIVATIVES
      CALL DERV
      RETURN
      END
      SUBROUTINE DERV
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(11)
     2      /F/ UT(11)
     3     /SD/ UX(11),UXX(11)
     4      /P/      N,     IP
C...
C...  BOUNDARY CONDITION (3)
      U(N)=1.
C...
C...  FOUR RUNS ARE PROGRAMMED AS OUTLINED IN SUBROUTINE INITAL
      GO TO(1,2,3,4),NORUN
C...
C...  EQUATION (1) WITH THREE-POINT CENTERED DIFFERENCES
1     CALL DSS002(0.,1.,N,U ,UX )
      UX(1)=0.
      DO 11 I=1,N
      UX(I)=(1.+U(I))*UX(I)
11    CONTINUE
      CALL DSS002(0.,1.,N,UX,UXX)
      DO 12 I=1,N
      UT(I)=(1./(1.+U(I)))*UXX(I)
12    CONTINUE
      RETURN
C...
C...  EQUATION (1) WITH FIVE-POINT CENTERED DIFFERENCES
2     CALL DSS004(0.,1.,N,U ,UX)
      UX(1)=0.
      DO 21 I=1,N
      UX(I)=(1.+U(I))*UX(I)
21    CONTINUE
      CALL DSS004(0.,1.,N,UX,UXX)
      DO 22 I=1,N
      UT(I)=(1./(1.+U(I)))*UXX(I)
22    CONTINUE
      RETURN
C...
C...  EQUATION (7) WITH THREE-POINT CENTERED DIFFERENCES
3     CALL DSS002(0.,1.,N,U ,UX )
      UX(1)=0.
      CALL DSS004(0.,1.,N,UX,UXX)
      DO 31 I=1,N
      UT(I)=(1./(1.+U(I)))*((1.+U(I))*UXX(I)+UX(I)**2)
31    CONTINUE
      RETURN
C...
C...  EQUATION (7) WITH FIVE-POINT CENTERED DIFFERENCES
4     CALL DSS004(0.,1.,N,U ,UX )
      UX(1)=0.
      CALL DSS004(0.,1.,N,UX,UXX)
      DO 41 I=1,N
      UT(I)=(1./(1.+U(I)))*((1.+U(I))*UXX(I)+UX(I)**2)
41    CONTINUE
      RETURN
      END
      SUBROUTINE PRINT(NI,NO)
      COMMON/T/      T,  NSTOP,  NORUN
     1      /Y/  U(11)
     2      /F/ UT(11)
     3     /SD/ UX(11),UXX(11)
     4      /P/      N,     IP
      IP=IP+1
C...
C...  PRINT A HEADING FOR THE NUMERICAL SOLUTION
      IF(IP.EQ.1)WRITE(NO,1)
1     FORMAT(1H ,9H        T,3X,9H U(T,0.0),3X,9H U(T,0.2),3X,
     1                          9H U(T,0.4),3X,9H U(T,0.6),3X,
     2                          9H U(T,0.8),3X,9H U(T,1.0))
C...
C...  PRINT THE NUMERICAL SOLUTION
      WRITE(NO,2)T,(U(I),I=1,N,2)
2     FORMAT(1X,F9.3,6(3X,F9.6))
      RETURN
      END

