Simpson 3/8 Rule | FORTRAN 95

Fortran Code

PROGRAM simpson38
IMPLICIT NONE
INTEGER::i,n
REAL::x0,xn,h,s,f

PRINT *,'====================================================' 
PRINT *,'Program for Simpson 3/8 rule [www.BottomScience.com]' 
PRINT *,'===================================================='

PRINT *,'Enter value of lower limit a?'
READ(*,*)x0

PRINT *,'Enter value of upper limit b?'
READ(*,*)xn

PRINT *,'Number of subintervals?'
READ(*,*)n

h=(xn-x0)/n
s=f(x0)+f(xn)
DO i=1,n-1
  IF (MOD(i,3)==0) THEN
    s=s+(2*f(x0+(i*h)))
  ELSE
    s=s+(3*f(x0+(i*h)))
  END IF
END DO

s=(3/8.)*h*s

PRINT *,"Value of integral is",s

END PROGRAM

REAL function f(x1)
REAL::x1
  !FUNCTION
  f=0.2+(25*x1)-200*(x1**2)+675*(x1**3)-900*(x1**4)+400*(x1**5)
return
end function