#include <stdio.h>
#include <math.h>

#define L     50.0
#define imx   50
#define bE    0.0
#define dx    1.0
#define dE    1e-7
#define e     1e-3
#define buf   100

#define _odd
#define even

static double fL[imx+1];
static double fR[imx+1];
static char   fn[buf+1];
static int    fg[imx+1];

double V(int);

int main()
{
  FILE *fp;
  double kp,E = bE;
  int i,j;
  
  do{
    for(j=0;j<=imx;j++)
        {
          fL[j] = fR[j] = 0.0;
          fg[j] = 0;
        }
      
    fL[1]     =  1.0;
#ifdef odd
    fR[imx-1] = -1.0; 
#endif
#ifdef even
    fR[imx-1] =  1.0; 
#endif

    for(j=1;j<imx-1;j++)
      {
        kp = E - V(j);
        fL[j+1] = (2.0 - kp*dx*dx)*fL[j] - fL[j-1];
      }
    for(j=imx-1;j>1;j--)
      {
        kp = E - V(j);
        fR[j-1] = (2.0 - kp*dx*dx)*fR[j] - fR[j+1];
      }
    
    for(j=0;j<=imx;j++)
      if(fabs(fL[j]-fR[j]) <= e)
        fg[j] = 1;
      
    i = 0;
    for(j=0;j<=imx;j++)
      i += fg[j];
    
    if(i == imx+1) break; else E += dE;
  }while(1);
  
  printf("Energy = %f\n",E);
  
  sprintf(fn,"result.dat");
  fp = fopen(fn,"w");
  if(fp == NULL) return -1;
  
  for(j=0;j<=imx;j++)
    fprintf(fp,"%f %f %f %f\n",j*dx,fL[j],fR[j],V(j));
  fclose(fp);
  
  return 0;
}

double V(int i)
{
  /* return -sin(3.141592*(double)(i*dx/L)); */
  return 0.0;
}