getPsi()

void getPsi(double Psi[], double E, int l,
            double V[], double r[], double dr[], int N)

Output:

• Psi[k]: values of normalized $\Psi(r)$ on all grid points r[k]

Input:

• E: electron energy (found from ``getallEs'')
• l: angular momentum of $\Psi$
• V[]: potential on grid points
• r[], dr[], N: standard grid information

Write a routine getPsi() of the above prototype in ``physics.c'', which produces the values on the grid r[] of the wave function with energy E and angular momentum l, properly normalized so that $\int_0^\infty \vert\Psi(r)\vert^2\, dr=1$. Your routine will consist of several steps:

First, you should construct the appropriate effective potential,

$$V_{eff}(r)=V(r)+\frac{\ell(\ell+1)}{2r^2}.$$

Next, you need to determine an appropriate point, kmatch, for matching the inward and outward solutions. An appropriate value is the boundary between the exponential and oscillatory solutions, the so-called ``classical turning point'' where $V(r)=E$. To find the outermost point where $V(r)<E$, you may use the following fragment:

                        .
                        .
  kmatch=0;
  for (k=0; k<=N; k++)
    if (Veff[k]<E) kmatch=k;
  if (kmatch==0) {
    printf("Veff never below E=%f in getPsi.\n",E);
    exit(1);
  }
                        .
                        .

(This is the so-called ``classical turning point''.) Then, you should perform the inward (k=0..kmatch) and outward (k=N..kmatch) integrations using schint(), scaling both solutions to match so that $\Psi(r[\mbox{kmatch}])=1$.

Next, interpolate your result to include values at the odd points using your new routine interpolate().

Finally, normalize the result using the formula

$$\Psi(r) \leftarrow \frac{\Psi(r)}{\sqrt{\int_0^\infty \vert\Psi(r)\vert^2\,dr}},$$

where you perform the integral using simpint().