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

double v[40000];
double x[40000];

int N = 20000;
double v0 = 10.;
double x0 = 1;
double h;
double dl;
double dr;
double tol = 1.e-8;

double bisect(double kl, double kr, double even);
double f(double e, double even);
double psi_r(double e, double &der);
double psi_l(double e, double &der);
void write_psi(double e, double even);

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

  h = 8*x0/N;
  for(int i = 0; i <= N; i++){
    x[i] = -4*x0+i*h;
    v[i] = 0.;
    if(fabs(x[i]) < x0) v[i] = -v0+v0/x0*fabs(x[i]);
  } 


  double e0;
  double step = 0.002;
  // even solution
  for(double e = -v0; e <= 0; e+=step){
    e0 = bisect(e,e+step,1.);
    if(e0 != 99999 && dl*dr > 0.){
      std::cout << "SOLUTION EVEN " << e0 << std::endl;
      write_psi(e0,1.);
    }
  }
  // odd solution
  for(double e = -v0; e <= 0; e+=step){
    e0 = bisect(e,e+step,-1);
    if(e0 != 99999 && dl*dr < 0.){
      std::cout << "SOLUTION ODD " << e0 << std::endl;
      write_psi(e0,-1);
    }
  }

}

double f(double e, double even)
{
  return (psi_l(e,dl)-even*psi_r(e,dr));
//  double yl = psi_l(e,dl);
//  double yr = psi_r(e,dr);
//  return (dl-factor*dr);
}

void write_psi(double e, double even)
{
  double k = sqrt(fabs(e));
  double hh = h*h;
  double c = hh/12.;
  double y0 = exp(k*x[0]);
  double y1 = exp(k*x[1]);
  double w0 = (1-c*(v[0]-e))*y0; 
  double w1 = (1-c*(v[1]-e))*y1; 
  double w2;
  for(int i = 2; i <= 3*N/8; i++){
/*
    w2 = 2*(1-5./12.*hh*(e-v[i-1]))*y1-(1+c*(e-v[i-2]))*y0;
    w2 = w2/(1+c*(e-v[i]));
    y0 = y1;
    y1 = w2; 
*/
    w2 = -w0+2*w1+hh*(v[i-1]-e)*y1;
    w0 = w1;
    w1 = w2;
    y1 =  w2/(1-c*(v[i]-e));
    std::cout << x[i] << " " << y1 << std::endl;
  }

  y0 = exp(k*x[0]);
  y1 = exp(k*x[1]);
  w0 = (1-c*(v[0]-e))*y0;
  w1 = (1-c*(v[1]-e))*y1;
  for(int i = N-2; i >= 3*N/8; i--){
/*
    w2 = 2*(1-5./12.*hh*(e-v[i+1]))*y1-(1+c*(e-v[i+2]))*y0;
    w2 = w2/(1+c*(e-v[i]));
    y0 = y1;
    y1 = w2; 
*/
    w2 = -w0+2*w1+hh*(v[i+1]-e)*y1;
    w0 = w1;
    w1 = w2;
    y1 =  w2/(1-c*(v[i]-e));
    std::cout << x[i] << " " << even*y1 << std::endl;
  }

}

double psi_l(double e, double &der)
{
  double k = sqrt(fabs(e));
  double hh = h*h;
  double c = hh/12.;
  double y0 = exp(k*x[0]);
  double y1 = exp(k*x[1]);
  double w0 = (1-c*(v[0]-e))*y0; 
  double w1 = (1-c*(v[1]-e))*y1; 
  double w2;
  for(int i = 2; i <= 3*N/8; i++){
    w2 = -w0+2*w1+hh*(v[i-1]-e)*y1;
    w0 = w1;
    w1 = w2;
    y1 =  w2/(1-c*(v[i]-e));
  }
  y0 =  w0/(1-c*(v[N/2-1]-e));
  der = (y1-y0)/h;
  return y1;
}

double psi_r(double e, double &der)
{
  double k = sqrt(fabs(e));
  double hh = h*h;
  double c = hh/12.;
  double y0 = exp(k*x[0]);
  double y1 = exp(k*x[1]);
  double w0 = (1-c*(v[N]-e))*y0; 
  double w1 = (1-c*(v[N-1]-e))*y1; 
  double w2;
  for(int i = N-2; i >= 3*N/8; i--){
    w2 = -w0+2*w1+hh*(v[i+1]-e)*y1;
    w0 = w1;
    w1 = w2;
    y1 =  w2/(1-c*(v[i]-e));
  }
  y0 =  w0/(1-h*h/12.*(v[N/2+1]-e));
  der = (y0-y1)/h;
  return y1;
}

double
bisect(double kl, double kr, double even)
{
  double ret_val = 99999;

  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;
}