/*
 * Schrodinger equation for infinite well
 */

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

#define N 80

/* 
 * 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 d[2*N+1], e[2*N+1], xstep;
  int i, j, m, info;
  /* Compute eigenvalues only */
  /*
  char jobz[1] = {'N'};          
  double *z = (double *) NULL;
  double *work = (double *) NULL;
  */
  /* Compute eigenvalues and eigenvectors */
  char jobz[1] = {'V'};          
  double z[(2*N+1)*(2*N+1)], work[4*N];

  m = 2*N + 1;
  xstep = 1./(N+1);

  for (i = 0; i < m; i++) {
    d[i] = 2.;
    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++) {
      printf("%5i    %14.6f    %14.6f\n", i, d[i]/(xstep*xstep*M_PI*M_PI), (i+1)*(i+1)/4.);
    }
    for (j = 0; j < m; j++ ) {
      printf("%5i    ", j);
      for (i = 0; i < 5; i++ ) {
        /*printf("%14.6f   ", z[i+j*m]);*/
        printf("%14.6f   ", z[i*m+j]);
      }
      printf("\n");
    }
  }
  
  return(0);
}