! 

subroutine RK(f, a, b, x0, y0, h0, y1, acc, eps, d)
use MY_STUFF
Implicit none

     integer :: i, N
     integer, intent(in) :: d
     real(KIND=dbl) :: tol, h0_old, err(d), h0, x0, y0(d)
     real(KIND=dbl), intent(in) :: a, b, acc, eps
     real(KIND=dbl), intent(out) :: y1(d)
     real, external :: f

N=1

do
   call RKSTEP(f, x0, y0, h0, d, y1, err)
   tol = MAX(SQRT(h0/(b-a))*acc, SQRT(h0/(b-a))*eps*(MAXVAL(y1)))
   if (ABS(MAXVAL(err)) > (2.0*tol)) then
      h0 = h0*((tol/ABS((MAXVAL(err))))**0.25)*0.95
      print *, "Step:", h0   
      cycle
      else
         N = N + 1
         write (1, FMT='(4f10.5, i4)') x0+h0, y1(1), y1(2), h0, N
         x0 = x0 + h0
         if (x0 == b) then
            exit
         endif
         y0 = y1
         h0_old = h0
         h0 = h0*((tol/(ABS(MAXVAL(err))))**0.25)*0.95
         if (h0 > 2*h0_old) then
            h0 = 2*h0_old
         endif
         if ((h0 + x0) > b) then
            h0 = b - x0
         endif
   endif
enddo

end subroutine RK