      subroutine GRSUB2(Y,W,IPS,WT,YP,C,ERROR,SAVER,FOURU)
      common/GEDATA/T,HOLD,H,EPS,AHINV,A(7),NEQN,K,NFE,
     &ICRASH,KOLD,ITST,IWEVAL,START
      common/GEDAT2/BND,HMIN,TOLD,IFAIL,KP1,CONVRG
      double precision T, HOLD, H, EPS, AHINV, A
      double precision HMIN, TOLD, BND
      double precision CNORM, FOURU, SUM, GFACT
      double precision Y, YP, WT, W, C, ERROR, SAVER
      integer NEQN, K, NFE, ICRASH, KOLD, ITST,IWEVAL
      integer IFAIL,KP1,J,J1,J2,I,ITER,IP,IM1,IBACK,IP1,IPS
      logical CONVRG, START
      dimension Y(NEQN,7), YP(NEQN), WT(NEQN), W(NEQN,NEQN), IPS(NEQN),
     &          C(NEQN), ERROR(NEQN), SAVER(NEQN,7)
      CONVRG = .false.
100   continue
      T = T+H
      do 110 J = 2, KP1
         do 110 J1 = J, KP1
            J2 = K+J-J1
            do 110 I = 1, NEQN
               Y(I,J2) = Y(I,J2)+Y(I,J2+1)
110   continue
      do 120 I = 1, NEQN
         ERROR(I) = 0.0
120   continue
      ITER = 0
      if (IWEVAL.lt.1) go to 150
      call FCT (T,Y,YP)
      NFE = NFE+1
      call FDER (T,Y,W,YP)
      IWEVAL = -1
      if (NEQN.eq.1) go to 140
      AHINV = 1./(A(1)*H)
      do 130 I = 1, NEQN
         W(I,I) = W(I,I)+AHINV
130   continue
      call DECOMP (NEQN,W,IPS,C,J1)
      if (J1.lt.0) go to 250
      go to 160
140   continue
      W(1,1) = 1.+W(1,1)*A(1)*H
      if (W(1,1).eq.0.0) go to 250
      go to 210
150   continue
      call FCT (T,Y,YP)
      NFE = NFE+1
      if (NEQN.eq.1) go to 210
160   continue
      IP = IPS(1)
      C(1) = (Y(IP,2)-YP(IP)*H)*AHINV
      do 180 I = 2, NEQN
         IP = IPS(I)
         IM1 = I-1
         SUM = 0.0
         do 170 J = 1, IM1
            SUM = SUM+W(IP,J)*C(J)
170      continue
         C(I) = (Y(IP,2)-YP(IP)*H)*AHINV-SUM
180   continue
      IP = IPS(NEQN)
      C(NEQN) = C(NEQN)/W(IP,NEQN)
      do 200 IBACK = 2, NEQN
         I = NEQN-IBACK+1
         IP = IPS(I)
         IP1 = I+1
         SUM = 0.0
         do 190 J = IP1, NEQN
            SUM = SUM+W(IP,J)*C(J)
190      continue
         C(I) = (C(I)-SUM)/W(IP,I)
200   continue
      go to 220
210   continue
      C(1) = (Y(1,2)-YP(1)*H)/W(1,1)
220   continue
      CNORM = 0.
      do 230 I = 1, NEQN
         Y(I,1) = Y(I,1)+A(1)*C(I)
         Y(I,2) = Y(I,2)-C(I)
         ERROR(I) = ERROR(I)+C(I)
         CNORM = dmax1(CNORM,dabs(C(I)/WT(I)))
230   continue
      if (CNORM.gt.BND) go to 240
      CONVRG = .true.
      go to 290
240   continue
      ITER = ITER+1
      if (ITER.lt.3) go to 150
250   continue
      T = TOLD
      if (IWEVAL.eq.0) go to 270
      H = 0.25*H
      if (dabs(H).lt.HMIN) go to 270
      do 260 I = 1, NEQN
         Y(I,1) = SAVER(I,1)
260   continue
      CONVRG = .false.
      go to 290
270   continue
      do 280 J = 1, KP1
         do 280 I = 1, NEQN
            Y(I,J) = SAVER(I,J)
280   continue
      IWEVAL = 1
      if (dabs(H).ge.HMIN) go to 100
      GFACT = CNORM/BND
      if (GFACT.lt.2.0) GFACT = 2.
      EPS = EPS*GFACT*(1.+FOURU)
      H = dsign(HMIN,H)
      ICRASH = 3
      go to 290
290   continue
      return
      end
