///////////////////////////////////////////////////////////////////////////
//                                                                       //
// Program file name: Sturm.java                                         //
//                                                                       //
// © Tao Pang 2006                                                       //
//                                                                       //
// Last modified: January 18, 2006                                       //
//                                                                       //
// (1) This Java program is part of the book, "An Introduction to        //
//     Computational Physics, 2nd Edition," written by Tao Pang and      //
//     published by Cambridge University Press on January 19, 2006.      //
//                                                                       //
// (2) No warranties, express or implied, are made for this program.     //
//                                                                       //
///////////////////////////////////////////////////////////////////////////

// An example of implementing the simplest method to
// solve the Sturm-Liouville problem.

import java.lang.*;
public class Sturm {
  static final int n = 100, ni = 10;
  public static void main(String argv[]) {
    double del = 1e-6, l = 0.5, dl = 0.1;
    l = secant(ni, del, l, dl);

 // Output the eigenvalue obtained
    System.out.println("The eigenvalue is: " + l);
  }

// Method to carry out the secant search.

  public static double secant(int n, double del,
    double x, double dx) {
    int k = 0;
    double x1 = x+dx;
    while ((Math.abs(dx)>del) && (k<n)) {
      double d = f(x1)-f(x);
      double x2 = x1-f(x1)*(x1-x)/d;
      x  = x1;
      x1 = x2;
      dx = x1-x;
      k++;
    }
    if (k==n) System.out.println("Convergence not" +
      " found after " + n + " iterations");
    return x1;
  }

// Method to provide the function for the root search.

  public static double f(double l) {
    double u[] = new double[n+1];
    double p[] = new double[n+1];
    double q[] = new double[n+1];
    double s[] = new double[n+1];
    double p1[] = new double[n+1];
    double h = 1.0/n;
    double u0 = 0;
    double u1 = h;

    for (int i=0; i<=n; ++i){
      double x = h*i;
      p[i] = 1-x*x;
      p1[i] = -2*x;
      q[i]  = l*(l+1);
      s[i]  = 0;
    }
    u = sturmLiouville(h, u0, u1, p, p1, q, s);
    return u[n]-1; 
  }

// Method to integrate the Sturm-Liouville problem.

  public static double[] sturmLiouville(double h,
    double u0, double u1, double p[], double p1[],
    double q[], double s[]) {
    int n = p.length-1;
    double u[] = new double[n+1];
    double h2 =  h*h;
    u[0] = u0;
    u[1] = u1;
    for (int i=1; i<n; ++i){
      double c2 = 2*p[i]+h*p1[i];
      double c1 = 4*p[i]-2*h2*q[i];
      double c0 = 2*p[i]-h*p1[i];
      double d  = 2*h2*s[i];
      u[i+1] = (c1*u[i]-c0*u[i-1]+d)/c2;
    }
    return u;
  }
}