// Résolution de l'équation de Schrödinger par une méthode spectrale
// On développe sur la base de fonctions propres de -d²/dx² = H_0
// Conventions: hbar = 1 et m = 0.5

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

#define sqr(x) ((x)*(x))
#define L 10.0
#define sigma (0.1*(L))
#define T 500.0 
#define PI acos(-1.0)

#define LPOINTS 100
#define TPOINTS 100

#define ax  (L)/(LPOINTS)
#define at  (T)/(TPOINTS)

double simpson( double (*f)(double), double A, double B, int N){
//Méthode de Simpson pour l'intégration numérique 
    double h=(B-A)/( (double)N +1.0);
    double temp=h*(f(A)+f(B))/3.0;
    int j;
    double x;
    for(j=1;j<N;j=j+2){
	x=A+j*h;
	temp=temp + (4.0*h*f(x)/3.0);}

    for(j=2;j<N;j=j+2){
	x=A+j*h;
	temp=temp + (2.0*h*f(x)/3.0);} 
	return(temp);
    }



double eigenvalue(int n){// Valeurs propres indexées par un entier n
    return(sqr(n*PI/L));}

double eigenfunction(int n, double x){// fonctions propres indexées par n
    return(sqrt(2.0/L)*sin((double)(n*PI*x/L)));}

double initialcondition(double x){// condition initiale
    //a compléter
  return( exp(-sqr(x-0.5*L)*0.5/sqr(sigma))/sqrt(2.0*PI*sqr(sigma)) );
}

double overlap(int n){/* calcule les coefficients du développement 
				   en fonctions propres de la solution */
  double f(double x){
    return(initialcondition(x)*eigenfunction(n,x));}
  
  return( simpson(f,0.0, L, 1025) ); // varier le nombre de points intermédiaires (ici 1025)
}

main(){
    complex somme;
    int n, i, j;
    double x=0.5*L,t=T*0.1, cn;
    int NMAX, N; 


    scanf("%d", &NMAX);
    
    

    for(i=0;i<=LPOINTS;i++){
	x = i*ax;
	for(j=0;j<=TPOINTS;j=j+(TPOINTS/10)){
	    t = j*at;
//    for(N=2;N<NMAX; N = 2*N){
	    somme = 0.0;
	    for(n=1;n<NMAX;n++){
		cn = overlap(n);
		somme = somme + cn*eigenfunction(n, x)*cexp(-I*eigenvalue(n)*t);}

	    //printf("%d %g\n",N, cabs(somme));
	    //}
	
	    	printf("%g\t%g\t%g\n",t,x,cabs(somme));
	}
    }
}