program rk_solve
 implicit none
 integer , parameter                 :: dp = kind(0.0_8)
 real(dp), allocatable, dimension(:) :: t, x1, x2
 real(dp)                            :: ti, tf, x10,x20
 integer                             :: Nt
 integer                             :: i
!--------------------------------------------------------
!Input:
 print *,'# Runge-Kutta Method for 2-ODEs Integration'
 print *,'# Enter   Nt,ti,tf,x10,x20:'
 read  *,           Nt,ti,tf,x10,x20
 print *,'# Nt  = ',Nt ,'  dt  = ',(tf-ti)/(Nt-1)
 print *,'# ti  = ',ti ,'  tf  = ',tf
 print *,'# x10 = ',x10,'  x20 = ',x20
!Allocate arrays:
 ALLOCATE(t(Nt),x1(Nt),x2(Nt))
!Calculate:
 call rk
!Output:
 open(unit=11,file='rk.dat')
 do i=1,Nt
  write(11,*)T(i),X1(i),X2(i)
 enddo
 close(11)
!--------------------------------------------------------
contains
!--------------------------------------------------------
 function       f1(t,x1,x2)
  real(dp) :: t, x1, x2, f1 !t,x1,x2 are scalars in f1
  f1 = x2
 end function   f1
!--------------------------------------------------------
 function       f2(t,x1,x2)
  real(dp) :: t, x1, x2, f2 !t,x1,x2 are scalars in f2
  f2 = -10_dp * x1          !harmonic oscillator
 end function   f2
!--------------------------------------------------------
 subroutine     rk
  real(dp) :: dt, ts, x1s, x2s

  dt     = (tf-ti)/(Nt-1)

  t(1)   = ti ; x1(1) = x10 ; x2(1) = x20
  ts     = ti ; x1s   = x10 ; x2s   = x20

  do i   = 2, Nt
   call rkstep(ts,x1s,x2s,dt)
   t (i) = ti  + (i-1) * dt
   x1(i) = x1s ; x2(i) = x2s
  end do

 end subroutine rk
!--------------------------------------------------------
 subroutine     rkstep(t,x1,x2,dt)
  real(dp) :: t,x1,x2,dt
  real(dp) :: k11,k12,k13,k14,k21,k22,k23,k24
  real(dp) :: h,h2,h6

  
  h  =  dt         !h =dt, integration step
  h2 =  0.5_dp * h !h2=h/2
  h6 =h/6.0_dp     !h6=h/6
      
  k11=f1(t,x1,x2)
  k21=f2(t,x1,x2)
  k12=f1(t+h2,x1+h2*k11,x2+h2*k21)
  k22=f2(t+h2,x1+h2*k11,x2+h2*k21)
  k13=f1(t+h2,x1+h2*k12,x2+h2*k22)
  k23=f2(t+h2,x1+h2*k12,x2+h2*k22)
  k14=f1(t+h ,x1+h *k13,x2+h *k23)
  k24=f2(t+h ,x1+h *k13,x2+h *k23)

  t  =t+h
  x1 =x1+h6*(k11+2.0_dp*(k12+k13)+k14)
  x2 =x2+h6*(k21+2.0_dp*(k22+k23)+k24)
 

 end subroutine rkstep
!--------------------------------------------------------
end  program    rk_solve
!--------------------------------------------------------

!  ---------------------------------------------------------------------
!  Copyright by Konstantinos N. Anagnostopoulos (2004-2021)
!  Physics Dept., National Technical University,
!  konstant@mail.ntua.gr, www.physics.ntua.gr/~konstant
!  
!  This program is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, version 3 of the License.
!  
!  This program is distributed in the hope that it will be useful, but
!  WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  General Public License for more details.
!  
!  You should have received a copy of the GNU General Public Liense along
!  with this program.  If not, see <http://www.gnu.org/licenses/>.
!  -----------------------------------------------------------------------
