      subroutine STIFF(Y,YP,IPS,YY,W,WORK,FOURU)
      common /SFDATA/ T,TOUT,RELERR,ABSERR,DTSGN,NEQN,IFLAG,
     &                 IOLD,INIT,KOP
      common /GEDATA/ X,HOLD,H,EPS,RVAR(8),NEQ,K,NFE,
     &                ICRASH,KOLD,ITST,IWEVAL,START
      double precision ABSDT, ABSEPS, ABSERR, DT, DTSGN,
     &                 EPS, H, HOLD, RELEPS, RELERR, RVAR,
     &                 T, TEND, TOUT, W, WORK, X, Y, YP, YY
      double precision FOURU
      integer NEQN,MAXNFE,MFLAG,IFLAG,NEQ,INIT,IOLD,NFE,KOP,L,IPS
      integer ICRASH,K,KOLD,ITST,IWEVAL
      logical START
      dimension Y(NEQN),YP(NEQN),IPS(NEQN),YY(NEQN,7),
     &W(NEQN,NEQN),WORK(NEQN,9)
      data MAXNFE/1000/
      MFLAG = iabs(IFLAG)
      if ((MFLAG.gt.7).or.(MFLAG.eq.0).or.(NEQN.lt.1)) go to 180
      if ((RELERR.lt.0.0).or.(ABSERR.lt.0.0)) go to 180
      EPS = dmax1(RELERR,ABSERR)
      if (EPS.le.0.0) go to 180
      NEQ = NEQN
      DT = TOUT-T
      ABSDT = dabs(DT)
      TEND = T+10.*DT
      if (IFLAG.lt.0) TEND = TOUT
      RELEPS = RELERR/EPS
      ABSEPS = ABSERR/EPS
      if (MFLAG.eq.1) go to 140
      if (T.eq.TOUT) go to 180
      if (INIT.eq.0) go to 170
      if (IOLD.lt.0) go to 150
      if (DTSGN*DT.lt.0.0) go to 150
100   continue
      if (dabs(X-T).lt.ABSDT) go to 110
      call GINTRP (TOUT,Y,YP,YY)
      IOLD = IFLAG
      IFLAG = 2
      T = TOUT
      go to 290
110   continue
      if ((IFLAG.lt.0).and.(dabs(TOUT-X).lt.FOURU*dabs(X))) go to 190
      if (NFE.gt.MAXNFE) go to 210
      if (dabs(H).lt.dabs(TEND-X)) go to 120
      KOP = KOP+1
      if (KOP.gt.50) go to 240
      H = TEND-X
120   continue
      do 130 L = 1, NEQN
         Y(L) = RELEPS*dabs(YY(L,1))+ABSEPS
         if (Y(L).eq.0.0) go to 230
130   continue
      call GEAR (YY,W,IPS,Y,YP,WORK(1,1),WORK(1,2),WORK(1,3),FOURU)
      if (ICRASH.ne.0) go to 220
      go to 100
140   continue
      NFE = 0
150   continue
      START = .true.
      INIT = 0
      KOP = 0
      X = T
      do 160 L = 1, NEQN
         YY(L,1) = Y(L)
160   continue
      if (T.ne.TOUT) go to 170
      call FCT (X,YY,YP)
      NFE = NFE+1
      IFLAG = 2
      go to 290
170   continue
      INIT = 1
      DTSGN = DT
      H = dsign(dmax1(dabs(TOUT-X),FOURU*dabs(X)),TOUT-X)
      go to 100
180   continue
      IFLAG = 9
      go to 290
190   continue
      H = TOUT-X
      call FCT (X,YY,YP)
      NFE = NFE+1
      do 200 L = 1, NEQN
         Y(L) = YY(L,1)+H*YP(L)
200   continue
      IFLAG = 2
      T = TOUT
      call FCT (T,Y,YP)
      NFE = NFE+1
      IOLD = -1
      go to 290
210   continue
      IFLAG = isign(5,IFLAG)
      NFE = NFE-MAXNFE
      go to 250
220   continue
      IFLAG = isign(4,IFLAG)
      if (ICRASH.eq.2) IFLAG = isign(3,IFLAG)
      if (ICRASH.eq.3) IFLAG = isign(6,IFLAG)
      RELERR = EPS*RELEPS
      ABSERR = EPS*ABSEPS
      if (ICRASH.eq.1) RELERR = RELERR*3.0
      if (ICRASH.eq.1) ABSERR = ABSERR*3.0
      go to 250
230   continue
      IFLAG = 8
      go to 250
240   continue
      KOP = 0
      IFLAG = isign(7,IFLAG)
250   continue
      T = X
      if (HOLD.eq.0.0) go to 270
      IOLD = 1
      do 260 L = 1, NEQN
         Y(L) = YY(L,1)
         YP(L) = YY(L,2)/HOLD
260   continue
      go to 290
270   continue
      do 280 L = 1, NEQN
         Y(L) = YY(L,1)
280   continue
      call FCT (T,Y,YP)
      IOLD = -1
290   continue
      return
      end
