#include<stdio.h>
#include<gsl/gsl_matrix.h>
#include<gsl/gsl_eigen.h>
const double pi=3.14159265358979323844;
int main(void){
	int n=400;
	gsl_matrix *H = gsl_matrix_calloc(n,n);
	for(int i=0;i<n-1;i++){
//		gsl_matrix_set(H,i,i,-2);
		gsl_matrix_set(H,i,i,-3);
		gsl_matrix_set(H,i,i+1,1);
		gsl_matrix_set(H,i+1,i,1);
		}
//	gsl_matrix_set(H,n-1,n-1,-1);
//	gsl_matrix_set(H,n-1,n-1,-2);
	gsl_matrix_set(H,n-1,n-1,-3);
	gsl_matrix_scale(H,-1);
	gsl_eigen_symmv_workspace *w =  gsl_eigen_symmv_alloc (n);
	gsl_vector *eval = gsl_vector_alloc(n);
	gsl_matrix *evec = gsl_matrix_calloc(n,n);
	gsl_eigen_symmv(H,eval,evec,w);
	gsl_eigen_symmv_sort(eval,evec,GSL_EIGEN_SORT_VAL_ASC);

	for(int k=0;k<n/3;k++){// print first n/3 eigenvalues (rescaled)
		double scale=gsl_vector_get(eval,k)*(n+1)*(n+1)/pi/pi/(k+1)/(k+1);
		fprintf(stderr,"%g\n",scale);
		}
	
	for(int k=0;k<3;k++){// print first three eigenvectors
//		printf("%f %f\n",0.0,0.0); // psi_0 = 0
		for(int i=0;i<n;i++)
			printf("%f %f\n",(i+1.0)/(n+1.0),gsl_matrix_get(evec,i,k));
//		printf("%f %f\n",1.0,0.0); // psi_{n+1} = 0
		printf("\n\n");
		}

	gsl_matrix_free(H);
	gsl_matrix_free(evec);
	gsl_vector_free(eval);
	gsl_eigen_symmv_free(w);

return 0;
}