[数值问题] 照抄的程序不知道哪里不对劲，结果竟然是0

module newton_cotes
contains

subroutine solve(func,s,a,b,n)
implicit real*8(a-z)
external func
integer :: n,i
hstep = (b-a)/n
s=0

do i=1,n
c=a+i*hstep
d=a+i*hstep
call cotes(func,s1,c,d)
s=s+s1
end do
end subroutine

subroutine cotes(func,s1,c,d)
implicit real*8(a-z)
external func
h=(d-c)/500
call func(f1,c)
call func(f2,c+h)
call func(f3,c+2d0*h)
call func(f4,c+3d0*h)
call func(f5,c+4d0*h)
call func(f6,d)
s1=19d0*f1+75d0*f2+50d0*f3+50d0*f4+75d0*f5+19d0*f6
s1=s1*5d0*h/288d0
end subroutine cotes

subroutine fun1(f,x)
implicit real*8(a-z)
f=x**2+dsin(x)
end subroutine fun1
end module newton_cotes

program main
use newton_cotes
implicit real*8(a-z)
integer :: n
open(unit=11,file='result.txt')
write(11,101)
call solve(fun1,s,-2d0,2d0,8)
write(11,102) s
101 format(/,T5,'复合高阶Newton_Cotes法计算数值积分',/)
102 format(T5,'积分结果',F15.8)
end program main

c=a+i*hstep
d=a+i*hstep

c=a+(i-1)*hstep
d=a+i*hstep