!Driver routine for integration of ODE's using the Runge-Kutta 3rd order
!method with adaptive step-size, now generalized for systems of differential
!equations.

program RUNGE_KUTTA_general
use ADAPT_STEP_SIZE
use COUNTER
implicit none

     integer :: N
     integer, parameter :: pi = 3.141592654, d=2
     real :: a, b, x0, x1, y0(d), y1(d), h0, acc, eps, dummy(d), h_new, &
             & x0_new, y0_new(d), h_new2, x0_new2, y0_new2(d), resultat
     real, external :: f

!y is found in the y(1)'s and y_prime is found in the y(2)'s.
!Initial conditions: y(0)=0, y'(0)=1.

a = 0.0 ; b = 4.5*pi
h0 = 0.1
x0 = 0 ; y0(1) = 0 ; y0(2) = 1
acc = 0.00001 ; eps = 0.00001

OPEN(1, file="RK_general.out", status="replace", action="write")
write (1, *) "# x, y, y', step length and iteration:"
write (1, *) x0, y0(1), y0(2), h0 

call RK_general(f, x0, y0, h0, x1, y1, d)
call STEP_SIZE_general(f, a, b, x0, y0, h0, x1, y1, acc, eps, &
                       & dummy, h_new, d, x0_new, y0_new)

N=0

do
   write (1, *) x0_new, y0_new(1), y0_new(2), h_new, N
   N=N+1
   if (N==10000) STOP "Too many iterations!"
   if (x0_new == b) then
      resultat = y0_new2(1) - y0(1)
      exit
   endif

   call RK_general(f, x0_new, y0_new, h_new, x1, y1, d)
   call STEP_SIZE_general(f, a, b, x0_new, y0_new, h_new, x1, y1, acc, eps, &
                  & dummy, h_new2, d, x0_new2, y0_new2)
   h_new = h_new2
   x0_new = x0_new2
   y0_new = y0_new2
   if (x0_new + h_new > b) then
      h_new = b - x0_new
   endif 
enddo

print *, "Absolute accuracy", acc
print *, "Relative accuracy", eps
print *, "Value of integral:", resultat 
print *, "Number of function calls:" , X

CLOSE (1)
end program RUNGE_KUTTA_general

!arg1=x, arg2=y1, arg3=y2, j decides whether F1 or F2 is called.

real function f(arg1, arg2, arg3, j)
use COUNTER   
implicit none
        real :: arg1, arg2, arg3
        integer :: j
   if (j==1) then
      f = arg3
      else if(j==2) then
         f = -arg2
         else
            STOP "Sorry. No such function!"
   endif
   X = X + 1
end function f