One Dimensional Schrodinger's Equation solution using Numerov Method

Question

I have been trying to solve Time Independent Schrodinger's equation in one dimension using Numerov Method as discussed in this excellent lecture notes I found on net.

The Numerov method can solve an equation of the following kind:

$$\frac{{d^2}y}{dx^2}=-g(x) y(x) +s(x) $$

We can compare this with out Time Independent Schrodinger Equation :

$$\frac{{d^2}\psi}{dx^2}=-\frac{2m}{\bar h^2} \big(E-V) \psi$$

Where $ g(x)= -\frac{2m}{\bar h^2}(E-V) $

We can discretize the equation using Taylor expansion and write the iterative formula as :

$$\psi_{i+1}= \frac{\psi_i\big[2-\frac{5}{6}\cdot h^2 g_i]-\psi_{i-1}\big[1+\frac{h^2}{12}\cdot g_{i-1}] }{(1+\frac{h^2}{12}\cdot g_i)}$$

Using the known Energy of Harmonic Oscillator I solved the equation and my solutions were working well.

I also used the shooting method to determine Energy value for a particular Eigen State and got a great result.

I was succesfully able to solve the schrodinger equation for Most of Even potential including Morse Potential and The Double Well Potential.

The Problem arrises if I try to take a non symmetric solution. Because I have taken the symmetry of my wave function into account and just solved for $\psi$ from $x=0 $ to $ x=xmax$, and then took output in a file by multiplying a parity to the $\psi$ based on weather my function is odd or not.

But When I solved the same problem from $x=-xmax $ to $x=xmax$ using same algorithm so that I could devise a general algorithm for Numerov Method. The output started getting bizarre. While the ground state solution is same , the odd state solutions are getting inverted somehow.

---

This is a part of my program for the first (working) Numerov Algorithm. (ignore the signchange counter I have used it to determine number of nodes for shooting method)

cout<<"Enter mesh number"<<endl;
 cin>>n;

 cout<<"Enter maximum value of x"<<endl; 
 cin>>xmax;
 dx=xmax/n;

 if (nodes % 2==0) //Even Nodes, function is even function , i.e. y(-1) =      y(1) and f(-1) = f (1); Thus using algorithm!
        {
      y[0]=1.0;
      y[1]=( y[0]*( 12.0 - 10.0 *f[0] ) )/( 2*f[1] );
    }
  else //nodes are odd, i.e there is a node at x=0;
    {
      y[0]=0;
      y[1]=dx;           //arbitrary small value
    }

    signchange=0;

    //outward integration using algorithm and counting number of sign changes//

 for(int i=1; i <= n; i++) 

    {
      y[ i+1 ]= ( ( 12.0 - 10.0*f[i] )*y[i] -( y[i-1]*f[i-1] ) ) / f[i+1]; 

     if (y[i] != copysign(y[i],y[i+1]))  
       ++signchange;
    }
//For Output :

if (nodes % 2 == 0) 
    p=1;
    else
    p=-1;
    for (int i = n; i > 0; --i)
    file<<setw(1)<<-1*x[i]<<setw(15)<<p*y[i]<<setw(15)<<y[i]*y[i]<<setw(15)<<V[i]<<endl;

    //-------x<0-------//
    for (int i = 0; i <= n; ++i)
    file<<setw(1)<<x[i]<<setw(15)<<y[i]<<setw(15)<<y[i]*y[i]<<setw(15)<<V[i]<<endl;

---

And this is part of the second program which is giving me inverted output. Where I am iterating across the whole mesh

cout<<"Enter maximum value of x"<<endl;
cin>>xmax;

xmin=-1*xmax;

h=(2*xmax)/n;

y[0]=0;
    y[1]=h;

    /*calculating the wave function psi or y at all points for Energy e*/
        for(int i=1; i <= n; i++) 

    {
      y[ i+1 ]= ( ( 12.0 - 10.0*f[i] )*y[i] -( y[i-1]*f[i-1] ) ) / f[i+1]; 

     if (y[i] != copysign(y[i],y[i+1]))  
       ++signcount;
    }

I would love to know why this is happening. It does not make sense to me. But i am also a beginner in Computational Physics. Also I would like to know what should be the best way to solve the equation for Non-Symmetric potentials, for example a Coulumb potential (inverse R).

EDIT: Since I am trying to solve Harmonic Oscillator problem the wave function should be bound on both sides. So i am taking $\psi_0$ as 0. And $\psi_1$ as some random small value which will be taken in to account when the wave function will be normalized later. This is for the case when I am solving the equatioin across the whole mesh, from $-xmax$ to $xmax$ .

But When I am solving the equation for symetric potential(i.e. from 0 to xmax), I am using a property of solution that I have already known, that if i am solving for odd Energy state, the wave funcion $\psi_{x=0}$ will always be zero. I have provided my declaration of initial value in code too. I hope this helps, in clarifying my problem a little.