program quadratura_trapezio

!Formula di quadratura del trapezio:
! [a,b]            intervallo di integrazione      
! h = (b-a)/n      passo di integrazione 
! x(i) = a + i*h,  i=1,...,n-1 
! res = h*[f(a)/2 + f(x(1)) + ... + f(x(n-1)) + f(b)/2]

 implicit none

 real*8     :: a,b,h,res,x 
 integer    :: n,i
 real*8     :: f; external f
 real*8     :: thr, ex 

 write(*,*) 'Integration extremal points'
 read(*,*) a,b
 write(*,*) 'Convergence threshold'
 read(*,*) thr
 write(*,*) 'Exact value'
 read(*,*) ex 

 n=10
 h=(b-a)/dble(n) 
 do 
    res=f(a)+f(b)
    res=0.5d0*res
    do i=1,n-1
       x = a+i*h
       res=res+f(x) 
    enddo
    res=res*h
    write(*,*) 'Integral', res, 'with', n ,'points'
    if (abs(res-ex).le.thr) then
       write(*,*) 'Convergence achieved!'
       exit
    endif   

    n = n + 5
    h=(b-a)/dble(n) 

 enddo

 stop

end program quadratura_trapezio    


real*8 function f(x)

 implicit none

 real*8, intent(in)   :: x

 f = x**2

 return

end function f