我的程序（龙哥库塔）
module m_rk_ods
integer,parameter::m=50000   
contains

    subroutine solve(func,t0,tt,y0,n)
    implicit real*8(a-z)

        external func
        integer::n
        real*8::y0(n),y(n)
        real*8::t0
        real*8::k1(n),k2(n),k3(n),k4(n)

        h=(tt-t0)/m
        t=t0
        y=y0


        do i=1,m
           call func(k1,t,y,n)
           call func(k2,t+h/2,y+h/2*k1,n)
           call func(k3,t+h/2,y+h/2*k2,n)
           call func(k4,t+h,y+h*k3,n)
           y=y+(k1+2*k2+2*k3+k4)*h/6
           t=t0+i*h
           write(11,101)t,y
        end do
          101 format(<m+1>f12.6)
    end subroutine solve

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    subroutine fun1(f,t,y,n)
    implicit real*8(a-z)
        integer::n
        real*8::f(n),y(n)

        f(1)=100*(y(1)-y(1)**3/3-y(2))
        f(2)=y(1)+1.02+0.05*y(3)
        f(3)=y(4)
        f(4)=-y(3)*(2*pi/5)**2

    end subroutine fun1


end module m_rk_ods

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
       program main
       use m_rk_ods
       implicit real*8(a-z)
           integer::n
           real*8::y(4),y0(4)

           open(unit=11,file='fout2.txt')
          
           n=4
           t0=0
           tt=50
           y0=(/-1.02,-0.67,0.00,0.34/) !第四个初始值随意设置的
           call solve(fun1,t0,tt,y0,n)
       end