	PROGRAM D11R5
C	Driver for routine BALANC
	PARAMETER(NP=5)
	DIMENSION A(NP,NP),R(NP),C(NP)
	DATA A/1.0,1.0,1.0,1.0,1.0,100.0,1.0,100.0,1.0,100.0,
     *		1.0,1.0,1.0,1.0,1.0,100.0,1.0,100.0,1.0,100.0,
     *		1.0,1.0,1.0,1.0,1.0/
C	Print norms
	DO 12 I=1,NP
		R(I)=0.0
		C(I)=0.0
		DO 11 J=1,NP
			R(I)=R(I)+ABS(A(I,J))
			C(I)=C(I)+ABS(A(J,I))
11		CONTINUE
12	CONTINUE
	WRITE(*,*) 'Rows:'
	WRITE(*,*) (R(I),I=1,NP)
	WRITE(*,*) 'Columns:'
	WRITE(*,*) (C(I),I=1,NP)
	WRITE(*,'(/1X,A/)') '***** Balancing Matrix *****'
	CALL BALANC(A,NP,NP)
C	Print norms
	DO 14 I=1,NP
		R(I)=0.0
		C(I)=0.0
		DO 13 J=1,NP
			R(I)=R(I)+ABS(A(I,J))
			C(I)=C(I)+ABS(A(J,I))
13		CONTINUE
14	CONTINUE
	WRITE(*,*) 'Rows:'
	WRITE(*,*) (R(I),I=1,NP)
	WRITE(*,*) 'Columns:'
	WRITE(*,*) (C(I),I=1,NP)
	END
	PROGRAM D3R10
C	Driver for routine BCUCOF
	DIMENSION C(4,4),Y(4),Y1(4),Y2(4)
	DIMENSION Y12(4),X1(4),X2(4)
	DATA X1/0.0,2.0,2.0,0.0/
	DATA X2/0.0,0.0,2.0,2.0/
	D1=X1(2)-X1(1)
	D2=X2(4)-X2(1)
	DO 11 I=1,4
		X1X2=X1(I)*X2(I)
		EE=EXP(-X1X2)
		Y(I)=X1X2*EE
		Y1(I)=X2(I)*(1.0-X1X2)*EE
		Y2(I)=X1(I)*(1.0-X1X2)*EE
		Y12(I)=(1.0-3.0*X1X2+X1X2**2)*EE
11	CONTINUE
	CALL BCUCOF(Y,Y1,Y2,Y12,D1,D2,C)
	WRITE(*,*) 'Coefficients for bicubic interpolation'
	DO 12 I=1,4
		WRITE(*,'(1X,4e15.6)') (C(I,J),J=1,4)
12	CONTINUE
	END
	PROGRAM D3R11
C	Driver for routine BCUINT
	DIMENSION Y(4),Y1(4),Y2(4),Y12(4),XX(4),YY(4)
	DATA XX/0.0,2.0,2.0,0.0/
	DATA YY/0.0,0.0,2.0,2.0/
	X1L=XX(1)
	X1U=XX(2)
	X2L=YY(1)
	X2U=YY(4)
	DO 11 I=1,4
		XXYY=XX(I)*YY(I)
		Y(I)=XXYY**2
		Y1(I)=2.0*YY(I)*XXYY
		Y2(I)=2.0*XX(I)*XXYY
		Y12(I)=4.0*XXYY
11	CONTINUE
	WRITE(*,'(/1X,T6,A,T14,A,T22,A,T28,A,T38,A,T44,A,T54,A,T60,A/)')
     *		'X1','X2','Y','EXPECT','Y1','EXPECT','Y2','EXPECT'
	DO 12 I=1,10
		X1=0.2*I
		X2=X1
		CALL BCUINT(Y,Y1,Y2,Y12,X1L,X1U,X2L,X2U,X1,X2,ANSY,ANSY1,ANSY2)
		X1X2=X1*X2
		EY=X1X2**2
		EY1=2.0*X2*X1X2
		EY2=2.0*X1*X1X2
		WRITE(*,'(1X,8F8.4)') X1,X2,ANSY,EY,ANSY1,EY1,ANSY2,EY2
12	CONTINUE
	END
	PROGRAM D6R26
C	Driver for routine BESSI
	CHARACTER TEXT*27
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Modified Bessel Function In') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T12,A,T20,A,T34,A)')
     *		'N','X','Actual','BESSI(N,X)'
	DO 11 I=1,NVAL
		READ(5,*) N,X,VALUE
		WRITE(*,'(1X,I4,F8.2,2E16.7)') N,X,VALUE,BESSI(N,X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R21
C	Driver for routine BESSI0
	CHARACTER TEXT*27
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Modified Bessel Function I0') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T13,A,T28,A)')
     *		'X','Actual','BESSI0(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,VALUE
		WRITE(*,'(F6.2,2E16.7)') X,VALUE,BESSI0(X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R23
C	Driver for routine BESSI1
	CHARACTER TEXT*27
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Modified Bessel Function I1') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T13,A,T28,A)')
     *		'X','Actual','BESSI1(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,VALUE
		WRITE(*,'(F6.2,2E16.7)') X,VALUE,BESSI1(X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R20
C	Driver for routine BESSJ
	CHARACTER TEXT*18
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Bessel Function Jn') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T12,A,T20,A,T33,A)')
     *		'N','X','Actual','BESSJ(N,X)'
	DO 11 I=1,NVAL
		READ(5,*) N,X,VALUE
		WRITE(*,'(1X,I4,F8.2,2E15.6)') N,X,VALUE,BESSJ(N,X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R15
C	Driver for routine BESSJ0
	CHARACTER TEXT*18
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Bessel Function J0') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A1,T12,A6,T22,A9)')
     *		'X','Actual','BESSJ0(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,VALUE
		WRITE(*,'(F6.2,2F12.7)') X,VALUE,BESSJ0(X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R17
C	Driver for routine BESSJ1
	CHARACTER TEXT*18
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Bessel Function J1') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A1,T12,A6,T22,A9)')
     *		'X','Actual','BESSJ1(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,VALUE
		WRITE(*,'(F6.2,2F12.7)') X,VALUE,BESSJ1(X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R25
C	Driver for routine BESSK
	CHARACTER TEXT*27
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Modified Bessel Function Kn') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T12,A,T20,A,T35,A)')
     *		'N','X','Actual','BESSK(N,X)'
	DO 11 I=1,NVAL
		READ(5,*) N,X,VALUE
		WRITE(*,'(1X,I4,F8.2,2E16.7)') N,X,VALUE,BESSK(N,X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R22
C	Driver for routine BESSK0
	CHARACTER TEXT*27
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Modified Bessel Function K0') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T13,A,T28,A)')
     *		'X','Actual','BESSK0(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,VALUE
		WRITE(*,'(F6.2,2E16.7)') X,VALUE,BESSK0(X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R24
C	Driver for routine BESSK1
	CHARACTER TEXT*27
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Modified Bessel Function K1') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T13,A,T28,A)')
     *		'X','Actual','BESSK1(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,VALUE
		WRITE(*,'(F6.2,2E16.7)') X,VALUE,BESSK1(X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R19
C	Driver for routine BESSY
	CHARACTER TEXT*18
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Bessel Function Yn') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T12,A,T20,A,T33,A)')
     *		'N','X','Actual','BESSY(N,X)'
	DO 11 I=1,NVAL
		READ(5,*) N,X,VALUE
		WRITE(*,'(1X,I4,F8.2,2E15.6)') N,X,VALUE,BESSY(N,X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R16
C	Driver for routine BESSY0
	CHARACTER TEXT*18
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Bessel Function Y0') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A1,T12,A6,T22,A9)')
     *		'X','Actual','BESSY0(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,VALUE
		WRITE(*,'(F6.2,2F12.7)') X,VALUE,BESSY0(X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R18
C	Driver for routine BESSY1
	CHARACTER TEXT*18
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Bessel Function Y1') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A1,T12,A6,T22,A9)')
     *		'X','Actual','BESSY1(X)'
	DO 11 I=1,NVAL
		READ(5,*) X,VALUE
		WRITE(*,'(F6.2,2F12.7)') X,VALUE,BESSY1(X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R5
C	Driver for routine BETA
	CHARACTER TEXT*13
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Beta Function') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A1,T11,A1,T24,A6,T43,A9)')
     *		'W','Z','Actual','BETA(W,Z)'
	DO 11 I=1,NVAL
		READ(5,*) W,Z,VALUE
		WRITE(*,'(2F6.2,2E20.6)') W,Z,VALUE,BETA(W,Z)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R13
C	Driver for routine BETAI,BETACF
	CHARACTER TEXT*24
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Incomplete Beta Function') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T5,A,T15,A,T27,A,T36,A,T47,A)')
     *		'A','B','X','Actual','BETAI(X)'
	DO 11 I=1,NVAL
		READ(5,*) A,B,X,VALUE
		WRITE(*,'(F6.2,4F12.6)') A,B,X,VALUE,BETAI(A,B,X)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D6R3
C	Driver for routine BICO
	CHARACTER TEXT*21
	OPEN(5,FILE='FNCVAL.DAT',STATUS='OLD')
10	READ(5,'(A)') TEXT
	IF (TEXT.NE.'Binomial Coefficients') GOTO 10
	READ(5,*) NVAL
	WRITE(*,*) TEXT
	WRITE(*,'(1X,T6,A1,T12,A1,T19,A6,T28,A9)')
     *		'N','K','Actual','BICO(N,K)'
	DO 11 I=1,NVAL
		READ(5,*) N,K,BINCO
		WRITE(*,'(2I6,2F12.0)') N,K,BINCO,BICO(N,K)
11	CONTINUE
	CLOSE(5)
	END
	PROGRAM D7R9
C	Driver for routine BNLDEV
	PARAMETER(N=20,NPTS=1000,ISCAL=200,NN=100)
	DIMENSION DIST(21)
	CHARACTER TEXT(50)*1
	IDUM=-133
	LLEN=50
10	DO 11 J=1,21
		DIST(J)=0.0
11	CONTINUE
	WRITE(*,*) 'Mean of binomial distribution (0 to 20) (neg to end)'
	READ(*,*) XM
	IF (XM.LT.0) GOTO 99
	PP=XM/NN
	DO 12 I=1,NPTS
		J=INT(BNLDEV(PP,NN,IDUM))
		IF ((J.GE.0).AND.(J.LE.20)) DIST(J+1)=DIST(J+1)+1
12	CONTINUE
	WRITE(*,'(1X,T5,A,T10,A,T18,A)') 'x','p(x)','graph:'
	DO 15 J=1,20
		DIST(J)=DIST(J)/NPTS
		DO 13 K=1,50
			TEXT(K)=' '
13		CONTINUE
		TEXT(1)='*'
		KLIM=INT(ISCAL*DIST(J))
		IF (KLIM.GT.LLEN) KLIM=LLEN
		DO 14 K=1,KLIM
			TEXT(K)='*'
14		CONTINUE
		WRITE(*,'(1X,F5.1,F8.4,3X,50A1)') FLOAT(J-1),DIST(J),
     *			(TEXT(K),K=1,50)
15	CONTINUE
	GOTO 10
99	END
	PROGRAM D10R3
C	Driver for routine BRENT
	PARAMETER(TOL=1.0E-6,EQL=1.E-4)
	DIMENSION AMIN(20)
	EXTERNAL FUNC
	NMIN=0
	WRITE(*,'(/1X,A)') 'Minima of the function BESSJ0'
	WRITE(*,'(/1X,T6,A,T19,A,T28,A,T40,A/)') 'Min. #','X',
     *		'BESSJ0(X)','BESSJ1(X)'
	DO 12 I=1,100				 
		AX=I
		BX=I+1.0
		CALL MNBRAK(AX,BX,CX,FA,FB,FC,FUNC)
		B=BRENT(AX,BX,CX,FUNC,TOL,XMIN)
		IF (NMIN.EQ.0) THEN
			AMIN(1)=XMIN
			NMIN=1
			WRITE(*,'(1X,5X,I2,3X,3F12.6)') NMIN,XMIN,
     *				BESSJ0(XMIN),BESSJ1(XMIN)
		ELSE
			IFLAG=0
			DO 11 J=1,NMIN
				IF (ABS(XMIN-AMIN(J)).LE.EQL*XMIN)
     *					IFLAG=1
11			CONTINUE
			IF (IFLAG.EQ.0) THEN
				NMIN=NMIN+1
				AMIN(NMIN)=XMIN
				WRITE(*,'(1X,5X,I2,3X,3F12.6)') NMIN,
     *					XMIN,BESSJ0(XMIN),BESSJ1(XMIN)
			ENDIF
		ENDIF
12	CONTINUE
	END

	FUNCTION FUNC(X)
	FUNC=BESSJ0(X)
	RETURN
	END
	PROGRAM D15R6
C	Driver for routine BSSTEP
	EXTERNAL DERIVS
	PARAMETER(N=4)
	DIMENSION Y(N),DYDX(N),YSCAL(N)
	X=1.0
	Y(1)=BESSJ0(X)
	Y(2)=BESSJ1(X)
	Y(3)=BESSJ(2,X)
	Y(4)=BESSJ(3,X)
	DYDX(1)=-Y(2)
	DYDX(2)=Y(1)-Y(2)
	DYDX(3)=Y(2)-2.0*Y(3)
	DYDX(4)=Y(3)-3.0*Y(4)
	DO 11 I=1,N
		YSCAL(I)=1.0
11	CONTINUE
	HTRY=1.0
	WRITE(*,'(/1X,T8,A,T19,A,T31,A,T43,A/)') 
     *		'eps','htry','hdid','hnext'
	DO 12 I=1,15
		EPS=EXP(-FLOAT(I))
		CALL BSSTEP(Y,DYDX,N,X,HTRY,EPS,YSCAL,HDID,HNEXT,DERIVS)
		WRITE(*,'(2X,E12.4,F8.2,2X,2F12.6)') EPS,HTRY,HDID,HNEXT
12	CONTINUE
	END

	SUBROUTINE DERIVS(X,Y,DYDX)
	DIMENSION Y(1),DYDX(1)
	DYDX(1)=-Y(2)
	DYDX(2)=Y(1)-(1.0/X)*Y(2)
	DYDX(3)=Y(2)-(2.0/X)*Y(3)
	DYDX(4)=Y(3)-(3.0/X)*Y(4)
	RETURN
	END
