/*************************    Program 3.C    ****************************/
/*                                                                      */
/************************************************************************/
/* 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 101

main()
/* Program for the eigenvalue problem with a combination of the
   bisection method and the Numerov algorithm for u" = -k**2*u
   with boundary conditions u(0)=u(1)=0.
   Copyright (c) Tao Pang 1997. */
{
int i,istep,n;
double dl=1e-6;
double h,ak,bk,dk,ek,f0,f1;
double q[NMAX],s[NMAX],u[NMAX];
void nmrv();

n = NMAX;
h = 1.0/(n-1);
ak = 2;
bk = 4;
dk = 1;
istep = 0;
u[0] = 0;
u[1] = 0.01;

ek = ak;
for (i = 0; i < n; ++i)
  {
  s[i] = 0;
  q[i] = ek*ek;
  }
nmrv (n,h,q,s,u);
f0 = u[n-1];
/* Bisection method for the root */

while (fabs(dk) > dl)
  {
  ek = (ak+bk)/2;
  for (i = 0; i < n; ++i)
    {
    q[i] = ek*ek;
    }
  nmrv (n,h,q,s,u);
  f1 = u[n-1];
  if ((f0*f1) < 0)
    {
    bk = ek;
    dk = bk-ak;
    }
  else
    {
    ak = ek;
    dk = bk-ak;
    f0 = f1;
    }
  istep = istep+1;
  }
printf("%4d %16.8lf %16.8lf %16.8lf\n", istep,ek,dk,f1);
}

void nmrv (n,h,q,s,u)
/* The Numerov algorithm for the equation u"(x)+q(x)u(x)=s(x)
   as given in Eqs. (3.77)-(3.80) 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]-q[i+1])/2+q[i]);
  c1 = 2-g*(q[i+1]+q[i-1]+8*q[i]);
  c2 = 1+g*((q[i+1]-q[i-1])/2+q[i]);
  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;
  }
}