#include <iostream>
#include <cmath>
#include <cstdlib>
using namespace std;

// Solve the analytic constraint equations resulting from the
// finite-well problem for a specific value of U0.

#define real	double
#define TOL	1.e-4

real U0;				// scaled exterior potential
int  even;

// We want to solve xi sin xi = eta cos xi (even), or -xi cos xi = eta
// sin xi (odd), subject to xi^2 + etai^2 = U0.

real g(real z)				// equation to solve, z is xi
{
    real eta = U0 - z*z;
    if (eta < 0) eta = 0;		// avoid problems with z^2 ~ U0
    eta = sqrt(eta);

    if (even)
	return z*sin(z) - eta*cos(z);
    else
	return z*cos(z) + eta*sin(z);
}

real solve(real z1, real z2, real epsilon)
{
    // Use bisection until the range is less than the specified
    // tolerance.

    real zm;
    while (z2 - z1 > epsilon) {
	zm = (z1 + z2)/2;
	if (g(zm)*g(z1) > 0)
	    z1 = zm;
	else				// note no check for equality
	    z2 = zm;
    }

    // Perform one final false-position iteration.

    return z1 + (z2 - z1) * (-g(z1)) / (g(z2) - g(z1));
}

main(int argc, char** argv)
{
    real zl, zr, dz;
    const char* oddeven[2] = {"odd", "even"};

    U0 = 4;				// scaled potential

    if (argc > 1)
	U0 = atof(argv[1]);		// specify U0 on the command line

    dz = 0.01*sqrt(U0);
    if (dz > 0.25) dz = 0.25;

    for (even = 1; even >= 0; even--) {

	// Systematically seek new roots between 0 and U0.

	zl = 0;

	while (zl < sqrt(U0)-TOL*dz) {
	    zr = zl + dz;
	    if (zr > sqrt(U0)) zr = sqrt(U0);

	    if (g(zl)*g(zr) <= 0 && (even || zl > 0)) {
		real z = solve(zl, zr, TOL);
		printf("Found %s energy eigenvalue %f (xi = %f)\n",
			oddeven[even], z*z, z);
	    }

	    zl = zr;
	}
    }
}