!Adaptive step-size (prediction-correction method)
!for integration of ODE's. 

module ADAPT_STEP_SIZE_pc
  interface
    recursive subroutine pc_STEP_SIZE(f, a, b, x1, y1, s, c, acc, eps, &
                                      & x2, y2, s_new, N)
       Implicit none

       integer, intent(in) :: N
       real :: tol, err
       real,intent(in) ::  a, b, x1, y1, s, acc, eps
       real,intent(out) :: x2, y2, s_new
       real, intent(inout) :: c
       real, external :: f
    end subroutine pc_STEP_SIZE
  end interface
end module ADAPT_STEP_SIZE_pc

recursive subroutine pc_STEP_SIZE(f, a, b, x1, y1, s, c, acc, eps, &
                                  & x2, y2, s_new, N)

Implicit none

     integer :: N
     real :: x1, y1, s, x2, y2, err, acc, eps, a, b, s_new, c, tol
     real, external :: f


call PC(f, x1, y1, s, c, x2, y2, err)

tol = acc*SQRT(s/(b-a)) + eps*ABS(y2)

if (ABS(err) > 2*tol) then 
   CALL pc_STEP_SIZE(f, a, b, x1, y1, s/2.0, c, acc, eps, x2, y2, s_new, N)
   else
      write (1, *) x2, y2, s, N
      s_new = s*((tol/err)**0.25)*0.95
      if (s_new > 2*s) then
          s_new = 2*s
      endif
endif

end subroutine pc_STEP_SIZE