C	THE PROGRAM RKPBD.F4 IS A SAMPLE PROGRAM,WHEN USED IN
C	CONJUNCTION WITH RKPBX(LOADED AND SAVED AS RKPBD.SAV)
C	WILL PRODUCE A RUNNING EXAMPLE OF THE USE OF RKPBX.
	SUBROUTINE RKPB2 (      TEMP,X,DX,Y,F,N)
       DIMENSION TEMP(1),Y(1),F(1),RKCON(4)
 401	M=N+1
       RKCON(1)=1.0
       RKCON(2)=2.0
       RKCON(3)=2.0
       RKCON(4)=1.0
       X=TEMP(2*M+1)
       DO 421 I=2,M
       IP2M=I+2*M
       Y(I-1)=TEMP(IP2M)
       IP3M=IP2M+M
       F(I-1)=TEMP(IP3M)
 421	TEMP(I)=0.0
       DO 491 J=1,4
       DO 441 I=2,M
       IPM=I+M
       TEMP(IPM)=F(I-1)*DX
 441	TEMP(I)=TEMP(I)+TEMP(IPM)*RKCON(J)
       IF(J-4) 451,501,5001
 451	X=TEMP(2*M+1)+DX/RKCON(J+1)
       DO 461 I=2,M
       IPM=I+M
       IP2M=IPM+M
 461	Y(I-1)=TEMP(IP2M)+TEMP(IPM)/RKCON(J+1)
       CALL DERIV
 491	CONTINUE
 501	X=TEMP(2*M+1)+DX
       DO 521 I=2,M
       IP2M=I+2*M
 521	Y(I-1)=TEMP(IP2M)+TEMP(I)/6.0
 5001	RETURN
       END
	SUBROUTINE RKPB1 (      TEMP,X,DX,Y,F,N)
       DIMENSION TEMP(1),Y(1),F(1)
 101	M=N+1
       TEMP(2*M+1)=X
       DO 121 I=2,M
       IP2M=I+2*M
 121	TEMP(IP2M)=Y(I-1)
       CALL DERIV
       TEMP(3*M+1)=DX
       DO 141 I=2,M
       IP3M=I+3*M
 141	TEMP(IP3M)=F(I-1)
 5001	RETURN
       END
CCCCCCCCCCCCCCCCCCCCCCCCCC

