
      program trapezoidalSimpson
      implicit real*8 (a-h,o-z)
      parameter (N=10)
      dimension x(0:N),f(0:N)
      
***********************
*     INSERT DATA     *
***********************
      x(0)=0d0
      x(N)=1d0
      h=(x(N)-x(0))/N
      print*,'N   =',N
      print*,'x(0)=',x(0)
      print*,'x(N)=',x(N)
      print*,'h   =',h
      print*
      
      do i=0,N
         x(i)=x(0)+i*h
         f(i)=dexp(x(i))
         print *,i,x(i),f(i)
      enddo
      print*
      
*****************************
*     SOLVE THE PROBLEM     *
*****************************
*     1.   TRAPEZOIDAL RULE
      sum=0d0
      do i=0,N-1,1
         sum=sum+h/2d0*(f(i)+f(i+1))
      enddo
      Etrap=sum
      
*     2.   SIMPSON'S RULE
      sum=0d0
      do i=0,0                  !CORRECT HERE
         sum=sum                !CORRECT HERE, USE EQ.(5.8)
      enddo
      ESimp=sum
      
*     3.   SIMPSON'S 3/8 RULE
      sum=0d0
      do i=0,0                  !CORRECT HERE
         sum=sum                !CORRECT HERE, USE EQ.(5.11)
      enddo
      ESimp38=sum
      
****************************************************
*     ANALYTIC (EXACT) SOLUTION FOR COMPARISON     *
****************************************************
      Eexact=dexp(1d0)-1d0
      
*************************
*     PRINT RESULTS     *
*************************
      print*,'E (exact )=',Eexact
      print*,'E (trap  )=',Etrap
      print*,'E (Simp  )=',ESimp
      print*,'E (Simp38)=',ESimp38
      print*
      
***************************
*     ESTIMATE ERRORS     *
***************************
*     For trapezoidal rule, Equation (5.15)
      print*,'Trap'
      error_min=-(x(N)-x(0))/12d0*h**2d0*dexp(x(0))
      error_max=-(x(N)-x(0))/12d0*h**2d0*dexp(x(N))
      error=Eexact-Etrap
      print*,'Minimum estimated error=',error_min
      print*,'Achieved          error=',error
      print*,'Maximum estimated error=',error_max
      print*
      
*     For Simpson's rule, Equation (5.17)
      print*,'Simp'
      error_min=0d0                     !CORRECT HERE
      error_max=0d0                     !CORRECT HERE
      error=0d0                         !CORRECT HERE
      print*,'Minimum estimated error=',error_min
      print*,'Achieved          error=',error
      print*,'Maximum estimated error=',error_max
      print*
      
*     For Simpson's 3/8 rule, Equation (5.19)
      print*,'Simp38'
      error_min=0d0                     !CORRECT HERE
      error_max=0d0                     !CORRECT HERE
      error=0d0                         !CORRECT HERE
      print*,'Minimum estimated error=',error_min
      print*,'Achieved          error=',error
      print*,'Maximum estimated error=',error_max
      print*
      
      pause
      end
