Simpson 1/3rd Rule | FORTRAN 95

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

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

PRINT *,'================'
PRINT *,'SIMPSON 1/3RD'
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

IF (MOD(n,2)==0) THEN
h=(xn-x0)/n
s=f(x0)+f(xn)+4*f(x0+h)
DO i=3,n-1,2
s=s+(4*f(x0+(i*h)))+(2*f(x0+(i-1)*h))
END DO
s=(h*s)/3
PRINT *,"Value of integral is",s
ELSE
PRINT *,"Number of interval is not even"
END IF

END PROGRAM