anchor.diff:  end function
   
   subroutine solve_radial_eigenproblem(n, l, Ein, eps, max_iter, &
  -    R, Rp, V, Z, c, relat, perturb, Emin_init, Emax_init, converged, E, Y)
  -! Solves the radial Dirac (Schrodinger) equation and returns the eigen value
  -! (E) and normalized eigenvector (Y) for the given "n" and "l".
  +    R, Rp, V, Z, c, relat, perturb, Emin_init, Emax_init, converged, E, P, Q)
  +! Solves the radial Dirac (Schroedinger) equation and returns the eigenvalue
  +! (E) and normalized eigenvector (P, Q) for the given "n" and "l".
   !
   !    Finds the wavefunction with defined "n" and "l". The potential is "V".
   !    relat ... 0 nonrelat (runge-kutta)
  
positive.code:     converged = 2
      return
  end if
  converged = 0
  Y = convert_PQ_to_R(R, Rp, P, Q, relat)
  end subroutine
  
  
  function convert_PQ_to_R(R, Rp, P, Q, relat) result(Y)
  ! Converts P, Q into R(r) depending on 'relat' and normalizes it
  real(dp), intent(in) :: R(:), Rp(:), P(:), Q(:)
  real(dp) :: Y(size(R))
  integer, intent(in) :: relat
  if (relat == 0) then
      Y = P
  else
      Y = sqrt(P**2 + Q**2)
  end if
  call normalize(Rp, Y)
  Y = Y / R
  end function
  
  integer function find_ctp(V, E) result(ctp)
  ! Finds the classical turning point for the potential 'V' and energy 'E'
  
negative.code: !        energy is not in the range, the "converged" variable will be equal to
  !        either 9 (Emin_init too big) or 10 (Emin_init too small).
  !
  !    Returns eigenvalue E and the normalized wavefunctions P, Q. If it doesn't
  !    converge, then converged /= 0, and E is undefined (P and Q then contains
  !    the latest integration of the shooting solver if available --- this is
  !    useful for debugging why it did not converged).
  !
  !    Conceptually, the algorithm is to use bisection to converge energy (using
  !    the number of nodes as the criterion) until we get close enough, so that
  !    the perturbation theory starts to converge, and then use perturbation to
  !    finish it. If perturbation does not converge, we need to fail over to
  !    bisection.
  integer, intent(in) :: n, l, relat, Z, max_iter
  real(dp), intent(in) :: R(:), Rp(:), V(:), eps, Ein, c
  logical, intent(in) :: perturb
  real(dp), intent(in) :: Emin_init, Emax_init
  integer, intent(out) :: converged
  real(dp), intent(out) :: P(:), Q(:), E
  
  
  real(dp) :: Emin, Emax, dE, Pr(size(R)), Qr(size(R)), factor, S
  integer :: minidx, ctp, iter
  logical :: isbig
  integer :: nnodes
  