Numerical Integration using simpson rules and recusion in fortran

 There is a simple integration using fortran. For the given function (x^n*e^(x-1)) to integrate we can use recursion. For other functions recursion may not support. So if you change the function make sure it will support integration and if it is make sure given integration formula is correct.
Compilation
           gfortran <filename>.f
it will create default executable. in windows a.exe and in linux a.out.
simply type 'a' to run in wndows and './a' to run in linux.
If you have any problem comment me.
--------------------------------------------------------------------------------------------------------------
         Do you want to publish source codes in your blog or web site as follows.
                              Visit Source Code Formatter

    program integration  
 ! declaring the variable values  
    real results,simp13_res,simp38_res,a,b,error  
    real f  
    integer n  
    external f  
      results=0.0  
            n=0  
            a=0.0  
            b=1.0  
      error=0.0  
     write(*,*)achar(10),achar(10)  
     write(*,*)"Intergration of (x^n*e^(x-1)) by dx range of 0 to 1"   
 10   write(*,*)achar(10)  
     write(*,*)"Enter 0 to exit from the program"  
     write(*,*)"Enter n : "   
     read(*,'(i10)')n   
 !to exit            
     if(n==0)then  
            goto 20  
     end if  
 !factorial function calling  
     results=factorial(n)  !calling the function factorial  
     write(*,*)achar(10),"Recursion Solution = ",results  ! print the value of the function  
 !get simpson 1/3 solution  
     call simpson13(f,a,b,simp13_res,real(n))  
     write(*,*)achar(10),"Simpson 1/3 Solution= ",simp13_res  
                 error=(results-simp13_res)*100/results  
     write(*,*)"Relative True Error in Simpson 1/3 (%) = ",abs(error)  
                 error=0.0  
 !get simpson 3/8 solution                  
     call simpson38(f,a,b,simp38_res,real(n))  
     write(*,*)achar(10),"Simpson 3/8 Solution = ",simp38_res  
                 error=(results-simp38_res)*100/results  
     write(*,*)"Relative True Error in Simpson 3/8 (%) = ",abs(error)  
 !continue                  
    goto 10  
 20  continue  
    end program integration  
 ! end of the main program, functions are below which is used in this program  
 !recursive algorithm  
    recursive function factorial(n) result(results)  
            real results,first,x  
      integer n  
             x=1  
             first=1/exp(x)  
       if(n<=0) then  
        results = 0  
        return  
             else if(n==1) then  
                    results=1/exp(x)  
              return  
             else  
                    results=1-n*factorial(n-1)  
                    return  
             end if  
    end function factorial  
 !simpson 1/3 algorithm  
    Subroutine simpson13(f,a,b,simp38_res,n)  
            real a,b,x,n,h,f,simp38_res  
                  h=(b-a)/2.0  
      simp38_res=h*(f(a,n)+f(a+h,n)+f(b,n))/3.0  
      return  
    end Subroutine simpson13  
 !simpson 3/8 algorithm  
    Subroutine simpson38(f,a,b,simp13_res,n)  
            real a,b,x,n,h,f,simp13_res  
                  h=(b-a)/3.0  
      simp13_res=h*3.0*(f(a,n)+3.0*f(a+h,n)+3.0*f(a+2*h,n)+f(b,n))/8.0  
      return  
    end Subroutine simpson38  
 !function to integrate  
    function f(x,n) result(res)  
            real x,n,res  
            res=(x**n)*exp(x-1)  
      return  
    end function f