      SUBROUTINE STEP1(F,NEQN,Y,X,H,EPS,WT,START,
     1 HOLD,K,KOLD,CRASH,PHI,P,YP,PSI,
     2 ALPHA,BETA,SIG,V,W,G,PHASE1,NS,NORND)
C
C     SANDIA MATHEMATICAL PROGRAM LIBRARY
C     APPLIED MATHEMATICS DIVISION 2642
C     SANDIA LABORATORIES
C     ALBUQUERQUE, NEW MEXICO  87115
C     CONTROL DATA 6600 VERSION 6.1  JANUARY 1976
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C   *                 ISSUED BY SANDIA LABORATORIES,
C   *                   A PRIME CONTRACTOR TO THE
C   *   UNITED STATES ENERGY RESEARCH AND DEVELOPMENT ADMINISTRATION
C   * * * * * * * * * * * * * *  NOTICE   * * * * * * * * * * * * * * *
C   * THIS REPORT WAS PREPARED AS AN ACCOUNT OF WORK SPONSORED BY THE
C   * UNITED STATES GOVERNMENT.  NEITHER THE UNITED STATES NOR THE
C   * UNITED STATES ENERGY RESEARCH AND DEVELOPMENT ADMINISTRATION,
C   * NOR ANY OF THEIR EMPLOYEES, NOR ANY OF THEIR CONTRACTORS,
C   * SUBCONTRACTORS, OR THEIR EMPLOYEES, MAKES ANY WARRANTY, EXPRESS
C   * OR IMPLIED, OR ASSUMES ANY LEGAL LIABILITY OR RESPONSIBILITY
C   * FOR THE ACCURACY, COMPLETENESS OR USEFULNESS OF ANY INFORMATION,
C   * APPARATUS, PRODUCT OR PROCESS DISCLOSED, OR REPRESENTS THAT ITS
C   * USE WOULD NOT INFRINGE PRIVATELY OWNED RIGHTS.
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C   * THE PRIMARY DOCUMENT FOR THE LIBRARY OF WHICH THIS ROUTINE IS
C   * IS A PART IS SAND75-0545, PRINTED FEB 1976.
C   * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
C
C   WRITTEN BY L. F. SHAMPINE AND M. K. GORDON
C
C   ABSTRACT
C
C
C   SUBROUTINE  STEP1  IS NORMALLY USED INDIRECTLY THROUGH SUBROUTINE
C   ODE .  BECAUSE  ODE  SUFFICES FOR MOST PROBLEMS AND IS MUCH EASIER
C   TO USE, USING IT SHOULD BE CONSIDERED BEFORE USING  STEP1  ALONE.
C
C   SUBROUTINE  STEP1  INTEGRATES A SYSTEM OF  NEQN  FIRST ORDER ORDINAR
C   DIFFERENTIAL EQUATIONS ONE STEP, NORMALLY FROM X TO X+H, USING A
C   MODIFIED DIVIDED DIFFERENCE FORM OF THE ADAMS PECE FORMULAS.  LOCAL
C   EXTRAPOLATION IS USED TO IMPROVE ABSOLUTE STABILITY AND ACCURACY.
C   THE CODE ADJUSTS ITS ORDER AND STEP SIZE TO CONTROL THE LOCAL ERROR
C   PER UNIT STEP IN A GENERALIZED SENSE.  SPECIAL DEVICES ARE INCLUDED
C   TO CONTROL ROUNDOFF ERROR AND TO DETECT WHEN THE USER IS REQUESTING
C   TOO MUCH ACCURACY.
C
C   THIS CODE IS COMPLETELY EXPLAINED AND DOCUMENTED IN THE TEXT,
C   COMPUTER SOLUTION OF ORDINARY DIFFERENTIAL EQUATIONS0  THE INITIAL
C   VALUE PROBLEM  BY L. F. SHAMPINE AND M. K. GORDON.
C
C
C   THE PARAMETERS REPRESENT0
C      F -- SUBROUTINE TO EVALUATE DERIVATIVES
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED
C      Y(*) -- SOLUTION VECTOR AT X
C      X -- INDEPENDENT VARIABLE
C      H -- APPROPRIATE STEP SIZE FOR NEXT STEP.  NORMALLY DETERMINED BY
C           CODE
C      EPS -- LOCAL ERROR TOLERANCE
C      WT(*) -- VECTOR OF WEIGHTS FOR ERROR CRITERION
C      START -- LOGICAL VARIABLE SET .TRUE. FOR FIRST STEP,  .FALSE.
C           OTHERWISE
C      HOLD -- STEP SIZE USED FOR LAST SUCCESSFUL STEP
C      K -- APPROPRIATE ORDER FOR NEXT STEP (DETERMINED BY CODE)
C      KOLD -- ORDER USED FOR LAST SUCCESSFUL STEP
C      CRASH -- LOGICAL VARIABLE SET .TRUE. WHEN NO STEP CAN BE TAKEN,
C           .FALSE. OTHERWISE.
C      YP(*) -- DERIVATIVE OF SOLUTION VECTOR AT  X  AFTER SUCCESSFUL
C           STEP
C   THE ARRAYS  PHI, PSI  ARE REQUIRED FOR THE INTERPOLATION SUBROUTINE
C   INTRP .  THE ARRAY  P  IS INTERNAL TO THE CODE.  THE REMAINING NINE
C   VARIABLES AND ARRAYS ARE INCLUDED IN THE CALL LIST ONLY TO ELIMINATE
C   LOCAL RETENTION OF VARIABLES BETWEEN CALLS.
C
C   INPUT TO STEP1
C
C      FIRST CALL --
C
C   THE USER MUST PROVIDE STORAGE IN HIS CALLING PROGRAM FOR ALL ARRAYS
C   IN THE CALL LIST, NAMELY
C
C     DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12),
C    1  ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13)
C
C   THE USER MUST ALSO DECLARE  START ,  CRASH ,  PHASE1  AND  NORND
C   LOGICAL VARIABLES AND  F  AN EXTERNAL SUBROUTINE, SUPPLY THE
C   SUBROUTINE  F(X,Y,YP)  TO EVALUATE
C      DY(I)/DX = YP(I) = F(X,Y(1),Y(2),...,Y(NEQN))
C   AND INITIALIZE ONLY THE FOLLOWING PARAMETERS0
C      NEQN -- NUMBER OF EQUATIONS TO BE INTEGRATED
C      Y(*) -- VECTOR OF INITIAL VALUES OF DEPENDENT VARIABLES
C      X -- INITIAL VALUE OF THE INDEPENDENT VARIABLE
C      H -- NOMINAL STEP SIZE INDICATING DIRECTION OF INTEGRATION
C           AND MAXIMUM SIZE OF STEP.  MUST BE VARIABLE
C      EPS -- LOCAL ERROR TOLERANCE PER STEP.  MUST BE VARIABLE
C      WT(*) -- VECTOR OF NON-ZERO WEIGHTS FOR ERROR CRITERION
C      START -- .TRUE.
C
C   STEP1  REQUIRES THAT THE L2 NORM OF THE VECTOR WITH COMPONENTS
C   LOCAL ERROR(L)/WT(L)  BE LESS THAN  EPS  FOR A SUCCESSFUL STEP.  THE
C   ARRAY  WT  ALLOWS THE USER TO SPECIFY AN ERROR TEST APPROPRIATE
C   FOR HIS PROBLEM.  FOR EXAMPLE,
C      WT(L) = 1.0  SPECIFIES ABSOLUTE ERROR,
C            = ABS(Y(L))  ERROR RELATIVE TO THE MOST RECENT VALUE OF THE
C                 L-TH COMPONENT OF THE SOLUTION,
C            = ABS(YP(L))  ERROR RELATIVE TO THE MOST RECENT VALUE OF
C                 THE L-TH COMPONENT OF THE DERIVATIVE,
C            = AMAX1(WT(L),ABS(Y(L)))  ERROR RELATIVE TO THE LARGEST
C                 MAGNITUDE OF L-TH COMPONENT OBTAINED SO FAR,
C            = ABS(Y(L))*RELERR/EPS + ABSERR/EPS  SPECIFIES A MIXED
C                 RELATIVE-ABSOLUTE TEST WHERE  RELERR  IS RELATIVE
C                 ERROR,  ABSERR  IS ABSOLUTE ERROR AND  EPS =
C                 AMAX1(RELERR,ABSERR) .
C
C      SUBSEQUENT CALLS --
C
C   SUBROUTINE  STEP1  IS DESIGNED SO THAT ALL INFORMATION NEEDED TO
C   CONTINUE THE INTEGRATION, INCLUDING THE STEP SIZE  H  AND THE ORDER
C   K , IS RETURNED WITH EACH STEP.  WITH THE EXCEPTION OF THE STEP
C   SIZE, THE ERROR TOLERANCE, AND THE WEIGHTS, NONE OF THE PARAMETERS
C   SHOULD BE ALTERED.  THE ARRAY  WT  MUST BE UPDATED AFTER EACH STEP
C   TO MAINTAIN RELATIVE ERROR TESTS LIKE THOSE ABOVE.  NORMALLY THE
C   INTEGRATION IS CONTINUED JUST BEYOND THE DESIRED ENDPOINT AND THE
C   SOLUTION INTERPOLATED THERE WITH SUBROUTINE  INTRP .  IF IT IS
C   IMPOSSIBLE TO INTEGRATE BEYOND THE ENDPOINT, THE STEP SIZE MAY BE
C   REDUCED TO HIT THE ENDPOINT SINCE THE CODE WILL NOT TAKE A STEP
C   LARGER THAN THE  H  INPUT.  CHANGING THE DIRECTION OF INTEGRATION,
C   I.E., THE SIGN OF  H , REQUIRES THE USER SET  START = .TRUE. BEFORE
C   CALLING  STEP1  AGAIN.  THIS IS THE ONLY SITUATION IN WHICH  START
C   SHOULD BE ALTERED.
C
C   OUTPUT FROM STEP1
C
C      SUCCESSFUL STEP --
C
C   THE SUBROUTINE RETURNS AFTER EACH SUCCESSFUL STEP WITH  START  AND
C   CRASH  SET .FALSE. .  X  REPRESENTS THE INDEPENDENT VARIABLE
C   ADVANCED ONE STEP OF LENGTH  HOLD  FROM ITS VALUE ON INPUT AND  Y
C   THE SOLUTION VECTOR AT THE NEW VALUE OF  X .  ALL OTHER PARAMETERS
C   REPRESENT INFORMATION CORRESPONDING TO THE NEW  X  NEEDED TO
C   CONTINUE THE INTEGRATION.
C
C      UNSUCCESSFUL STEP --
C
C   WHEN THE ERROR TOLERANCE IS TOO SMALL FOR THE MACHINE PRECISION,
C   THE SUBROUTINE RETURNS WITHOUT TAKING A STEP AND  CRASH = .TRUE. .
C   AN APPROPRIATE STEP SIZE AND ERROR TOLERANCE FOR CONTINUING ARE
C   ESTIMATED AND ALL OTHER INFORMATION IS RESTORED AS UPON INPUT
C   BEFORE RETURNING.  TO CONTINUE WITH THE LARGER TOLERANCE, THE USER
C   JUST CALLS THE CODE AGAIN.  A RESTART IS NEITHER REQUIRED NOR
C   DESIRABLE.
C
      LOGICAL START,CRASH,PHASE1,NORND
      DIMENSION Y(NEQN),WT(NEQN),PHI(NEQN,16),P(NEQN),YP(NEQN),PSI(12),
     1  ALPHA(12),BETA(12),SIG(13),V(12),W(12),G(13)
      DIMENSION TWO(13),GSTR(13)
      EXTERNAL F
C
      DATA TWO(1)/2.0/, TWO(2)/4.0/, TWO(3)/8.0/, TWO(4)/16.0/,
     1     TWO(5)/32.0/, TWO(6)/64.0/, TWO(7)/128.0/, TWO(8)/256.0/,
     2     TWO(9)/512.0/, TWO(10)/1024.0/, TWO(11)/2048.0/,
     3     TWO(12)/4096.0/, TWO(13)/8192.0/
      DATA GSTR(1)/0.500/, GSTR(2)/0.0833/, GSTR(3)/0.0417/,
     1     GSTR(4)/0.0264/, GSTR(5)/0.0188/, GSTR(6)/0.0143/,
     2     GSTR(7)/0.0114/, GSTR(8)/0.00936/, GSTR(9)/0.00789/,
     3     GSTR(10)/0.00679/, GSTR(11)/0.00592/, GSTR(12)/0.00524/,
     4     GSTR(13)/0.00468/
C
C   ****** U IS A MACHINE DEPENDENT PARAMETER. IT IS THE SMALLEST
C          POSITIVE NUMBER FOR WHICH 1.0 + U .GT. 1.0.
C
      U = SPMPAR(1)
      TWOU = 2.0*U
      FOURU = 4.0*U
C
C
C       ***     BEGIN BLOCK 0     ***
C   CHECK IF STEP SIZE OR ERROR TOLERANCE IS TOO SMALL FOR MACHINE
C   PRECISION.  IF FIRST STEP, INITIALIZE PHI ARRAY AND ESTIMATE A
C   STARTING STEP SIZE.
C                   ***
C
C   IF STEP SIZE IS TOO SMALL, DETERMINE AN ACCEPTABLE ONE
C
      CRASH = .TRUE.
      IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 5
      H = SIGN(FOURU*ABS(X),H)
      RETURN
 5    P5EPS = 0.5*EPS
C
C   IF ERROR TOLERANCE IS TOO SMALL, INCREASE IT TO AN ACCEPTABLE VALUE
C
      ROUND = 0.0
      DO 10 L = 1,NEQN
 10     ROUND = ROUND + (Y(L)/WT(L))**2
      ROUND = TWOU*SQRT(ROUND)
      IF(P5EPS .GE. ROUND) GO TO 15
      EPS = 2.0*ROUND*(1.0 + FOURU)
      RETURN
 15   CRASH = .FALSE.
      G(1) = 1.0
      G(2) = 0.5
      SIG(1) = 1.0
      IF(.NOT.START) GO TO 99
C
C   INITIALIZE.  COMPUTE APPROPRIATE STEP SIZE FOR FIRST STEP
C
      CALL F(X,Y,YP)
      SUM = 0.0
      DO 20 L = 1,NEQN
        PHI(L,1) = YP(L)
        PHI(L,2) = 0.0
 20     SUM = SUM + (YP(L)/WT(L))**2
      SUM = SQRT(SUM)
      ABSH = ABS(H)
      IF(EPS .LT. 16.0*SUM*H*H) ABSH = 0.25*SQRT(EPS/SUM)
      H = SIGN(AMAX1(ABSH,FOURU*ABS(X)),H)
      HOLD = 0.0
      K = 1
      KOLD = 0
      START = .FALSE.
      PHASE1 = .TRUE.
      NORND = .TRUE.
      IF(P5EPS .GT. 100.0*ROUND) GO TO 99
      NORND = .FALSE.
      DO 25 L = 1,NEQN
 25     PHI(L,15) = 0.0
 99   IFAIL = 0
C       ***     END BLOCK 0     ***
C
C       ***     BEGIN BLOCK 1     ***
C   COMPUTE COEFFICIENTS OF FORMULAS FOR THIS STEP.  AVOID COMPUTING
C   THOSE QUANTITIES NOT CHANGED WHEN STEP SIZE IS NOT CHANGED.
C                   ***
C
 100  KP1 = K+1
      KP2 = K+2
      KM1 = K-1
      KM2 = K-2
C
C   NS IS THE NUMBER OF STEPS TAKEN WITH SIZE H, INCLUDING THE CURRENT
C   ONE.  WHEN K.LT.NS, NO COEFFICIENTS CHANGE
C
      IF(H .NE. HOLD) NS = 0
      IF (NS.LE.KOLD) NS = NS+1
      NSP1 = NS+1
      IF (K .LT. NS) GO TO 199
C
C   COMPUTE THOSE COMPONENTS OF ALPHA(*),BETA(*),PSI(*),SIG(*) WHICH
C   ARE CHANGED
C
      BETA(NS) = 1.0
      REALNS = NS
      ALPHA(NS) = 1.0/REALNS
      TEMP1 = H*REALNS
      SIG(NSP1) = 1.0
      IF(K .LT. NSP1) GO TO 110
      DO 105 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
        REALI = I
 105    SIG(I+1) = REALI*ALPHA(I)*SIG(I)
 110  PSI(K) = TEMP1
C
C   COMPUTE COEFFICIENTS G(*)
C
C   INITIALIZE V(*) AND SET W(*).
C
      IF(NS .GT. 1) GO TO 120
      DO 115 IQ = 1,K
        TEMP3 = IQ*(IQ+1)
        V(IQ) = 1.0/TEMP3
 115    W(IQ) = V(IQ)
      GO TO 140
C
C   IF ORDER WAS RAISED, UPDATE DIAGONAL PART OF V(*)
C
 120  IF(K .LE. KOLD) GO TO 130
      TEMP4 = K*KP1
      V(K) = 1.0/TEMP4
      NSM2 = NS-2
      IF(NSM2 .LT. 1) GO TO 130
      DO 125 J = 1,NSM2
        I = K-J
 125    V(I) = V(I) - ALPHA(J+1)*V(I+1)
C
C   UPDATE V(*) AND SET W(*)
C
 130  LIMIT1 = KP1 - NS
      TEMP5 = ALPHA(NS)
      DO 135 IQ = 1,LIMIT1
        V(IQ) = V(IQ) - TEMP5*V(IQ+1)
 135    W(IQ) = V(IQ)
      G(NSP1) = W(1)
C
C   COMPUTE THE G(*) IN THE WORK VECTOR W(*)
C
 140  NSP2 = NS + 2
      IF(KP1 .LT. NSP2) GO TO 199
      DO 150 I = NSP2,KP1
        LIMIT2 = KP2 - I
        TEMP6 = ALPHA(I-1)
        DO 145 IQ = 1,LIMIT2
 145      W(IQ) = W(IQ) - TEMP6*W(IQ+1)
 150    G(I) = W(1)
 199    CONTINUE
C       ***     END BLOCK 1     ***
C
C       ***     BEGIN BLOCK 2     ***
C   PREDICT A SOLUTION P(*), EVALUATE DERIVATIVES USING PREDICTED
C   SOLUTION, ESTIMATE LOCAL ERROR AT ORDER K AND ERRORS AT ORDERS K,
C   K-1, K-2 AS IF CONSTANT STEP SIZE WERE USED.
C                   ***
C
C   CHANGE PHI TO PHI STAR
C
      IF(K .LT. NSP1) GO TO 215
      DO 210 I = NSP1,K
        TEMP1 = BETA(I)
        DO 205 L = 1,NEQN
 205      PHI(L,I) = TEMP1*PHI(L,I)
 210    CONTINUE
C
C   PREDICT SOLUTION AND DIFFERENCES
C
 215  DO 220 L = 1,NEQN
        PHI(L,KP2) = PHI(L,KP1)
        PHI(L,KP1) = 0.0
 220    P(L) = 0.0
      DO 230 J = 1,K
        I = KP1 - J
        IP1 = I+1
        TEMP2 = G(I)
        DO 225 L = 1,NEQN
          P(L) = P(L) + TEMP2*PHI(L,I)
 225      PHI(L,I) = PHI(L,I) + PHI(L,IP1)
 230    CONTINUE
      IF(NORND) GO TO 240
      DO 235 L = 1,NEQN
        TAU = H*P(L) - PHI(L,15)
        P(L) = Y(L) + TAU
 235    PHI(L,16) = (P(L) - Y(L)) - TAU
      GO TO 250
 240  DO 245 L = 1,NEQN
 245    P(L) = Y(L) + H*P(L)
 250  XOLD = X
      X = X + H
      ABSH = ABS(H)
      CALL F(X,P,YP)
C
C   ESTIMATE ERRORS AT ORDERS K,K-1,K-2
C
      ERKM2 = 0.0
      ERKM1 = 0.0
      ERK = 0.0
      DO 265 L = 1,NEQN
        TEMP3 = 1.0/WT(L)
        TEMP4 = YP(L) - PHI(L,1)
        IF(KM2)265,260,255
 255    ERKM2 = ERKM2 + ((PHI(L,KM1)+TEMP4)*TEMP3)**2
 260    ERKM1 = ERKM1 + ((PHI(L,K)+TEMP4)*TEMP3)**2
 265    ERK = ERK + (TEMP4*TEMP3)**2
      IF(KM2)280,275,270
 270  ERKM2 = ABSH*SIG(KM1)*GSTR(KM2)*SQRT(ERKM2)
 275  ERKM1 = ABSH*SIG(K)*GSTR(KM1)*SQRT(ERKM1)
 280  TEMP5 = ABSH*SQRT(ERK)
      ERR = TEMP5*(G(K)-G(KP1))
      ERK = TEMP5*SIG(KP1)*GSTR(K)
      KNEW = K
C
C   TEST IF ORDER SHOULD BE LOWERED
C
      IF(KM2)299,290,285
 285  IF(AMAX1(ERKM1,ERKM2) .LE. ERK) KNEW = KM1
      GO TO 299
 290  IF(ERKM1 .LE. 0.5*ERK) KNEW = KM1
C
C   TEST IF STEP SUCCESSFUL
C
 299  IF(ERR .LE. EPS) GO TO 400
C       ***     END BLOCK 2     ***
C
C       ***     BEGIN BLOCK 3     ***
C   THE STEP IS UNSUCCESSFUL.  RESTORE  X, PHI(*,*), PSI(*) .
C   IF THIRD CONSECUTIVE FAILURE, SET ORDER TO ONE.  IF STEP FAILS MORE
C   THAN THREE TIMES, CONSIDER AN OPTIMAL STEP SIZE.  DOUBLE ERROR
C   TOLERANCE AND RETURN IF ESTIMATED STEP SIZE IS TOO SMALL FOR MACHINE
C   PRECISION.
C                   ***
C
C   RESTORE X, PHI(*,*) AND PSI(*)
C
      PHASE1 = .FALSE.
      X = XOLD
      DO 310 I = 1,K
        TEMP1 = 1.0/BETA(I)
        IP1 = I+1
        DO 305 L = 1,NEQN
 305      PHI(L,I) = TEMP1*(PHI(L,I) - PHI(L,IP1))
 310    CONTINUE
      IF(K .LT. 2) GO TO 320
      DO 315 I = 2,K
 315    PSI(I-1) = PSI(I) - H
C
C   ON THIRD FAILURE, SET ORDER TO ONE.  THEREAFTER, USE OPTIMAL STEP
C   SIZE
C
 320  IFAIL = IFAIL + 1
      TEMP2 = 0.5
      IF(IFAIL - 3) 335,330,325
 325  IF(P5EPS .LT. 0.25*ERK) TEMP2 = SQRT(P5EPS/ERK)
 330  KNEW = 1
 335  H = TEMP2*H
      K = KNEW
      IF(ABS(H) .GE. FOURU*ABS(X)) GO TO 340
      CRASH = .TRUE.
      H = SIGN(FOURU*ABS(X),H)
      EPS = EPS + EPS
      RETURN
 340  GO TO 100
C       ***     END BLOCK 3     ***
C
C       ***     BEGIN BLOCK 4     ***
C   THE STEP IS SUCCESSFUL.  CORRECT THE PREDICTED SOLUTION, EVALUATE
C   THE DERIVATIVES USING THE CORRECTED SOLUTION AND UPDATE THE
C   DIFFERENCES.  DETERMINE BEST ORDER AND STEP SIZE FOR NEXT STEP.
C                   ***
 400  KOLD = K
      HOLD = H
C
C   CORRECT AND EVALUATE
C
      TEMP1 = H*G(KP1)
      IF(NORND) GO TO 410
      DO 405 L = 1,NEQN
        RHO = TEMP1*(YP(L) - PHI(L,1)) - PHI(L,16)
        Y(L) = P(L) + RHO
 405    PHI(L,15) = (Y(L) - P(L)) - RHO
      GO TO 420
 410  DO 415 L = 1,NEQN
 415    Y(L) = P(L) + TEMP1*(YP(L) - PHI(L,1))
 420  CALL F(X,Y,YP)
C
C   UPDATE DIFFERENCES FOR NEXT STEP
C
      DO 425 L = 1,NEQN
        PHI(L,KP1) = YP(L) - PHI(L,1)
 425    PHI(L,KP2) = PHI(L,KP1) - PHI(L,KP2)
      DO 435 I = 1,K
        DO 430 L = 1,NEQN
 430      PHI(L,I) = PHI(L,I) + PHI(L,KP1)
 435    CONTINUE
C
C   ESTIMATE ERROR AT ORDER K+1 UNLESS0
C     IN FIRST PHASE WHEN ALWAYS RAISE ORDER,
C     ALREADY DECIDED TO LOWER ORDER,
C     STEP SIZE NOT CONSTANT SO ESTIMATE UNRELIABLE
C
      ERKP1 = 0.0
      IF(KNEW .EQ. KM1  .OR.  K .EQ. 12) PHASE1 = .FALSE.
      IF(PHASE1) GO TO 450
      IF(KNEW .EQ. KM1) GO TO 455
      IF(KP1 .GT. NS) GO TO 460
      DO 440 L = 1,NEQN
 440    ERKP1 = ERKP1 + (PHI(L,KP2)/WT(L))**2
      ERKP1 = ABSH*GSTR(KP1)*SQRT(ERKP1)
C
C   USING ESTIMATED ERROR AT ORDER K+1, DETERMINE APPROPRIATE ORDER
C   FOR NEXT STEP
C
      IF(K .GT. 1) GO TO 445
      IF(ERKP1 .GE. 0.5*ERK) GO TO 460
      GO TO 450
 445  IF(ERKM1 .LE. AMIN1(ERK,ERKP1)) GO TO 455
      IF(ERKP1 .GE. ERK  .OR.  K .EQ. 12) GO TO 460
C
C   HERE ERKP1 .LT. ERK .LT. AMAX1(ERKM1,ERKM2) ELSE ORDER WOULD HAVE
C   BEEN LOWERED IN BLOCK 2.  THUS ORDER IS TO BE RAISED
C
C   RAISE ORDER
C
 450  K = KP1
      ERK = ERKP1
      GO TO 460
C
C   LOWER ORDER
C
 455  K = KM1
      ERK = ERKM1
C
C   WITH NEW ORDER DETERMINE APPROPRIATE STEP SIZE FOR NEXT STEP
C
 460  HNEW = H + H
      IF(PHASE1) GO TO 465
      IF(P5EPS .GE. ERK*TWO(K+1)) GO TO 465
      HNEW = H
      IF(P5EPS .GE. ERK) GO TO 465
      TEMP2 = K+1
      R = (P5EPS/ERK)**(1.0/TEMP2)
      HNEW = ABSH*AMAX1(0.5,AMIN1(0.9,R))
      HNEW = SIGN(AMAX1(HNEW,FOURU*ABS(X)),H)
 465  H = HNEW
      RETURN
C       ***     END BLOCK 4     ***
      END
