/*************************    Program 3.4    ****************************/
/*                                                                      */
/************************************************************************/
/* 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

main()
/* Main program for solving the Legendre equation with
   the simplest algorithm for the Sturm-Liouville
   equation and the bisection method for the root search.
   Copyright (c) Tao Pang 1997. */
{
int i,n,istep;
double dl=1e-6;
double h,ak,bk,ek,dk,f0,f1;
double u[NMAX];
void smpl();

/* Initialization of the problem */

n = NMAX;
h = 2.0/(n-1);
ak = 0.5;
bk = 1.5;
dk = 0.5;
istep = 0;
u[0] = -1;
u[1] = -1+h;
ek = ak;
smpl (n,h,ek,u);
f0 = u[n-1]-1;

/* Bisection method for the root */

while (fabs(dk) > dl)
  {
  ek = (ak+bk)/2;
  smpl (n,h,ek,u);
  f1 = u[n-1]-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 smpl (n,h,ek,u)
/* The simplest algorithm for the Sturm-Liouville equation.
   Copyright (c) Tao Pang 1997. */
int n;
double h,ek;
double u[];
{
int i;
double h2,q,p,p1,x;

h2 = 2*h*h;
q = ek*(1+ek);
for (i = 1; i <n-1; ++i)
  {
  x = -1+i*h;
  p = 2*(1-x*x);
  p1 = -2*x*h;
  u[i+1] = ((2*p-h2*q)*u[i]+(p1-p)*u[i-1])/(p1+p);
  }
}