Debugging Physics code is by far the worst code to debug!

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

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;
}

---

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:

---

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;
        }
}

---

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'.

Thanks in advance.

At higher eigenvalues the energy is greater and the wavelength is shorter. Have you tried decreasing the step size? What are the stability criteria for the Numerov algorithm?