      subroutine STEP(Y,YP,WT,PHI,P,FOURU)
      double precision ABSH,ALPHA,BETA,EPS,
     &                 ERK,ERKM1,ERKM2,ERKP1,err,G,GSTR,
     &                 H,HNEW,HOLD,P,PHI,PSI,R,REALNS,RHO,ROUND,
     &                 SIG,SUM,TAU,TEMP1,TEMP2,TEMP3,
     &                 TEMP4,TEMP5,TEMP6,TWO,V,W,WT,X,XOLD,Y,YP
      double precision FOURU
      integer NEQN,ICRASH,L,IFAIL,NFE,K,KOLD,KP1,KP2,KM1,KM2,NS,NSP1
     &,I,IM1,IQ,NSM2,J,LIMIT1,NSP2,LIMIT2,IP1,KNEW
      logical START,PHASE1,NORND
      common /STDATA/ X,H,HOLD,EPS,ALPHA(12),BETA(12),G(13),
     &                SIG(13),V(12),W(12),PSI(12),NEQN,KOLD,
     &                ICRASH,K,NFE,NS,NORND,PHASE1,START
      dimension Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),
     &          YP(NEQN),GSTR(13),TWO(12)
      data TWO/8.,16.,32.,64.,128.,256.,512.,1024.,2048.,4096.,
     &         8192.,16384./
      data GSTR/.5,.0833,.0417,.0264,.0188,.0143,.0114,
     &          .00936,.00789,.00697,.00592,.00524,.00468/
      if (dabs(H).ge.FOURU*dabs(X)) go to 100
      H = dsign(FOURU*dabs(X),H)
      ICRASH = 1
      return
100   continue
      ROUND = 0.
      do 110 L = 1, NEQN
         ROUND = dmax1(ROUND,dabs(Y(L)/WT(L)))
110   continue
      ROUND = FOURU*ROUND
      if (EPS.ge.ROUND) go to 120
      EPS = ROUND*(1.+FOURU)
      ICRASH = 2
      return
120   continue
      ICRASH = 0
      G(1) = 1.0
      G(2) = 0.5
      SIG(1) = 1.0
      IFAIL = 0
      if (.not.START) go to 150
      call FCT (X,Y,YP)
      NFE = NFE+1
      SUM = 0.
      do 130 L = 1, NEQN
         PHI(L,1) = YP(L)
         PHI(L,2) = 0.
         SUM = dmax1(SUM,dabs(YP(L)/WT(L)))
130   continue
      ABSH = dabs(H)
      if (EPS.lt.16.*SUM*H*H) ABSH = 0.25*dsqrt(EPS/SUM)
      H = dsign(dmax1(ABSH,FOURU*dabs(X)),H)
      HOLD = 0.
      K = 1
      KOLD = 0
      START = .false.
      PHASE1 = .true.
      NORND = .true.
      if (EPS.gt.200.*ROUND) go to 150
      NORND = .false.
      do 140 L = 1, NEQN
         PHI(L,15) = 0.
140   continue
150   continue
      KP1 = K+1
      KP2 = K+2
      KM1 = K-1
      KM2 = K-2
      if (H.ne.HOLD) NS = 0
      NS = min0(NS+1,KOLD+1)
      NSP1 = NS+1
      if (K.lt.NS) go to 260
      BETA(NS) = 1.
      REALNS = NS
      ALPHA(NS) = 1./REALNS
      TEMP1 = H*REALNS
      SIG(NSP1) = 1.
      if (K.lt.NSP1) go to 170
      do 160 I = NSP1, K
         IM1 = I-1
         TEMP2 = PSI(IM1)
         PSI(IM1) = TEMP1
         BETA(I) = BETA(IM1)*PSI(IM1)/TEMP2
         TEMP1 = TEMP2+H
         ALPHA(I) = H/TEMP1
         SIG(I+1) = float(I)*ALPHA(I)*SIG(I)
160   continue
170   continue
      PSI(K) = TEMP1
      if (NS.gt.1) go to 190
      do 180 IQ = 1, K
         TEMP3 = IQ*(IQ+1)
         V(IQ) = 1./TEMP3
         W(IQ) = V(IQ)
180   continue
      go to 230
190   continue
      if (K.le.KOLD) go to 210
      TEMP4 = K*KP1
      V(K) = 1./TEMP4
      NSM2 = NS-2
      if (NSM2.lt.1) go to 210
      do 200 J = 1, NSM2
         I = K-J
         V(I) = V(I)-ALPHA(J+1)*V(I+1)
200   continue
210   continue
      LIMIT1 = KP1-NS
      TEMP5 = ALPHA(NS)
      do 220 IQ = 1, LIMIT1
         V(IQ) = V(IQ)-TEMP5*V(IQ+1)
         W(IQ) = V(IQ)
220   continue
      G(NSP1) = W(1)
230   continue
      NSP2 = NS+2
      if (KP1.lt.NSP2) go to 260
      do 250 I = NSP2, KP1
         LIMIT2 = KP2-I
         TEMP6 = ALPHA(I-1)
         do 240 IQ = 1, LIMIT2
            W(IQ) = W(IQ)-TEMP6*W(IQ+1)
240      continue
         G(I) = W(1)
250   continue
260   continue
      if (K.lt.NSP1) go to 290
      do 280 I = NSP1, K
         TEMP1 = BETA(I)
         do 270 L = 1, NEQN
            PHI(L,I) = TEMP1*PHI(L,I)
270      continue
280   continue
290   continue
      do 300 L = 1, NEQN
         PHI(L,KP2) = PHI(L,KP1)
         PHI(L,KP1) = 0.
         P(L) = 0.
300   continue
      do 320 J = 1, K
         I = KP1-J
         IP1 = I+1
         TEMP2 = G(I)
         do 310 L = 1, NEQN
            P(L) = P(L)+TEMP2*PHI(L,I)
            PHI(L,I) = PHI(L,I)+PHI(L,IP1)
310      continue
320   continue
      if (NORND) go to 340
      do 330 L = 1, NEQN
         TAU = H*P(L)-PHI(L,15)
         P(L) = Y(L)+TAU
         PHI(L,16) = (P(L)-Y(L))-TAU
330   continue
      go to 360
340   continue
      do 350 L = 1, NEQN
         P(L) = Y(L)+H*P(L)
350   continue
360   continue
      XOLD = X
      X = X+H
      ABSH = dabs(H)
      call FCT (X,P,YP)
      NFE = NFE+1
      ERKM2 = 0.
      ERKM1 = 0.
      ERK = 0.
      do 390 L = 1, NEQN
         TEMP3 = 1./WT(L)
         TEMP4 = YP(L)-PHI(L,1)
         if (KM2) 390, 380, 370
370      continue
         ERKM2 = dmax1(ERKM2,dabs((PHI(L,KM1)+TEMP4)*TEMP3))
380      continue
         ERKM1 = dmax1(ERKM1,dabs((PHI(L,K)+TEMP4)*TEMP3))
         ERK = dmax1(ERK,dabs(TEMP4*TEMP3))
390   continue
      if (KM2) 420, 410, 400
400   continue
      ERKM2 = ABSH*SIG(KM1)*GSTR(KM2)*ERKM2
410   continue
      ERKM1 = ABSH*SIG(K)*GSTR(KM1)*ERKM1
420   continue
      TEMP5 = ABSH*ERK
      err = TEMP5*(G(K)-G(KP1))
      ERK = TEMP5*SIG(KP1)*GSTR(K)
      KNEW = K
      if (KM2) 450, 440, 430
430   continue
      if (dmax1(ERKM1,ERKM2).le.ERK) KNEW = KM1
      go to 450
440   continue
      if (ERKM1.le.0.5*ERK) KNEW = KM1
450   continue
      if (err.le.EPS) go to 530
      PHASE1 = .false.
      X = XOLD
      do 470 I = 1, K
         TEMP1 = 1./BETA(I)
         IP1 = I+1
         do 460 L = 1, NEQN
            PHI(L,I) = TEMP1*(PHI(L,I)-PHI(L,IP1))
460      continue
470   continue
      if (K.lt.2) go to 490
      do 480 I = 2, K
         PSI(I-1) = PSI(I)-H
480   continue
490   continue
      IFAIL = IFAIL+1
      TEMP2 = 0.5
      if (IFAIL-3) 520, 510, 500
500   continue
      if (EPS.lt.0.5*ERK) TEMP2 = dsqrt(0.5*EPS/ERK)
510   continue
      KNEW = 1
520   continue
      H = TEMP2*H
      K = KNEW
      if (dabs(H).ge.FOURU*dabs(X)) go to 150
      ICRASH = 3
      H = dsign(FOURU*dabs(X),H)
      EPS = EPS+EPS
      return
530   continue
      KOLD = K
      HOLD = H
      TEMP1 = H*G(KP1)
      if (NORND) go to 550
      do 540 L = 1, NEQN
         RHO = TEMP1*(YP(L)-PHI(L,1))-PHI(L,16)
         Y(L) = P(L)+RHO
         PHI(L,15) = (Y(L)-P(L))-RHO
540   continue
      go to 570
550   continue
      do 560 L = 1, NEQN
         Y(L) = P(L)+TEMP1*(YP(L)-PHI(L,1))
560   continue
570   continue
      call FCT (X,Y,YP)
      NFE = NFE+1
      do 580 L = 1, NEQN
         PHI(L,KP1) = YP(L)-PHI(L,1)
         PHI(L,KP2) = PHI(L,KP1)-PHI(L,KP2)
580   continue
      do 590 I = 1, K
         do 590 L = 1, NEQN
            PHI(L,I) = PHI(L,I)+PHI(L,KP1)
590   continue
      ERKP1 = 0.
      if ((KNEW.eq.KM1).or.(K.eq.12)) PHASE1 = .false.
      if (PHASE1) go to 620
      if (KNEW.eq.KM1) go to 630
      if (KP1.gt.NS) go to 640
      do 600 L = 1, NEQN
         ERKP1 = dmax1(ERKP1,dabs(PHI(L,KP2)/WT(L)))
600   continue
      ERKP1 = ABSH*GSTR(KP1)*ERKP1
      if (K.gt.1) go to 610
      if (ERKP1.ge.0.5*ERK) go to 640
      go to 620
610   continue
      if (ERKM1.le.dmin1(ERK,ERKP1)) go to 630
      if ((ERKP1.ge.ERK).or.(K.eq.12)) go to 640
620   continue
      K = KP1
      ERK = ERKP1
      go to 640
630   continue
      K = KM1
      ERK = ERKM1
640   continue
      HNEW = H+H
      if (PHASE1) go to 650
      if (EPS.ge.ERK*TWO(K)) go to 650
      HNEW = H
      if (EPS.ge.2.*ERK) go to 650
      R = (0.5*EPS/ERK)**(1./float(K+1))
      if (R.gt.0.9) R = 0.9
      if (R.lt.0.5) R = 0.5
      HNEW = dsign(dmax1(ABSH*R,FOURU*dabs(X)),H)
650   continue
      H = HNEW
      return
      end
