      SUBROUTINE HSTART (F, NEQ, A, B, Y, YPRIME, ETOL, MORDER, SMALL,
     +   BIG, SPY, PV, YP, SF, RPAR, IPAR, H)
C
      DIMENSION Y(*),YPRIME(*),ETOL(*),SPY(*),PV(*),YP(*),SF(*),
     1   RPAR(*),IPAR(*)
      EXTERNAL F
C
C.......................................................................
C
C***FIRST EXECUTABLE STATEMENT  HSTART
      DX = B - A
      ABSDX = ABS(DX)
      RELPER = SMALL**0.375
      YNORM = HVNRM(Y,NEQ)
C
C.......................................................................
C
C     COMPUTE A WEIGHTED APPROXIMATE BOUND (DFDXB) ON THE PARTIAL
C     DERIVATIVE OF THE EQUATION WITH RESPECT TO THE
C     INDEPENDENT VARIABLE. PROTECT AGAINST AN OVERFLOW. ALSO
C     COMPUTE A WEIGHTED BOUND (FBND) ON THE FIRST DERIVATIVE LOCALLY.
C
      DA = SIGN(MAX(MIN(RELPER*ABS(A),ABSDX),100.*SMALL*ABS(A)),DX)
      IF (DA .EQ. 0.) DA = RELPER*DX
      CALL F(A+DA,Y,SF,RPAR,IPAR)
C
      IF (MORDER .EQ. 1) GO TO 20
      POWER = 2./(MORDER+1)
      DO 10 J=1,NEQ
        WTJ = ETOL(J)**POWER
        SPY(J) = SF(J)/WTJ
        YP(J) = YPRIME(J)/WTJ
   10   PV(J) = SPY(J) - YP(J)
      GO TO 40
C
   20 DO 30 J=1,NEQ
        SPY(J) = SF(J)/ETOL(J)
        YP(J) = YPRIME(J)/ETOL(J)
   30   PV(J) = SPY(J) - YP(J)
C
   40 DELF = HVNRM(PV,NEQ)
      DFDXB = BIG
      IF (DELF .LT. BIG*ABS(DA)) DFDXB = DELF/ABS(DA)
      YPNORM = HVNRM(YP,NEQ)
      FBND = MAX(HVNRM(SPY,NEQ),YPNORM)
C
C.......................................................................
C
C     COMPUTE AN ESTIMATE (DFDUB) OF THE LOCAL LIPSCHITZ CONSTANT FOR
C     THE SYSTEM OF DIFFERENTIAL EQUATIONS. THIS ALSO REPRESENTS AN
C     ESTIMATE OF THE NORM OF THE JACOBIAN LOCALLY.
C     THREE ITERATIONS (TWO WHEN NEQ=1) ARE USED TO ESTIMATE THE
C     LIPSCHITZ CONSTANT BY NUMERICAL DIFFERENCES. THE FIRST
C     PERTURBATION VECTOR IS BASED ON THE INITIAL DERIVATIVES AND
C     DIRECTION OF INTEGRATION. THE SECOND PERTURBATION VECTOR IS
C     FORMED USING ANOTHER EVALUATION OF THE DIFFERENTIAL EQUATION.
C     THE THIRD PERTURBATION VECTOR IS FORMED USING PERTURBATIONS BASED
C     ONLY ON THE INITIAL VALUES. COMPONENTS THAT ARE ZERO ARE ALWAYS
C     CHANGED TO NON-ZERO VALUES (EXCEPT ON THE FIRST ITERATION). WHEN
C     INFORMATION IS AVAILABLE, CARE IS TAKEN TO ENSURE THAT COMPONENTS
C     OF THE PERTURBATION VECTOR HAVE SIGNS WHICH ARE CONSISTENT WITH
C     THE SLOPES OF LOCAL SOLUTION CURVES.
C     ALSO CHOOSE THE LARGEST BOUND (FBND) FOR THE FIRST DERIVATIVE.
C     NO ATTEMPT IS MADE TO KEEP THE PERTURBATION VECTOR SIZE CONSTANT.
C
      IF (YPNORM .EQ. 0.) GO TO 60
C                       USE INITIAL DERIVATIVES FOR FIRST PERTURBATION
      ICASE = 1
      DO 50 J=1,NEQ
        SPY(J) = YPRIME(J)
   50   YP(J) = YPRIME(J)
      GO TO 80
C                       CANNOT HAVE A NULL PERTURBATION VECTOR
   60 ICASE = 2
      DO 70 J=1,NEQ
        SPY(J) = YPRIME(J)
   70   YP(J) = ETOL(J)
C
   80 DFDUB = 0.
      LK = MIN(NEQ+1,3)
      DO 260 K=1,LK
C                       SET YPNORM AND DELX
        YPNORM = HVNRM(YP,NEQ)
        IF (ICASE .EQ. 1  .OR.  ICASE .EQ. 3) GO TO 90
        DELX = SIGN(1.0,DX)
        GO TO 120
C                       TRY TO ENFORCE MEANINGFUL PERTURBATION VALUES
   90   DELX = DX
        IF (ABS(DELX)*YPNORM .GE. RELPER*YNORM) GO TO 100
        DELXB = BIG
        IF (RELPER*YNORM .LT. BIG*YPNORM) DELXB = RELPER*YNORM/YPNORM
        DELX = SIGN(DELXB,DX)
  100   DO 110 J=1,NEQ
          IF (ABS(DELX*YP(J)) .GT. ETOL(J)) DELX=SIGN(ETOL(J)/YP(J),DX)
  110     CONTINUE
C                       DEFINE PERTURBED VECTOR OF INITIAL VALUES
  120   DO 130 J=1,NEQ
  130     PV(J) = Y(J) + DELX*YP(J)
        IF (K .EQ. 2) GO TO 150
C                       EVALUATE DERIVATIVES ASSOCIATED WITH PERTURBED
C                       VECTOR  AND  COMPUTE CORRESPONDING DIFFERENCES
        CALL F(A,PV,YP,RPAR,IPAR)
        DO 140 J=1,NEQ
  140     PV(J) = YP(J) - YPRIME(J)
        GO TO 170
C                       USE A SHIFTED VALUE OF THE INDEPENDENT VARIABLE
C                                             IN COMPUTING ONE ESTIMATE
  150   CALL F(A+DA,PV,YP,RPAR,IPAR)
        DO 160 J=1,NEQ
  160     PV(J) = YP(J) - SF(J)
C                       CHOOSE LARGEST BOUND ON THE WEIGHTED FIRST
C                                                   DERIVATIVE
  170   IF (MORDER .EQ. 1) GO TO 190
        DO 180 J=1,NEQ
  180     YP(J) = YP(J)/ETOL(J)**POWER
        GO TO 210
  190   DO 200 J=1,NEQ
  200     YP(J) = YP(J)/ETOL(J)
  210   FBND = MAX(FBND,HVNRM(YP,NEQ))
C                       COMPUTE BOUND ON A LOCAL LIPSCHITZ CONSTANT
        DELF = HVNRM(PV,NEQ)
        IF (DELF .EQ. 0.) GO TO 220
        DELY = ABS(DELX)*YPNORM
        IF (DELF .GE. BIG*DELY) GO TO 270
        DFDUB = MAX(DFDUB,DELF/DELY)
C
  220   IF (K .EQ. LK) GO TO 280
C                       CHOOSE NEXT PERTURBATION VECTOR
        DO 250 J=1,NEQ
          IF (K .EQ. LK-1) GO TO 230
          ICASE = 3
          DY = ABS(PV(J))
          IF (DY .EQ. 0.) DY = MAX(DELF,ETOL(J))
          GO TO 240
  230     ICASE = 4
          DY = MAX(RELPER*ABS(Y(J)),ETOL(J))
  240     IF (SPY(J) .EQ. 0.) SPY(J) = YP(J)
          IF (SPY(J) .NE. 0.) DY = SIGN(DY,SPY(J))
  250     YP(J) = DY
  260   CONTINUE
C
C                       PROTECT AGAINST AN OVERFLOW
  270 DFDUB = BIG
C
C.......................................................................
C
C     COMPUTE A BOUND (YDPB) ON THE NORM OF THE SECOND DERIVATIVE
C
  280 YDPB = DFDXB + DFDUB*FBND
C
C.......................................................................
C
C     COMPUTE A STARTING STEP SIZE BASED ON THE ABOVE FIRST AND SECOND
C     DERIVATIVE INFORMATION
C
C                       RESTRICT THE STEP LENGTH TO BE NOT BIGGER THAN
C                       ABS(B-A).   (UNLESS  B  IS TOO CLOSE TO  A)
      H = ABSDX
C
      IF (YDPB .NE. 0.  .OR.  FBND .NE. 0.) GO TO 290
C
C                       BOTH FIRST DERIVATIVE TERM (FBND) AND SECOND
C                                    DERIVATIVE TERM (YDPB) ARE ZERO
      GO TO 310
C
  290 IF (YDPB .NE. 0.) GO TO 300
C
C                       ONLY SECOND DERIVATIVE TERM (YDPB) IS ZERO
      IF (1.0 .LT. FBND*ABSDX) H = 1./FBND
      GO TO 310
C
C                       SECOND DERIVATIVE TERM (YDPB) IS NON-ZERO
  300 SRYDPB = SQRT(0.5*YDPB)
      IF (1.0 .LT. SRYDPB*ABSDX) H = 1./SRYDPB
C
C                       FURTHER RESTRICT THE STEP LENGTH TO BE NOT
C                                                 BIGGER THAN  1/DFDUB
  310 IF (H*DFDUB .GT. 1.) H = 1./DFDUB
C
C                       FINALLY, RESTRICT THE STEP LENGTH TO BE NOT
C                       SMALLER THAN  100*SMALL*ABS(A).  HOWEVER, IF
C                       A=0. AND THE COMPUTED H UNDERFLOWED TO ZERO,
C                       THE ALGORITHM RETURNS  SMALL*ABS(B)  FOR THE
C                                                       STEP LENGTH.
      H = MAX(H,100.*SMALL*ABS(A))
      IF (H .EQ. 0.) H = SMALL*ABS(B)
C
C                       NOW SET DIRECTION OF INTEGRATION
      H = SIGN(H,DX)
C
      RETURN
      END
