        subroutine GRSUB4(Y,ERROR,WT,D,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 BND, T, HOLD, H, EPS, AHINV, A, HMIN, TOLD
      double precision D, D1, FOURU, R, RDWN, RUP
      double precision Y, ERROR, WT, ERUP, ERDWN
      integer NEQN, K, NFE, ICRASH, KOLD, ITST, IWEVAL, IFAIL, KP1
      integer J,I,KD
      logical START, CONVRG
      dimension Y(NEQN,7), ERROR(NEQN), WT(NEQN), ERUP(5), ERDWN(6)
      data ERUP/3.0,6.0,9.167,12.5,15.98/
      data ERDWN/1.0,1.0,0.5,0.1667,0.04133,0.008267/
      HOLD = H
      if (KP1.lt.3) go to 110
      do 100 J = 3, KP1
         do 100 I = 1, NEQN
            Y(I,J) = Y(I,J)+A(J)*ERROR(I)
100   continue
110   continue
      ITST = ITST-1
      if (ITST.eq.0) go to 130
      if ((ITST.gt.1).or.(K.eq.6)) go to 190
      do 120 I = 1, NEQN
         Y(I,7) = ERROR(I)
120   continue
      go to 190
130   continue
      KD = 0
      R = 1.2*(2.*D)**(1./float(KP1))
      if (K.eq.1) go to 150
      D1 = 0.0
      do 140 I = 1, NEQN
         D1 = dmax1(D1,dabs(Y(I,KP1)/WT(I)))
140   continue
      RDWN = 1.3*(2.*D1/(ERDWN(K)*EPS))**(1./float(K))
      if (RDWN.gt.R) go to 150
      R = RDWN
      KD = -1
150   continue
      if (K.eq.6) go to 180
      D1 = 0.
      do 160 I = 1, NEQN
         D1 = dmax1(D1,dabs((ERROR(I)-Y(I,7))/WT(I)))
160   continue
      RUP = 1.4*(2.*D1/(ERUP(K)*EPS))**(1./float(K+2))
      if (RUP.gt.R) go to 180
      KD = 1
      R = RUP
      D1 = A(KP1)/float(KP1)
      do 170 I = 1, NEQN
         Y(I,KP1+1) = ERROR(I)*D1
170   continue
180   continue
      ITST = KP1
      if ((KD.eq.0).and.(dabs(R-0.95).lt.0.06)) go to 190
      K = K+KD
      if (R.lt.0.5) R = 0.5
      H = dsign(dmax1(dabs(H/R),FOURU*dabs(T)),H)
190   continue
      return
      end
