implicit none

 integer :: f,n
 real(8) :: p,x1,x2,s
 real(8), external :: power,exponential

 1 print*,'Function (1=power, 2=exponential):'; read*,f; 
 if (f/=1.and.f/=2) goto 1
 print*,'Give power p (in x**p or exp(p*x)):'; read*,p
 print*,'Integrate between values (x1,x2):'; read*,x1,x2
 print*,'Number n of integration intervals:'; read*,n
 if (mod(n,2)/=0) then
   n=n+1; print*,'Using n = ',n
 end if
 if (f==1) then
    call simpson(x1,x2,n,power,p,s)
 else
    call simpson(x1,x2,n,exponential,p,s)
 endif
 print*,'The Simpson integral is: ',s

 end

!------------------------------------!
 subroutine simpson(x1,x2,n,func,p,s)

 implicit none
 
 integer :: i,n
 real(8) :: x1,x2,dx,p,s
 real(8), external :: func
 
 dx=(x2-x1)/dble(n)
 s=func(x1,p)
 do i=1,n-1,2
   s=s+func(x1+dx*i,p)*4.d0+func(x1+dx*(i+1),p)*2.d0
 end do
 s=(s-func(x2,p))/(3.d0*n)

 end subroutine simpson

!-------------------!
 function power(x,p)
 implicit none
 real(8) :: power,x,p
 power=x**p
 end function power

!-------------------------!
 function exponential(x,p)
 implicit none
 real(8) :: exponential,x,p
 exponential=exp(p*x)
 end function exponential