#include <iostream>
#include "math.h"

#define ERROR_VAL -99999

double a0[1000];
double a[1000];

// NOTICE: We take hbar^2/2m = 1
int N = 100; // max power in the expansion
double v0 = 20; // depth of the well
double x0 = 1;  // half-width of the well

double bisect(double kl, double kr, double even);
double f(double e, double even);
double psi(double x);
double psi_0(double x,double e);
double dpsi(double x);
double dpsi_0(double x, double e);

void
calc_coef_even(double e)
{
  a[0] = 1;
  a[1] = 0;
  for(int i = 2; i < N; i++){
    a[i] = -(e+v0)*a[i-2]/i/(i-1);
  }
}

void
calc_coef_odd(double e)
{
  a[0] = 0;
  a[1] = 1;
  for(int i = 2; i < N; i++){
    a[i] = -(e+v0)*a[i-2]/i/(i-1);
  }
}

int
main()
{
  double pi = acos(-1);

  double e0;
  double step = 0.002;
  // even solution
  for(double e = -v0; e <= 0; e+=step){
    e0 = bisect(e,e+step,1.);
    if(e0 != ERROR_VAL){
      std::cout << "SOLUTION EVEN " << e0 << std::endl;
      for(double x = -20; x < -x0; x+=0.02)
        std::cout << x << " " << psi_0(-x,e0) << std::endl;
      for(double x = -x0; x < x0; x+=0.02)
        std::cout << x << " " << psi(x) << std::endl;
      for(double x = x0; x < 20; x+=0.02)
        std::cout << x << " " << psi_0(x,e0) << std::endl;
    }
  }
  // odd solution
  for(double e = -v0; e <= 0; e+=step){
    e0 = bisect(e,e+step,-1);
    if(e0 != ERROR_VAL){
      std::cout << "SOLUTION ODD " << e0 << std::endl;
      for(double x = -20; x < -x0; x+=0.02)
        std::cout << x << " " << -psi_0(-x,e0) << std::endl;
      for(double x = -x0; x < x0; x+=0.02)
        std::cout << x << " " << psi(x) << std::endl;
      for(double x = x0; x < 20; x+=0.02)
        std::cout << x << " " << psi_0(x,e0) << std::endl;
    }
  }

}

double f(double e, double even)
{
  if(even > 0.) calc_coef_even(e);
  else calc_coef_odd(e);

  return dpsi_0(x0,e)-dpsi(x0);
}

double psi(double x)
{
  double y = 0.;
  for(int i = 0; i < N; i++){
    y += a[i]*pow(x,double(i));
  }
  return y;
}

double psi_0(double x, double e)
{
  double k = sqrt(fabs(e));
  return exp(-k*x)/exp(-k*x0)*psi(x0);
}

double dpsi(double x)
{
  double y = 0.;
  for(int i = 1; i < N; i++){
    y += i*a[i]*pow(x,double(i-1));
  }
  return y;
}

double dpsi_0(double x, double e)
{
  double k = sqrt(fabs(e));
  return -k*exp(-k*x)/exp(-k*x0)*psi(x0);
}

double
bisect(double kl, double kr, double even)
{
  double ret_val = ERROR_VAL;
  double tol = 1.e-8;

  if(f(kl,even)*f(kr,even) > 0) return ret_val;

  while(true){
    double km = (kl+kr)/2.;
    double fm = f(km,even);
    double fl = f(kl,even);
    if(fm*fl < 0)
      kr = km;
    else
      kl = km;

//    std::cout << kl << " " << kr << " " << fm << " " << fl << " " << fabs(kl-kr) << std::endl;

    if(fabs(kl-kr) < tol) break;
  }

  ret_val = kl;
  return ret_val;
}