Giải số tích phân trên ngôn ngữ lập trình FORTRAN 95

Code tính tích phân bằng phương pháp hình thang mở rộng (expanded trapezoidal rule):

      PROGRAM Integral_Calculation

      REAL::a,b,h,TP,t
      INTEGER::n,i
     
      WRITE(*,*)'Nhap can duoi a = '
      READ(*,*)a
      WRITE(*,*)'Nhap can tren b = '
      READ(*,*)b
      WRITE(*,*)'Nhap so khoang chia n ='
      READ(*,*)n

      h=(b-a)/n
      TP=h*(f(a)+f(b))/2  

      DO i=1,n-1,1
        t=(b-a)*i/n
        TP=TP+h*f(t)
      END DO

      WRITE(*,*)'Gia tri tich phan la : ',TP

      END PROGRAM
     
      FUNCTION f(x)
      REAL::f,x
      f=1/(1+x)
      END FUNCTION f