SUBROUTINE PDEONE(T, U, UDOT, NPDE, NPTS)                         PDE   10
      DIMENSION U(NPDE,NPTS), UDOT(NPDE,NPTS)
C PDEONE IS AN INTERFACE SUBROUTINE WHICH USES CENTERED DIF-
C FERENCE APPROXIMATIONS TO CONVERT ONE-DIMENSIONAL SYSTEMS
C OF PARTIAL DIFFERENTIAL EQUATIONS INTO A SYSTEM OF ORDINARY
C DIFFERENTIAL EQUATIONS, UDOT = F(T,X,U).  THIS ROUTINE IS
C INTENDED TO BE USED WITH A ROBUST ODE INTEGRATOR.
C INPUT..
C  NPDE = NUMBER OF PARTIAL DIFFERENTIAL EQUATIONS.
C  NPTS = NUMBER OF SPATIAL GRID POINTS.
C     T = CURRENT VALUE OF TIME.
C     U = AN NPDE BY NPTS ARRAY CONTAINING THE COMPUTED
C         SOLUTION AT TIME T.
C OUTPUT..
C  UDOT = AN NPDE BY NPTS ARRAY CONTAINING THE RIGHT HAND
C         SIDE OF THE RESULTING SYSTEM OF ODE*S, F(T,X,U),
C         OBTAINED BY DISCRETIZING THE GIVEN PDE*S.
C THE USER MUST INSERT A DIMENSION STATEMENT OF THE FOLLOW-
C ING FORM..
C     DIMENSION DVAL(**,**,2),UX(**),UAVG(**),ALPHA(**),
C    *          BETA(**),GAMMA(**)
C THE SYMBOLS ** DENOTE THE ACTUAL NUMERICAL VALUE OF NPDE
C FOR THE PROBLEM BEING SOLVED.
C COMMON BLOCK MESH CONTAINS THE USER SPECIFIED SPATIAL
C GRID POINTS.
C COMMON BLOCK COORD CONTAINS 0,1, OR 2 DEPENDING ON WHETHER
C THE PROBLEM IS IN CARTESIAN, CYLINDRICAL, OR SPHERICAL
C COORDINATES, RESPECTIVELY.
      COMMON /MESH/ X(1)
      COMMON /COORD/ ICORD
      ICORD1 = ICORD + 1
C UPDATE THE LEFT BOUNDARY VALUES
      CALL BNDRY(T, X(1), U, ALPHA, BETA, GAMMA, NPDE)
      ITEST = 0
      DXI = 1./(X(2)-X(1))
      DO 10 K=1,NPDE
        IF (BETA(K).NE.0.0) GO TO 10
        U(K,1) = GAMMA(K)/ALPHA(K)
        ITEST = ITEST + 1
   10 CONTINUE
      IF (ITEST.EQ.NPDE) GO TO 50
      IF (ITEST.EQ.0) GO TO 20
      CALL BNDRY(T, X(1), U, ALPHA, BETA, GAMMA, NPDE)
C EVALUATE D COEFFICIENTS AT THE LEFT BOUNDARY
   20 CALL D(T, X(1), U, DVAL, NPDE)
C FORM APPROXIMATIONS TO DU/DX AT THE LEFT BOUNDARY
      DO 40 K=1,NPDE
        IF (BETA(K).NE.0.0) GO TO 30
        UX(K) = DXI*(U(K,2)-U(K,1))
        GO TO 40
   30   UX(K) = (GAMMA(K)-ALPHA(K)*U(K,1))/BETA(K)
   40 CONTINUE
C EVALUATE U-AVERAGE IN THE FIRST INTERVAL
   50 DO 60 K=1,NPDE
        UAVG(K) = .5*(U(K,2)+U(K,1))
   60 CONTINUE
C EVALUATE D COEFFICIENTS IN THE FIRST INTERVAL
      XAVGR = .5*(X(2)+X(1))
      CALL D(T, XAVGR, UAVG, DVAL(1,1,2), NPDE)
      DXIL = 1.
      DXIR = DXI
      IF (ICORD.EQ.0) GO TO 70
      DXIL = X(1)**ICORD
      DXIR = DXIR*XAVGR**ICORD
C EVALUATE DUXX AT THE LEFT BOUNDARY
   70 IF (ITEST.EQ.NPDE) GO TO 100
      DXIC = FLOAT(ICORD1)/(XAVGR**ICORD1-X(1)**ICORD1)
      DO 90 L=1,NPDE
        DO 80 K=1,NPDE
          DVAL(K,L,1) = DXIC*(DVAL(K,L,2)*(U(L,2)-U(L,1))*
     *     DXIR-DVAL(K,L,1)*UX(L)*DXIL)
   80   CONTINUE
   90 CONTINUE
C EVALUATE RIGHTHAND SIDE OF PDE*S AT THE LEFT BOUNDARY
      CALL F(T, X(1), U, UX, DVAL, UDOT, NPDE)
C SET UDOT = 0 FOR KNOWN LEFT BOUNDARY VALUES
  100 DO 110 K=1,NPDE
        IF (BETA(K).EQ.0.0) UDOT(K,1) = 0.0
  110 CONTINUE
C UPDATE THE RIGHT BOUNDARY VALUES
      CALL BNDRY(T, X(NPTS), U(1,NPTS), ALPHA, BETA, GAMMA, NPDE)
      ITEST = 0
      DO 120 K=1,NPDE
        IF (BETA(K).NE.0.0) GO TO 120
        U(K,NPTS) = GAMMA(K)/ALPHA(K)
        ITEST = ITEST + 1
  120 CONTINUE
      IBCK = 1
      IFWD = 2
      ILIM = NPTS - 1
C MAIN LOOP TO FORM ORDINARY DIFFERENTIAL EQUATIONS AT THE
C GRID POINTS
      DO 160 I=2,ILIM
        K = IBCK
        IBCK = IFWD
        IFWD = K
        XAVGL = XAVGR
        XAVGR = .5*(X(I+1)+X(I))
        DXI = 1./(X(I+1)-X(I-1))
        DXIL = DXIR
        DXIR = 1./(X(I+1)-X(I))
        IF (ICORD.NE.0) DXIR = DXIR*XAVGR**ICORD
        DXIC = FLOAT(ICORD1)/(XAVGR**ICORD1-XAVGL**ICORD1)
C EVALUATE DU/DX AND U-AVERAGE AT THE I-TH GRID POINT AND
C INTERVAL RESPECTIVELY
        DO 130 K=1,NPDE
          UX(K) = DXI*(U(K,I+1)-U(K,I-1))
          UAVG(K) = .5*(U(K,I+1)+U(K,I))
  130   CONTINUE
C EVALUATE D COEFFICIENTS IN THE I-TH INTERVAL
        CALL D(T, XAVGR, UAVG, DVAL(1,1,IFWD), NPDE)
C EVALUATE DUXX AT THE I-TH GRID POINT
        DO 150 L=1,NPDE
          DO 140 K=1,NPDE
            DVAL(K,L,IBCK) = (DVAL(K,L,IFWD)*(U(L,I+1)-U(L,I))*
     *       DXIR-DVAL(K,L,IBCK)*(U(L,I)-U(L,I-1))*DXIL)*DXIC
  140     CONTINUE
  150   CONTINUE
C EVALUATE RIGHTHAND SIDE OF PDE*S AT THE I-TH GRID POINT
        CALL F(T, X(I), U(1,I), UX, DVAL(1,1,IBCK), UDOT(1,I), NPDE)
  160 CONTINUE
C FINISH UPDATING THE RIGHT BOUNDARY IF NECESSARY
      IF (ITEST.EQ.NPDE) GO TO 220
      IF (ITEST.EQ.0) GO TO 170
      CALL BNDRY(T, X(NPTS), U(1,NPTS), ALPHA, BETA, GAMMA, NPDE)
C EVALUATE D COEFFICIENTS AT THE RIGHT BOUNDARY
  170 CALL D(T, X(NPTS), U(1,NPTS), DVAL(1,1,IBCK), NPDE)
C FORM APPROXIMATIONS TO DU/DX AT THE RIGHT BOUNDARY
      DXI = 1./(X(NPTS)-X(ILIM))
      DO 190 K=1,NPDE
        IF (BETA(K).NE.0.0) GO TO 180
        UX(K) = DXI*(U(K,NPTS)-U(K,ILIM))
        GO TO 190
  180   UX(K) = (GAMMA(K)-ALPHA(K)*U(K,NPTS))/BETA(K)
  190 CONTINUE
      DXIL = DXIR
      DXIR = 1.
      IF (ICORD.NE.0) DXIR = X(NPTS)**ICORD
      DXIC = FLOAT(ICORD1)/(X(NPTS)**ICORD1-XAVGR**ICORD1)
C EVALUATE DUXX AT THE RIGHT BOUNDARY
      DO 210 L=1,NPDE
        DO 200 K=1,NPDE
          DVAL(K,L,IBCK) = DXIC*(DVAL(K,L,IBCK)*UX(L)*DXIR-DVAL(K,L,
     *     IFWD)*(U(L,NPTS)-U(L,ILIM))*DXIL)
  200   CONTINUE
  210 CONTINUE
C EVALUATE RIGHTHAND SIDE OF PDE*S AT THE RIGHT BOUNDARY
      CALL F(T, X(NPTS), U(1,NPTS), UX, DVAL(1,1,IBCK),
     * UDOT(1,NPTS), NPDE)
C SET UDOT = 0 FOR KNOWN RIGHT BOUNDARY VALUES
  220 DO 230 K=1,NPDE
        IF (BETA(K).EQ.0.0) UDOT(K,NPTS) = 0.0
  230 CONTINUE
      RETURN
      END

.