/*
 * Schrodinger equation for infinite well
 */

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

#define N 100
#define XMAX 10.

/* 
 * Computes all eigenvalues, and optionally, eigenvectors of a real
 * symmetric tridiagonal matrix.
 */
void dstev_(char *jobz, int *m1, double *d, double *e, double *z, int *m2, double *work, int *info);

int main() {
  double xstep, xstep2, xi;
  int i, j, m=2*N+1, info;
  char jobz[1] = {'V'};   /* Compute eigenvalues and eigenvectors */       
  double d[m], e[m];
  double z[m*m], work[4*m];

  xstep = XMAX/N;
  xstep2 = xstep*xstep;

  for (i = 0; i < m; i++) {
    xi = (i - N)*xstep;
    d[i] = 2. + xi*xi*xstep2;
    e[i] = -1.;
  }

  dstev_ (jobz, &m, d, e, z, &m, work, &info);

  if (info) {
    printf("Problems finding eigenvalues ...\n");
  } else {
    for (i = 0; i < 10; i++) {
      fprintf(stderr, "%5i    %14.6f\n", i, d[i]/(2*xstep2));
    }
    for (j = 0; j < m; j++ ) {
      printf("%5i    ", j);
      for (i = 0; i < 5; i++ ) {
        fprintf(stdout, "%14.6f   ", z[i*m+j]);
      }
      printf("\n");
    }
  }
  
  return(0);
}