PROGRAM INTEGRAL
      external g 
      a2=0.0
      b2=1.0  
	print*,a2,b2,g(a2),g(b2)
      call simp2(g,a2,b2,s2)
      print*, s2, abs(s2-exp(1.0)+2.0)     

      STOP
      END

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	function g(yy)
	external f
	common/xy/x,y
	real g,x,y,yy
	y=yy
	a1=0.0
	b1=1.0
	call simp1(f,a1,b1,s1)
	g=s1
	return
	end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	function f(xx)
	common/xy/x,y
	real f,x,y,xx
	x=xx
	f=y*exp(x*y)
	return
	end
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SIMP1(F,A,B,S)
      EXTERNAL F
 
      PARAMETER (N =   9)
      H  = (B-A)/(N-1)
      S  = 0.
      S0 = 0.
      S1 = 0.
      S2 = 0.
      DO  x=a,b-1.0*h,2.0*h    
        S1 = S1 + F(x)  
        S0 = S0 + F(x+h)
        S2 = S2 + F(x+2.0*h)
      END DO   
      S = H*(S1 + 4.*S0 + S2)/3.
      RETURN
      END
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      SUBROUTINE SIMP2(g,A,B,S)
      EXTERNAL g
 
      PARAMETER (N =   9)
      H  = (B-A)/(N-1)
      S  = 0.
      S0 = 0.
      S1 = 0.
      S2 = 0.
      DO  x=a,b-1.0*h,2.0*h    
        S1 = S1 + g(x)  
        S0 = S0 + g(x+h)
        S2 = S2 + g(x+2.0*h)
      END DO   
      S = H*(S1 + 4.*S0 + S2)/3.
      RETURN
      END