      SUBROUTINE QSIMP(ainteg1(X),A,B,S)
      PARAMETER (EPS=1.E-6, JMAX=20)
      save OST, OS
	external ainteg1
C
	print *, 'inside qsimp'

      OST=-1.E30
      OS= -1.E30
      DO 11 J=1,JMAX
        CALL TRAPZD(ainteg1(X),A,S,ST,J)
        S=(4.*ST-OST)/3.
        IF (ABS(S-OS).LT.EPS*ABS(OS)) RETURN
        OS=S
        OST=ST
11    CONTINUE
	print *, 'At end of qsimp',	J,A,B,S
      PAUSE 'Too many steps.'
      END
