Hello all, I've been writing a little physics code recently to solve the radial TISE for Hydrogen-like potentials using the Numerov Algorithm, basing my code on source code found on the internet. The program seems to work almost exactly as it should, however, if the radius is increased or the eigenvalue number is large then the algorithm seems to break down. I've included some of the key functions used in the technique.

The 'Solve' Method - For a given energy finds the nearest eigenvalue - Which in the context of Schrodinger's Equation is the energy

This is the 'Iterate' Method - For a given energy it calculates the wavefunction using Numerov's Method. The function 'Q' is the effective potential function. I know there are no problems with the 'Q' method.Code:double Numerov::Solve(bool dir) { double dE = 1.0e-02; double fE; double fE1; double E1, E2; E1 = E + dE; int it = 0; // Secant method to improve energy eigenvalue while (fabs(dE) > epsilon && it < 1000) { u1 = u0*exp(sqrt(-2.0*m_ec2*E)/h_bar_c*h); iterate(u0, u1, E, dir); fE = u[0]; u1 = u0*exp(sqrt(-2.0*m_ec2*E1)/h_bar_c*h); iterate(u0, u1, E1, dir); fE1 = u[0]; E2 = E1 - fE1*(E1 - E)/(fE1 - fE); E = E1; E1 = E2; dE = E1 - E; it++; } // calculate u with best energy value iterate(u0, u1, E, dir); Norm(); cout << endl; cout << "==> Energy difference: Rydberg - E = " << Ry - E << " after step: " << it << endl; return E; }

The array 'u' and 'v' represent the wavefunction and probability function respectively. 'h' is the step size, 'a' and 'b' are the starting and ending points considered (a is generally considered to be zero and so by 'large radii', I am referring to large 'b'). The data which it generates which is clearly wrong is given as output: '-1.#IND00000000'.Code:long int Numerov::iterate(double y0, double y1, double E, bool direction) { double hh12 = h*h/12.0; double x_i, x_im, x_ip; // Prograde recursion: This method is not working ! if (direction == true) { u[0] = 0.0; u[1] = y0; u[2] = y1; for (long int i=2; i < NGrid-2; i++) { x_i = a + i*h; x_im = x_i - h; x_ip = x_i + h; u[i+1] = ( 2.0*(1.0 - 5.0*hh12*Q(x_i, E))*u[i] - (1.0 + hh12*Q(x_im, E))*u[i-1] ) / (1.0 + hh12*Q(x_ip, E)); // Test condition: increasing wave function after decreasing if ((fabs(u[i-1]) > fabs(u[i])) && (fabs(u[i+1]) > fabs(u[i]))) return i; } } // Retrograde recursion else if (direction == false) { u[NGrid-1] = y0; u[NGrid-2] = y1; for (long int i=NGrid-2; i>=1; i--) { x_i = a + i*h; x_im = x_i - h; x_ip = x_i + h; u[i-1] = ( 2*(1.0 - 5.0*hh12*Q(x_i, E))*u[i] - (1.0 + hh12*Q(x_ip, E))*u[i+1] ) / (1.0 + hh12*Q(x_im, E)); } return 0; } }

Thanks in advance.