The subroutines DINT_BODE and QINT_BODE calculate integrals of the one independent variable function by the Bode quadrature formula with a given order (from 2 up to 12) by steps h on a uniform grid.

Structure:

Type:	-	FUNCTION
User Entry Names:	-	DINT_BODE, QINT_BODE
External References:	-	FUNC - user-supplied function

Usage:

S=DINT_BODE(FUNC,A,B,NN,IPOINTS) double precision
S=QINT_BODE(FUNC,A,B,NN,IPOINTS) quadruple precision, where:

A,B	-	are limits of integration;
NN	-	is number of auxiliary subintervals;
IPOINTS	-	is number of nodes of the Bode quadrature formula; IPOINTS can accept values from 2 up to 11, for double precision calculation - up to 6;
FUNC	-	is name of the subprogram-function, constituent users for calculation of integrand function. The name FUNC - should be declared EXTERNAL in the calling subroutine.

NN,IPOINTS - are of the type INTEGER;
A,B,FUNC,S - are for DINT_BODE - REAL*8, for QINT_BODE - REAL*16.

Method:

The Bode quadrature formulas will be used. The calculation accuracy by step h
up to O (h ^ (2 [(IPOINTS+1) /2]).

Notes:

The version of the program with quadruple precision QINT_BODE exists only for the operating system UNIX.

References:

1. M. Abramowitz and I. Stegun, Handbook of mathematical functions,
   National Bureau of Standarts, NY, 1964.

Example:

       . . .
       IMPLICIT REAL*8 (A-H,O-Z)
       DIMENSION DINTVAL(3)
       EXTERNAL FUNC
       A=1.D0
       B=3.D0
       NN=2
       IPOINTS=6 
       WRITE(*,*) 'INTERVAL','  ','[',A,B,']'
       WRITE(*,*) 'NUMBER OF POINTS QUADRATURE RULE OF BODE=',IPOINTS
       DO I=1,3
          DINTVAL(I)=DINT_BODE(FUNC,A,B,NN,IPOINTS)
          WRITE(*,*) 'NN=',NN,'   ','INTVAL=',DINTVAL(I)
          NN=NN*2
       ENDDO
       SIGMA=(DINTVAL(1)-DINTVAL(2))/(DINTVAL(2)-DINTVAL(3))
       EXVAL=(DINTVAL(1)-SIGMA*DINTVAL(2))/(1.D0-SIGMA)
       WRITE(*,*) 'RUNGE COEFFICIENT=',SIGMA
       WRITE(*,*) 'EXTRAPOLATION OF INTEGRAL=',EXVAL
       EXACT=DEXP(-A)-DEXP(-B)
       WRITE(*,*) 'EXACT OF INTEGRAL=        ',EXACT
       WRITE(*,*) 'ERROR=        ',EXACT-EXVAL
       . . .
       FUNCTION FUNC(X)
       IMPLICIT REAL*8 (A-H,O-Z) 
       FUNC=DEXP(-X)
       RETURN
       END

Result:

       INTERVAL  [       1.000000000000000        3.000000000000000]
       NUMBER OF POINTS QUADRATURE RULE OF BODE =          6
       NN = 2   INTVAL =  3.180924625011772E-001
       NN = 4   INTVAL =  3.180923742385016E-001
       NN = 8   INTVAL =  3.180923728261328E-001
       RUNGE COEFFICIENT =      62.492653158681330
       EXTRAPOLATION OF INTEGRAL =  3.180923728031647E-001
       EXACT OF INTEGRA =          3.180923728035784E-001
       ERROR =          4.136690989753333E-013

---