      SUBROUTINE PSET(Y,N0,CON,MITER,IER,YMAX,SAVE1,SAVE2,PW,IPIV,
     1                  DIFFUN,PEDERV)
C-----------------------------------------------------------------------
C PSET IS CALLED BY TSTEP TO COMPUTE AND TO PROCESS THE MATRIX
C P = I - (H/EL(2))*J, WHERE J IS AN APPROXIMATION TO THE
C JACOBIAN.  J IS COMPUTED BY EITHER THE USER SUPPLIED
C SUBROUTINE PEDERV, WHEN MITER = 1, OR BY FINITE DIFFERENCES,
C WHEN MITER = 2.  J IS STORED IN PW AND REPLACED BY P, USING
C CON = -H/EL(2).  THEN P IS SUBJECTED TO AN LU DECOMPOSITION
C FOR LATER SOLUTION OF LINEAR ALGEBRAIC SYSTEMS WITH P AS THE
C COEFFICIENT MATRIX.
C
C IN ADDITION TO VARIABLES DESCRIBED PREVIOUSLY, COMMUNICATION
C WITH PSET USES THE FOLLOWING..
C   EPSJ = SQRT(UROUND), USED IN THE NUMERICAL JACOBIAN INCREMENTS.
C   NSQ  = N0**2.
C-----------------------------------------------------------------------
C CAUTION   NOT ALL EPCOM1 VARIABLES ARE USED INTHIS SUBROUTINE.
C-----------------------------------------------------------------------
      INTEGER IER, MITER, N0
      INTEGER IPIV, JSTART, KFLAG, MF, N, NSQ
      INTEGER I, IPJ1, J, J1
      REAL CON, Y
      REAL EPS, EPSJ, H, HMAX, HMIN, PW, SAVE1, SAVE2,
     1     SS, T, UROUND, YMAX
      REAL D, R, R0, YJ
      REAL ONE, REP, ZERO
      DIMENSION Y(N0,1)
      DIMENSION YMAX(1),SAVE1(1),SAVE2(1),PW(1),IPIV(1)
      EXTERNAL DIFFUN,PEDERV
C
      COMMON /EPCOM1/ T,H,HMIN,HMAX,EPS,SS,UROUND,N,MF,KFLAG,JSTART
      COMMON /EPCOM8/ EPSJ,NSQ
      DATA ONE /1.0E0/, REP /1.0E-3/, ZERO /0.0E0/
C
      IF (MITER .EQ. 2) GO TO 20
C IF MITER = 1, CALL PEDERV AND MULTIPLY BY A SCALAR.  -----------------
      CALL PEDERV (N, T, Y, PW, N0)
      DO 10 I = 1,NSQ
  10    PW(I) = PW(I)*CON
      GO TO 60
C IF MITER = 2, MAKE N CALLS TO DIFFUN TO APPROXIMATE J. ---------------
  20  D = ZERO
      DO 30 I = 1,N
  30    D = D + SAVE2(I)**2
      R0 = ABS(H)*SQRT(D)*UROUND/REP
      J1 = 0
      DO 50 J = 1,N
        YJ = Y(J,1)
        R = EPSJ*YMAX(J)
        R = AMAX1(R,R0)
        Y(J,1) = Y(J,1) + R
        D = CON/R
        CALL DIFFUN (N, T, Y, SAVE1)
        DO 40 I = 1,N
          IPJ1 = I + J1
  40      PW(IPJ1) = (SAVE1(I) - SAVE2(I))*D
        Y(J,1) = YJ
        J1 = J1 + N0
  50    CONTINUE
C ADD ON THE IDENTITY MATRIX.  -----------------------------------------
  60  J = 1
      DO 70 I = 1,N
        PW(J) = PW(J) + ONE
  70    J = J + (N0 + 1)
C GET LU DECOMPOSITION OF P. -------------------------------------------
      CALL DEC (N, N0, PW, IPIV, IER)
      RETURN
C-----------------------  END OF SUBROUTINE PSET -----------------------
      END
