What's the name of this “shooting”-like method for numerical solving of PDE eigenproblems?

Consider a PDE of the following form:

$$\left(\partial_r^2+\frac5r\partial_r+\frac4{r^2} \hat L\right)\Psi(r,p)+(E- V(p)U(r))\Psi(r,p)=0,\tag1$$

where $$\hat L$$ is a differential operator differentiating only with respect to $$p$$ (i.e. ignoring $$r$$) and $$V$$, $$U$$ are some functions.

Suppose that boundary conditions are given such that $$|\Psi(0,p)|<\infty$$ and $$|\Psi(\infty,p)|<\infty$$, and some appropriate homogeneous conditions along the $$p$$ variable, independent from $$r$$.

I've read about a method for solving it numerically, in particular finding the eigenvalues, which I'll outline here.

First, a the solution $$\Psi(r,p)$$ is expanded into series over a set of orthogonal functions $$b_k(y)$$:

$$\Psi(r,p)=\sum\limits_{k=1}^\infty R_k(r) b_k(p),\tag2$$

where $$R_k$$ are $$r$$-dependent coefficients of expansion. These coefficients are collected into a column vector $$R(r)$$.

Expanding the $$p$$-dependent operators $$V(p)$$ and $$\hat L$$ in the same basis functions $$b_k$$, thus obtaining respectively the matrices $$\tilde V$$ and $$\tilde L$$, we get the following system of ODEs:

$$\left(\partial_r^2+\frac5r\partial_r+\frac4{r^2} \tilde L\right)R(r)+(I E-\tilde V U(r))R(r)=0,\tag3$$

where $$I$$ is the unit matrix

Then we introduce a matrix $$F(r)$$, which is defined as follows:

$$r R'(r)=F(r)R(r).\tag4$$

After inserting it into $$(3)$$, we get the following nonlinear equation for $$F(r)$$:

$$rF'(r)+F(r)^2+4F(r)+4\tilde L+(IE-\tilde V)r^2=0.\tag5$$

Now the equation $$(5)$$ automatically gives us an initial condition:

$$F(0)^2+4F(0)=-4\tilde L,\tag6$$

from which we can take (with some good definition of square root of a matrix):

$$F(0)=-2I+2\sqrt{I-L}.\tag7$$

Now (after truncating the series so that we work with finite-dimensional matrices) we can choose some value of $$E$$, and see if the solution of $$(5)$$ reaches a pole before we reach the second boundary. If it does, then we've found a possible zero of $$R(r)$$, and we have to lower $$E$$ to move the zero farther to the second boundary. Otherwise, $$E$$ is too small, and we should increase it to move the zero of $$R(r)$$ closer to the second boundary. Bisecting, we can find the eigenvalue as we do in the shooting method for ODE.

Other eigenvalues could be found by switching to an equation for $$G(r)=F(r)^{-1}$$ near the pole, and then resuming propagation for $$F(r)$$ itself.

My question: what is the name of this method? Where can I read more about it, its properties, pitfalls etc.? The source where I first encountered it simply described it similarly to my above description, without any justification of e.g. the choice of initial condition $$(7)$$ among all possible solutions of $$(6)$$, or any other explanations.