! sim_difeq.f90  solve simultaneous differential equations
!                by Runge-Kutta by Gills method
!
! solve system of equations Y_i = Y'_i(X, Y_1, ..., Y_n) for i=1,n
!
! user sets N number of equations (in N dependent variables)
!           Y array of N dependent variables with initial values
!           F array of N first derivatives of Y's (computed in loop)
!           Q array of N for use by Runge
!           X independent variable with initial value
!           XLIM limit value of X with value
!           H step size, delta X with value
!           M index used in Runge, initial value zero
!           K returned 1=compute derivatives, 2=have Y values at X+H 

program sim_difeq
  implicit none
  integer, parameter :: N = 4          ! number of equations
  double precision, dimension(N) :: Y  ! dependent variables
  double precision, dimension(N) :: F  ! computed derivatives, Y'
  double precision, dimension(N) :: Q  ! temporaries for Runge
  double precision :: X, XLIM, H
  integer :: I, M, K

  interface
    subroutine Runge(N, Y, F, Q, X, H, M, K) ! Gills method
      integer, intent(in) :: N
      double precision, dimension(N), intent(inout) :: Y
      double precision, dimension(N), intent(in) :: F
      double precision, dimension(N), intent(inout) :: Q
      double precision, intent(inout) :: X
      double precision, intent(in) :: H
      integer, intent(inout) :: M, K
    end subroutine Runge
  end interface

  print *, 'sim_difeq.f90 solving ',N,' equations'
  X = 0.0D0      ! initial independent variable
  XLIM = 1.0D0   ! final value of independent variable
  H = 0.1D0      ! step size
  M = 0          ! initial for Runge
  do I=1,N       ! initial condition for dependent variables
    Y(I) = 0.0D0
  end do

  print *, 'simple numeric check, answer 1, 2, 3, 4'
  do while(X .lt. XLIM) ! simple numeric check
    K = 1
    do while(K .eq. 1)
      call Runge(N, Y, F, Q, X, H, M, K)
      F(1) = 1.0D0
      F(2) = 2.0D0
      F(3) = 3.0D0
      F(4) = 4.0D0
    end do ! X incremented by H
    print *, 'X=',X,'Y(1)=',Y(1),' Y(2)=',Y(2)
    print *, '                    Y(3)=',Y(3),' Y(4)=',Y(4)
  end do ! X at XLIM (done)

  print *, ' '
  print *, 'fourth order check 1, 1/2, 1/6, 1/24'
  X = 0.0D0      ! initial independent variable
  XLIM = 1.0D0   ! final value of independent variable
  H = 0.01D0     ! step size
  M = 0          ! initial for Runge
  do I=1,N       ! initial condition for dependent variables
    Y(I) = 0.0D0
  end do
  do while(X .lt. XLIM) ! fourth order check
    K = 1
    do while(K .eq. 1)
      call Runge(N, Y, F, Q, X, H, M, K)
      F(1) = 1.0D0
      F(2) = Y(1)
      F(3) = Y(2)
      F(4) = Y(3)
    end do ! X incremented by H
  end do ! X at XLIM (done)
  print *, 'X=',X,'Y(1)=',Y(1),'  Y(2)=',Y(2)
  print *, '                    Y(3)=',Y(3),' Y(4)=',Y(4)
  print *, 'Y(4)*24.0 =', Y(4)*24.0D0
  print *, ' '
  print *, 'sim_difeq.f90 finished'
end program sim_difeq

subroutine Runge(N, Y, F, Q, X, H, M, K) ! Gills method
  implicit none
  integer, intent(in) :: N
  double precision, dimension(N), intent(inout) :: Y
  double precision, dimension(N), intent(in) :: F
  double precision, dimension(N), intent(inout) :: Q
  double precision, intent(inout) :: X
  double precision, intent(in) :: H
  integer, intent(inout) :: M, K
  double precision :: A
  integer :: I

  M = M+1 ! integration step
  K = 1   ! exit condition
  if(M .eq. 1) then
    do I=1,N
      Q(I) = 0.0D0
    end do
  else if(M .eq. 2) then
    A = 0.5D0
    X = X+0.5D0*H
    do I=1,N
      Y(I) = Y(I) + A*F(I)*H
      Q(I) = 2.0D0*A*H*F(I)
    end do
  else if(M .eq. 3) then
    A = 0.2928932188134524756D0
    do I=1,N
      Y(I) = Y(I) + A*(F(I)*H-Q(I))
      Q(I) = 2.0D0*A*H*F(I) + (1.0D0-3.0D0*A)*Q(I)
    end do
  else if(M .eq. 4) then
    A = 1.7071067811865475244D0
    X = X+0.5D0*H
    do I=1,N
      Y(I) = Y(I) + A*(F(I)*H-Q(I))
      Q(I) = 2.0D0*A*H*F(I) + (1.0D0-3.0D0*A)*Q(I)
    end do
  else if(M .eq. 5) then
    A = 0.2928932188134524756D0
    do I=1,N
      Y(I) = Y(I) + H*F(I)/6.0D0 - Q(I)/3.0D0
    end do
    M = 0
    K = 2
  else
    print *, 'ERROR in Runge, bad M=',M
    M = 0
    K = 2
  end if
end subroutine Runge

! end sim_difeq.f90