        PROGRAM D4R2
C       Driver for routine QTRAP
        PARAMETER(PIO2=1.5707963)
        A=0.0
        B=PIO2
        WRITE(*,'(1X,A)') 'Integral of FUNC computed with QTRAP'
        WRITE(*,'(1X,A,F10.6)') 'Actual value of integral is',
     *          FINT(B)-FINT(A)
        CALL QTRAP(A,B,S)
        WRITE(*,'(1X,A,F10.6)') 'Result from routine QTRAP is',S
        END
C
        function FUNC(X)
        FUNC=(X**2)*(X**2-2.0)*SIN(X)
        END
C
        FUNCTION FINT(X)
C       Integral of FUNC
        FINT=4.0*X*((X**2)-7.0)*SIN(X)
     *          -((X**4)-14.0*(X**2)+28.0)*COS(X)
        END
C
      SUBROUTINE QTRAP(A,B,S)
      PARAMETER (EPS=0.5e-5, JMAX=100)
      OLDS=-1.E30
      DO 11 J=1,JMAX
        CALL TRAPZD(A,B,S,J)
        IF (ABS(S-OLDS).LT.EPS) RETURN
        OLDS=S
11    CONTINUE
      PAUSE 'Too many steps.'
      END
C
      SUBROUTINE TRAPZD(A,B,S,N)
      IF (N.EQ.1) THEN
        S=0.5*(B-A)*(FUNC(a)+FUNC(b))
      ELSE
        DEL=(B-A)/n
        SUM=0.5*(func(a)+func(b))
        x=a+del
        DO 11 J=2,n
          SUM=SUM+FUNC(x)
          X=X+DEL
11      CONTINUE
        S=del*sum
      ENDIF
      RETURN
      END
