! Programa #07
! Física Computacional (2:00-3:00 pm)
! Elaborado por Nalleli Iridian Ávila García
! Este programa calcula mediante el método de Runge Kutta la ecuación diferencial dy/dx=2xy
! Bajo las condiciones y(0)=1, en donde x [0,1], con n=5.

program MetodoRungeKutta
    implicit none
    double precision :: x0, y0, xn, h, x(6), y(6), k1, k2, k3, k4
    integer :: n, i

    ! Parámetros
    x0 = 0.0
    y0 = 1.0
    xn = 1.0
    n = 5
    h = (xn-x0)/dble(n)

    ! Inicializar arrays
    x = (/ (x0 + i*h, i=0,n) /)
    y(1) = y0

    ! Método de Runge-Kutta de cuarto orden
    do i = 1, n
        k1 = h * f(x(i), y(i))
        k2 = h * f(x(i) + 0.5*h, y(i) + 0.5*k1)
        k3 = h * f(x(i) + 0.5*h, y(i) + 0.5*k2)
        k4 = h * f(x(i) + h, y(i) + k3)
        y(i+1) = y(i) + (k1 + 2*k2 + 2*k3 + k4)/6
    end do

    ! Imprimir los resultados
    do i = 1, n+1
        print*, "x = ", x(i), ", y = ", y(i)
    end do
end program MetodoRungeKutta

! Definir la función
double precision function f(x, y)
    implicit none
    double precision, intent(in) :: x, y
    f = 2*x*y
end function f