/*************************    Program 3.D    ****************************/
/*                                                                      */
/************************************************************************/
/* Please Note:                                                         */
/*                                                                      */
/* (1) This computer program is written by Tao Pang in conjunction with */
/*     his book, "An Introduction to Computational Physics," published  */
/*     by Cambridge University Press in 1997.                           */
/*                                                                      */
/* (2) No warranties, express or implied, are made for this program.    */
/*                                                                      */
/************************************************************************/

#include <stdio.h>
#include <math.h>
#define NMAX 501
#define v(x) (3-a*a*b*(b-1)/(2*pow(cosh(a*x),2)))
double a=1,b=4;

main()
/* Main program for solving the eigenvalue problem of the
   one-dimensional Schroedinger equation.
   Copyright (c) Tao Pang 1997. */
{
FILE *fout1,*fout2;
FILE *fopen();
int i,istep,im,n,m,imax,nl,nr;
double dl=1e-6;
double h,h2m,c,ea,eb,e,e1,de,xl0,xr0,xl,xr,f0,f1,fact,sum;
double ul[NMAX],ur[NMAX],ql[NMAX],qr[NMAX],s[NMAX];
void nmrv2 ();

fout1 = fopen("fort.15","a");
fout2 = fopen("fort.16","a");
n = NMAX;
m = 5;
imax = 100;
h2m = 0.5;
c   =  1/h2m;
xl0 = -10;
xr0 =  10;
h   =  (xr0-xl0)/(n-1);
ea  = 2.4;
eb  = 2.7;
e   = ea;
de  = 0.1;
ul[0] = 0;
ul[1] =  0.01;
ur[0] =  0;
ur[1] =  0.01;

/* Set up the potential q(x) and source s(x) */

for (i = 0; i < n; ++i)
  {
  xl = xl0+i*h;
  xr = xr0-i*h;
  ql[i] = c*(e-v(xl));
  qr[i] = c*(e-v(xr));
  s[i]  = 0;
  }

/* Find the matching point at the right turning point */

for (i = 0; i < n-1; ++i)
  {
  if (((ql[i]*ql[i+1]) < 0) && (ql[i] > 0)) im = i;
  }

/* Numerov algorithm from left to right and vice versa */

nl = im+2;
nr = n-im+1;
nmrv2 (nl,h,ql,s,ul);
nmrv2 (nr,h,qr,s,ur);

/* Rescale the left solution */
 
fact = ur[nr-2]/ul[im];
for (i = 0; i < nl; ++i)
  {
  ul[i] = fact*ul[i];
  }
f0 = ur[nr-1]+ul[nl-1]-ur[nr-3]-ul[nl-3];
f0 = f0/(2*h*ur[nr-2]);
 
/* Bisection method for the root */

istep = 0;
while ((fabs(de) > dl) && (istep < imax))
  {
  e1 = e;
  e  = (ea+eb)/2;
  for (i = 0; i < n; ++i)
    {
    ql[i] = ql[i]+c*(e-e1);
    qr[i] = qr[i]+c*(e-e1);
    }

/* Find the matching point at the right turning point */

  for (i = 0; i < n-1; ++i)
    {
    if (((ql[i]*ql[i+1]) < 0) && (ql[i] > 0)) im = i;
    }

/* Numerov algorithm from left to right and vice versa */

  nl = im+2;
  nr = n-im+1;
  nmrv2 (nl,h,ql,s,ul);
  nmrv2 (nr,h,qr,s,ur);

/* Rescale the left solution */
 
  fact = ur[nr-2]/ul[im];
  for (i = 0; i < nl; ++i)
    {
    ul[i] = fact*ul[i];
    }
  f1 = ur[nr-1]+ul[nl-1]-ur[nr-3]-ul[nl-3];
  f1 = f1/(2*h*ur[nr-2]);

  if ((f0*f1) < 0)
    {
    eb = e; 
    de = eb-ea;
    }
  else
    {
    ea = e;
    de = eb-ea;
    f0 = f1;
    }
  istep = istep+1;
  }
 
sum = 0;
for (i = 0; i < n; ++i)
  {
  if (i > im) ul[i] = ur[n-i-1];
  sum = sum+ul[i]*ul[i];
  }
sum = sqrt(h*sum);

printf("%4d %4d\n", istep,imax);
printf("%16.8lf %16.8lf %16.8lf %16.8lf\n", e,de,f0,f1);

for (i = 0; i < n; i = i+m)
  {
  xl = xl0+i*h;
  ul[i] = ul[i]/sum;
  fprintf(fout1,"%20.8lf %20.8lf\n", xl,ul[i]);
  fprintf(fout2,"%20.8lf %20.8lf\n", xl,v(xl));
  }
}

void nmrv2 (n,h,q,s,u)
/* The Numerov algorithm for the equation u"(x)+q(x)u(x)=s(x)
   as given in Eqs. (3.82)-(3.85) in the book.
   Copyright (c) Tao Pang 1997. */
int n;
double h;
double q[],s[],u[];
{
int i;
double g,c0,c1,c2,di,ut;

g = h*h/12;
for (i = 1; i <n-1; ++i)
  {
  c0 = 1+g*q[i-1];
  c1 = 2-10*g*q[i];
  c2 = 1+g*q[i+1];
  di = g*(s[i+1]+s[i-1]+10*s[i]);
  ut = c1*u[i]-c0*u[i-1]+di;
  u[i+1] = ut/c2;
  }
}