一个算微分方程程序求助

yifanxing1992

module m_rk4
    !四阶龙格库塔方法计算微分方程
    contains
    subroutine rk4(t0,tt,y0,N)
    implicit real*8(a-z)
    integer::N,i
    !步长
    H=(tt-t0)/N
    t=t0
    y=y0
    do i=1,N
        call fun1(k1,t,y)
        call fun1(k2,t+H/2,y+H/2*k1)
        call fun1(k3,t+H/2,y+H/2*k2)
        call fun1(k4,t+H,y+H*k3)
        y=y+(k1+2*k2+2*k3+k4)*H/6
        t=t0+i*H
        write(11,21)t,y
    end do
21  format (f15.5,2f20.8)
    end subroutine rk4

    subroutine fun1(f,t,y)
    !待计算的函数模块，其中f(t,y)
    implicit real*8(a-z)
    real*8::pi
    real*8::S1,S2,p1,V1,l,h,Tk1,Tk2,T1,T2,rou,g
    pi=3.141592653

    S1=0.36*pi*10**(-6)
    S2=2.31*pi*10**(-5)
    p1=1.01*10**5
    l=0.345
    h=0.15
    V1=pi*(4.81*10**(-3))**2*(l-h)
    T1=20
    T2=80
    Tk1=T1+273.15
    Tk2=T2+273.15
    rou=10**3
    g=9.8

    f=rou**-1*((S1/S2)-1)**-1*((p1*V1*Tk2*S1)/(Tk1*(V1+S1*y))+rou*g*(1-S1/S2))
   ! f=sqrt(rou**-1*((S1/S2)**2-1)**-1*(2*(p1-(p1*V1*Tk2)/(Tk1*(V1+S1*y)))+2*rou*g*(h1+(1-S1/S2)*y)))
    end subroutine fun1   
end module m_rk4

program rspq
!龙格库塔方法主程序

use m_rk4
implicit real*8(a-z)
integer::N=1000
open(unit=11,file='result.txt')
t0=0
t1=1
y0=0
write(11,66)
66 format(/,T5,'计算结果为:',/)
call rk4(t0,t1,y0,N)

end program rspq

以上是原程序，初始条件和方程验证无书写错误，请问为何最后的数值结果一直显示y=NAN？？？？

li913

NAN = not a number.
计算过程中出现了数学错误，比如，我是说比如，对负数开方。
修改设置，使之报错，这样就知道是哪里出问题了。