In trying to solve a problem occurring in QM calculations I've encountered the following pickle, with which I hope you could help me. I am trying to solve an inhomogeneous eigenvalue differential equation: $$ y''(x) + (\lambda+f(x))y(x)=g(x) $$ with the boundary conditions $y(0) = y(\infty) = 0$.
I am using Numerov's algorithm to integrate $y(x)$, by interpolating with $y_{n+1} = F(y_n, y_{n-1})$, and I'm estimating the eigenvalue using the shooting method (that is by examining the tail of the function, until it matches the boundary condition at infinity).
The question is: I have to impose a second, more strenuous restriction on $y(x)$. As often occurs in QM, the norm of the function has to be unity. In other words: $$ \int y(x)^2 \textrm{d}\omega(x)=1 $$
Unfortunately, I couldn't figure out how to incorporate this condition, and retain the shooting method. Firstly, I tried adjusting the parameter that commences the integration; that is, for Numerov's method I have to specify two initial values for $y_0, y_1$. $y_1$ is generally arbitrary and I can calibrate it. But, if the trial eigenvalue itself is very distant from an actual eigenvalue, this diff. eq. admits an exponentially exploding solution which is not influenced by the parameter $y_1$. In trying to optimise the norm of the function in this case, my code stalls. Furthermore, I haven't exactly figured out how to make optimizing $y_1$ automatic, if I do somehow manage to place the initial eigenvalue close to an actual stationary state. In a nutshell, what is the best strategy you would recommend to add this second constraint on the BVP? What would a general code implementing it look like?
My first instinct here, which is rather biased based on my own experiences, would be to try a spectral method. That is, cut off the domain at some large $L$, artificially require $y(L)=0$, and then approximate $y$ as a truncated Fourier sine series. This seems nice for this application because the problem is linear, or at least is linear for each fixed $\lambda$, and the coefficient on the derivative is constant.
When you write the equation in Fourier space, you get
$$K \hat{y} + \widehat{fy} + \lambda \hat{y} = \hat{g} \\ \| \hat{y} \|^2_2 = 1.$$
Here $K$ is a diagonal matrix with $K_{jj}=-\frac{j^2 \pi^2}{L^2}$ and the norm in the second equation is just the Euclidean norm. This is technically an inhomogeneous matrix eigenvalue problem, because the map $\hat{y} \to \widehat{fy}$ is linear. It should really not be posed as a matrix eigenvalue problem, however. This is because the matrix defining $\hat{y} \mapsto \widehat{fy}$ is probably dense, which would mean that it takes $O(N^2)$ time to multiply it with a vector. By contrast the time to do the procedure $\hat{y} \to y \to fy \to \widehat{fy}$ is $O(N \log(N))$ using the fast Fourier transform. This means that good methods for solving it are black box methods which only need to evaluate $Ax$.
I actually expected that methods of this type would be quite standard, but on a quick search I did not find very much, especially among recent work. I did find this paper: http://www.sciencedirect.com/science/article/pii/0024379587901236 which deals with the problem by power iteration. (The first algorithm is stated on page 11.) One nice thing about their approach for this problem is that the normalization of the vector is essentially "built in".