program Ecuaciones_Diferenciales_Foco
    ! Tarea #14
    ! Física Computacional (2:00-3:00)
    ! Programa elaborado por: Nalleli Iridian Avila Garcia (1941603).
    implicit none

    real(8) :: t0, tf, h
    real(8), dimension(2) :: u0, u, k1, k2, k3, k4
    integer :: i

    ! Imprime una descripción del programa
    print*, 'Considerando los sistemas de ecuaciones diferenciales: dx/dt=-y-x**2. dy/dt=x, dx/dt=-y-x**3. dy/dt=x.'
    print*, 'Utilizando el algoritmo de Runge-Kutta con diferentes condiciones iniciales, averiguamos si (0,0) es un centro o un foco para los dos sistemas dados.'

    ! Define las condiciones iniciales
    u0 = [0.0d0, 0.0d0]

    ! Define los parámetros de tiempo
    t0 = 0.0d0
    tf = 10.0d0
    h = 0.01d0

    ! Muestra las condiciones iniciales
    write(*,*) u0(1), u0(2)

    ! Realiza la integración de Runge-Kutta para el primer sistema de ecuaciones
    t = t0
    u = u0
    do i = 1, int((tf - t0) / h)
        k1 = h * sistema1(t, u)
        k2 = h * sistema1(t + 0.5d0 * h, u + 0.5d0 * k1)
        k3 = h * sistema1(t + 0.5d0 * h, u + 0.5d0 * k2)
        k4 = h * sistema1(t + h, u + k3)
        u = u + (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4) / 6.0d0
        t = t + h
        write(*,*) "Sistema 1:", u(1), u(2)
    end do

    ! Realiza la integración de Runge-Kutta para el segundo sistema de ecuaciones
    t = t0
    u = u0
    do i = 1, int((tf - t0) / h)
        k1 = h * sistema2(t, u)
        k2 = h * sistema2(t + 0.5d0 * h, u + 0.5d0 * k1)
        k3 = h * sistema2(t + 0.5d0 * h, u + 0.5d0 * k2)
        k4 = h * sistema2(t + h, u + k3)
        u = u + (k1 + 2.0d0 * k2 + 2.0d0 * k3 + k4) / 6.0d0
        t = t + h
        write(*,*) "Sistema 2:", u(1), u(2)
    end do

contains

    ! Define el primer sistema de ecuaciones diferenciales
    real(8) function sistema1(t, u)
        real(8), intent(in) :: t
        real(8), dimension(2), intent(in) :: u

        sistema1 = [-u(2) - u(1)**2, -u(1)]
    end function

    ! Define el segundo sistema de ecuaciones diferenciales
    real(8) function sistema2(t, u)
        real(8), intent(in) :: t
        real(8), dimension(2), intent(in) :: u

        sistema2 = [-u(2) - u(1)**3, -u(1)]
    end function

end program Ecuaciones_Diferenciales_Foco