
      FUNCTION AFunction (y)
      REAL y
        AFunction = EXP(y)
      END !FUNCTION



      SUBROUTINE IntegrateFunction (numseg, xl, xh, ans)
      INTEGER numseg
      REAL xl, xh, ans
      REAL strtpnt, endpnt, area1, area2, area, segwidth

      segwidth = (xh - xl) / (numseg * 1.0)

      endpnt = xh
      area1 = 0.0
      area2 = 0.0
      area = 0.0

      IF (MOD(numseg, 2) .NE. 0) THEN
      !!  SIMPSONS 3/8 rule
        area1 = 3.0 / 8.0 * segwidth *
     +          (AFunction(endpnt - 3 * segwidth) + 3 *
     +           AFunction(endpnt - 2 * segwidth) + 3 *
     +            AFunction(endpnt - segwidth) + AFunction(endpnt))
        endpnt = endpnt - 3 * segwidth
      ELSE
        area1 = 0.0
      END IF
      IF (numseg .NE. 3) THEN
      !!  SIMPSONS 1/3 rule
        strtpnt = xl
        IF  (strtpnt .GE. endpnt - segwidth) THEN
          area2 = area2 + 1.0 / 3.0 * (segwidth * (
     +            AFunction(strtpnt) + 4 *
     +            AFunction(strtpnt + segwidth) +
     +            AFunction(strtpnt + 2 * segwidth)))
        strtpnt = strtpnt + 2 * segwidth
        ELSE
          DO WHILE (strtpnt .LT. endpnt - segwidth)
            area2 = area2 + 1.0 / 3.0 * (segwidth * (
     +              AFunction(strtpnt) + 4 *
     +              AFunction(strtpnt + segwidth) +
     +              AFunction(strtpnt + 2 * segwidth)))
            strtpnt = strtpnt + 2 * segwidth
          END DO
        END IF
      END IF

      area = area1 + area2
      ans = area

      END !SUBROUTINE  Integrate Function




      SUBROUTINE IntegrateVector (y, samper, xl, xh, ans)
      INCLUDE 'STDHDR.FOR'
      REAL y(0:maxv), samper, ans
      INTEGER xl, xh

      REAL strtpnt, endpnt, area1, area2, area, segwidth


      numseg = xh - xl
      segwidth = samper
      endpnt = xh
      area1 = 0.0
      area2 = 0.0
      area = 0.0
      IF (MOD(numseg, 2) .NE. 0) THEN
      !! SIMPSONS 3/8 rule
        area1 = 3.0 / 8.0 * segwidth * (y(endpnt - 3) +
     +          3 * y(endpnt - 2) + 3 * y(endpnt - 1) + y(endpnt))
        endpnt = endpnt - 3
      ELSE
        area1 = 0.0
      END IF

      IF (numseg .NE. 3) THEN
      !!SIMPSONS 1/3 rule
        strtpnt = xl
        IF (strtpnt .GE. endpnt) THEN
          area2 = area2 + 1.0 / 3.0 * segwidth *
     +              (y(strtpnt) + 4 * y(strtpnt + 1) + y(strtpnt + 2))
            strtpnt = strtpnt + 2
        ELSE
          DO WHILE (strtpnt .LT. endpnt)
             area2 = area2 + 1.0 / 3.0 * segwidth *
     +              (y(strtpnt) + 4 * y(strtpnt + 1) + y(strtpnt + 2))
            strtpnt = strtpnt + 2
          END DO
        END IF
      END IF

      area = area1 + area2
      ans = area
      END !SUBROUTINE

      SUBROUTINE RungeKutta(XKnown, YKnown, x, et, h,
     +     HMin, YEst, status)
      REAL XKnown, YKnown, x, et, h, HMin, YEst
      INTEGER status

      PARAMETER ( NearZero = 1E-20)
      REAL  xc2,xc3,xc4,xc5,xc6,dc21,dc31,dc32,dc41,dc42,dc43
      REAL  dc51,dc52,dc53,dc54,dc61,dc62,dc63,dc64,dc65
      REAL  yc1,yc3,yc4,yc5,ec1,ec3,ec4,ec5,ec6,d1,d2,d3,d4,d5,d6
      REAL  errr,XEst,EtDiv32

      !Calculate the coefficients for the RKF equations

      xc2 = .25
      xc3 = 3.0 / 8.0
      xc4 = 12.0 / 13.0
      xc5 = 1.0
      xc6 = .5
      dc21 = .25
      dc31 = 3.0 / 32.0
      dc32 = 9.0 / 32.0
      dc41 = 1932.0 / 2197.0
      dc42 = 7200.0 / 2197.0
      dc43 = 7296.0 / 2197.0
      dc51 = 439.0 / 216.0
      dc52 = 8.0
      dc53 = 3680.0 / 513.0
      dc54 = 845.0 / 4104.0
      dc61 = 8.0/ 27.0
      dc62 = 2.0
      dc63 = 3544.0 / 2565.0
      dc64 = 1859.0 / 4104.0
      dc65 = 11.0 / 40.0
      yc1 = 25.0 / 216.0
      yc3 = 1408.0 / 2565.0
      yc4 = 2197.0 / 4104.0
      yc5 = 0.2
      ec1 = 1.0 / 360.0
      ec3 = 128.0 / 4275.0
      ec4 = 2197.0 / 75240.0
      ec5 = 0.02
      ec6 = 2.0 / 55.0

      EtDiv32 = et / 32.0
      XEst = XKnown
      YEst = YKnown
      st = 0
      DO WHILE (st .EQ. 0 )
        d1 = YPrime(XEst, YEst)
        d2 = YPrime(XEst + h * xc2, YEst + h * dc21 * d1)
        d3 = YPrime(XEst + h * xc3, YEst + h * (dc31 * d1 + dc32 * d2))
        d4 = YPrime(XEst + h * xc4, YEst + h *
     +     (dc41 * d1 - dc42 * d2 + dc43 * d3))
        d5 = YPrime(XEst + h * xc5, YEst + h *
     +     (dc51 * d1 - dc52 * d2 + dc53 * d3 - dc54 * d4))
        d6 = YPrime(XEst + h * xc6, YEst + h *
     +     (-dc61 * d1 + dc62 * d2 - dc63 * d3 + dc64 * d4 - dc65*d5))
        errr = h * (ec1 * d1 - ec3 * d3 - ec4 *
     +       d4 + ec5 * d5 + ec6 * d6)

        IF (errr .LT. (h * et)) THEN
          XEst = XEst + h
          YEst = YEst + h * (yc1*d1+yc3 * d3 + yc4 * d4 - yc5 * d5)
          IF (ABS(x - XEst) .GT. NearZero) THEN
            IF (errr .LE. (h * EtDiv32))  h = h * 2
            IF (h .GT. (x - XEst))   h = x - XEst
          ELSE
            st = 1
          END IF
        ELSE
          h = h / 2.0
          IF (h .LT. HMin)  st = 2
        END IF
        IF (st .EQ. 2) THEN
          status = 0
        ELSE
          status = 1
        END IF
      END DO !while loop

      END !SUBROUTINE {procedure}

      FUNCTION YPrime (x, y)
      REAL x, y
        YPrime = EXP(2.0 * x) + y
      END !FUNCTION

