给出一段使用runge-kutta方法的Fortran程序代码

下面是一个使用经典四阶Runge-Kutta方法求解常微分方程组的Fortran程序代码示例：

program runge_kutta_example
    implicit none
    integer, parameter :: n = 2 ! 方程组中的未知量个数
    real, parameter :: t0 = 0.0, tmax = 10.0, h = 0.01 ! 初始时间、最终时间和步长
    integer :: i, j, k
    real :: t, y(n), k1(n), k2(n), k3(n), k4(n), f(n)

    ! 初始条件
    y(1) = 1.0
    y(2) = 0.0

    ! 循环求解
    t = t0
    do while (t < tmax)
        ! 计算四个斜率
        do i = 1, n
            f(i) = dydt(i, t, y)
            k1(i) = h * f(i)
            k2(i) = h * dydt(i, t + 0.5*h, y + 0.5*k1)
            k3(i) = h * dydt(i, t + 0.5*h, y + 0.5*k2)
            k4(i) = h * dydt(i, t + h, y + k3)
        end do

        ! 更新解
        do i = 1, n
            y(i) = y(i) + (k1(i) + 2*k2(i) + 2*k3(i) + k4(i)) / 6.0
        end do

        ! 更新时间
        t = t + h

        ! 输出结果
        write(*, "(F8.4, 2(1X, F10.6))") t, y(1), y(2)
    end do

contains

    ! 定义常微分方程组
    function dydt(i, t, y) result(f)
        integer, intent(in) :: i
        real, intent(in) :: t, y(n)
        real :: f

        select case (i)
        case (1)
            f = y(2)
        case (2)
            f = -y(1)
        end select
    end function dydt

end program runge_kutta_example

在上面的代码中，`dydt`函数定义了常微分方程组，`k1`、`k2`、`k3`和`k4`分别表示四个斜率，`y`数组存储未知量的值，`f`数组存储`dy/dt`的值。在每个时间步长内，根据四个斜率求解新的解，并将时间向前推进一个步长`h`。程序最终输出的是`y(1)`和`y(2)`随时间的变化。