Fortran program to calculate an integral using the method of trapezes

The following program calculates the integral I

using the trapezes method. In this case, n=3, x_ia=0, x_ib=1 and f(x)=exp(a·x).

  Program inttr
  double precision isum,x,xa,xb,h,f,I,a,IA,error,e
  integer j,p,k1,k2,k3,n,l
  dimension x(3),xa(3),xb(3),h(3),a(3)
  n=3
  p=300
  do j=1,n
      a(j)=rand()
      xa(j)=0
      xb(j)=1
      h(j)=(xb(j)-xa(j))/p
  enddo
  isum=0
  error=0
c we have to use n do's one inside the other. In this case, n=3.
c to calculate the error we need to compute the second partial
c derivate of F
  do k1=1,p
      x(1)=xa(1)+(k1-1/2)*h(1)
      do k2=1,p                      
          x(2)=xa(2)+(k2-1/2)*h(2)
          do k3=1,p
              x(3)=xa(3)+(k3-1/2)*h(3)
              call func(a,x,f,n)
              isum=isum+f
              do l=1,n
                  error=error+f*(a(l)**2)
              enddo
          enddo
      enddo
  enddo
  I=isum*(h(1)*h(2)*h(3))
  e=(p**(-2/n))*error/(24*isum)
  IA=1
  do j=1,n
      IA=IA*((exp(a(j))-1)/a(j))
  enddo
  write(*,*) I,e,IA        
  stop
  end



  SUBROUTINE func(a,x,f,n)
      double precision a,x,f,ax
      integer n,j
      dimension a(n),x(n)
      ax=0
      do j=1,n
          ax=ax+a(j)*x(j)
      enddo
      f=exp(ax)
  end