[数值问题] 辛普森积分出现错误，error #6637

external f
    double precision fc,a,b,s
    a=0.0
    b=1.0
    eps=0.000001
    call simp(a,b,fc,eps,s)
    write(*, '(1x, "s=", d15.6)')s
    end
    function fc(x)
    double precision fc,x
    fc=log(1.0d0+x)/(1.0d0+x*x)
    return
    end
    
    subroutine simp(a,b,fc,eps,t)
    double precision a,b,fc,t,h,t1,s1,p,x,t2,s2
      
      n=1
      h=b-a
      t1=h*(fc(a)+fc(b))/2.0d0
      s1=t1
10    p=0.0d0
      do  20 k=0,n-1
          x=a+(k+0.5d0)*h
          p=p+fc(x)
20    continue
      t2=(t1+h*p)/2.0d0
      s2=(4.0d0*t2-t1)/3.0d0
      if (abs(s2-s1) .ge. eps)then
          t1=t2
          n=n+n
          h=h/2.0d0
          s1=s2
          goto 10
      endif
      t=s2
      return
    end