program main

 implicit none

 real*8 :: a,b,eps
 real*8 :: f,x
 real*8 :: vv,ex,h
 real*8, allocatable :: m(:,:), c(:)

 integer :: i,j,n,np 
  
!1)
 write(*,*) 'Extremal points for integration'
 read(*,*) a,b

 write(*,*) 'Threshold value for convergence'
 read(*,*) eps

 write(*,*) 'Exact value'
 read(*,*) ex

!2)
 np=2
 h=(b-a)/dble(n)
 do
    vv=f(a)+f(b)
    vv=0.5d0*vv
    do i=1,n-1
       x = a+i*h
       vv=vv+f(x)
    enddo
    vv=vv*h
    if (abs(vv-ex).le.eps) then
       write(*,*) 'Convergence achieved!'
       exit
    endif

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

 enddo

!3)
  if (vv.gt.0.d0) then
    n=abs(int(b))+3
    allocate(m(n,n))
    do i=1,n
       do j=1,n
          m(i,j)=exp(-dble(i))*cos(dble(i*j))
       enddo
    enddo    
  else
    allocate(c(np))
    do i=1,np
       c(i)=vv*i**2/(i+1)
    enddo
  endif  

!4)
 write(*,*) 'Approximated integral value', vv
 write(*,*) 'Exact integral value', ex
 write(*,*) 'Absolute difference', abs(vv-ex)
 write(*,*) 'Needed number of points', np

 if (vv.gt.0.d0) then
    open(11,file='m.dat')
    do i=1,n
       do j=1,n
          write(11,100) i,j,m(i,j)
       enddo
    enddo      
    close(11) 
    deallocate(m) 
 else  
    open(11,file='c.dat')
    do i=1,np
       write(11,101) i, c(i)
    enddo   
    deallocate(c)
 endif

100 format(I5,I5,F12.5)
101 format(I5,F12.5)

 stop

end program main 


real*8 function f(x)

!Integral is 1/4*x^4 = 1/3*x^3 + x

  implicit none 

  real*8, intent(in) :: x

  f = x**3 - x**2 + 1

  return

end function f