program trapezoidal_demo
C
C   this program illustrates integration of the exponential function
C   using the trapezoidal rule
C
C   CWJ SDSU 2/3/2005
C

      implicit none

C....... define lattice................

      real a,b	! endpoints
      real dx    ! 
      integer N  ! # of points

C....... define array.................

      integer size
      parameter (size=500)
      real f(0:size)

C........ dummy indices..............

      integer i
      real x    
C......... integral................
 
      real sum
      real exactsum   ! for comparison

C............... NOW BEGIN..........


      print*,' Demonstration of trapezoidal rule '
      print*,' on function exp(x) on interval (a,b) '
      print*,' '
    
C...............ENTER PARAMETERS.......

      print*,' enter a,b of interval '
      read*,a,b
    
      print*,' Enter # of points to be used '
      read*,N
 
      dx = (b-a)/float(n)

C.............. NOW FILL OUT ARRAY...........

      do i = 0,N
  	x = a + i*dx
        f(i) = exp(x)
      enddo
      
C..............INITIALIZE INTEGRAL .........

      sum = 0.0    

C..............apply trapezoidal rule......

      do i = 1,N
	sum = sum + 0.5*dx*(f(i)+f(i-1))
      enddo

C..............print results

      exactsum = exp(b)-exp(a)
      print*,' exact = ',exactsum,' approx = ',sum,
     & ' % error = ',100.*abs(exactsum-sum)/exactsum


      end