Modifying ``schint''

int schint(double *Psi, double *Psip, double Psiout[],
           int k1, int k2, double Veff[], double E,
           double r[], double dr[], int N);

Return value:

• number of nodes in solution between r[k1] and r[k2]

Output:

• Psi, Psip (passed by reference): values of and
• Psiout[]: values of generated by ``rk4p480''

Input:

• Psi, Psip (passed by reference): values of and
• k1, k2: indices of initial and final integration points (r[k1] r[k2])
• Veff[]: effective potential on grid
• E: proposed electron energy
• r[], dr[], N: standard grid information

Modify your ``schint'' routine to return the solution for for all points k=k1..k2. This involves two steps.

• First, if Psiout is not ``NULL'' we should record the values of the solution in Psiout:

                             .
                             .
    for (k=k1; k!=k2; k+=dk) {
      x=0.5*k*fabs(h);
      if (Psiout!=NULL) Psiout[k]=*Psi;
      derivs_Schrodinger(x,y,dydx,k,T,r,dr,N);
      rk4p480(y,dydx,2,x,h,y,derivs_Schrodinger,k,dk,T,r,dr,N);
                             .
                             .
    }
    if (Psiout!=NULL) Psiout[k]=*Psi;
                             .
                             .

  Note that because we record the initial value of Psi before ``rk4p480'' makes each step, we require one final assignment after the completion of the loop.
• Second, each time we scale the solution to keep Psi within 1e-150 to 1e150 (See Section 1.), we must be sure to scale all of the results stored in Psiout. But, to avoid losses in efficiency, this should only be done if we are returning Psioutp[]. A suitable code fragment, which works for both k1<k2 and k1>k2, is

      if (fabs(*Psi)>SQRTBIG) {
        y[1]=SMALL*y[1];
        y[2]=SMALL*y[2];
        if (Psiout!=NULL) { /* If filling array, need to repair */
  	for (kp=k; kp!=k1; kp-=dk)
  	  Psiout[kp]=SMALL*Psiout[kp];
  	Psiout[kp]=SMALL*Psiout[kp];
        }
      }