/* Equation de diffusion à d=1 */




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

#define sqr(x) (x)*(x)
#define delta(z) ( (fabs((z))<1e-8) ? 1 : 0 )
#define MIN(a,b) (a) > (b) ? (b) : (a)
#define N         // Nombre de points le long la direction spatiale 
#define L         // Longueur physique 
#define a (L)/(N) // spatial lattice spacing
#define M         // Nombre de pas de temps
#define T         // Temps physique total
#define h (T)/(M) // temporal lattice spacing
#define D        //   Coefficient de diffusion 
#define a2 sqr((a))
#define sx 


double width(double *X){

// A COMPLETER

}


double laplacian_1d(double *X,int n){//terme de diffusion
    double temp;
    temp=(X[n+1]-2.0*X[n]+X[n-1])/a2;
    return(temp);
    }


double hamiltonian(double *X, int n){
    double temp;
    temp=D*laplacian_1d(X,n);
    return(temp);
}

void init(double *X){
    int i;
    double z;
    for(i=1;i<N;i++){
	z=(double)i*a;
	X[i]=(double)exp(-(double)sqr(z-0.5*L)/(double)sx);
    }
}

void bc(double *X){
    X[0]=0.0;
    X[N]=0.0; 
}


void output(int i,double *X){
    int k;
    for(k=0;k<(N+1);k++){
	printf("%g %g %g\n",i*h,k*a,X[k]);
    }
}



void rk4(double *X){
    int k;
    double temp_X[N+1], temp_K[4][N+1];
    for(k=1;k<N;k++){
	temp_K[0][k]=h*hamiltonian(X,k);
	temp_X[k]=X[k]+0.5*temp_K[0][k];
    }
    for(k=1;k<N;k++){
	temp_K[1][k]=h*hamiltonian(temp_X,k);
    }
    for(k=1;k<N;k++){
	temp_X[k]=X[k]+0.5*temp_K[1][k];
    }
    for(k=1;k<N;k++){
	temp_K[2][k]=h*hamiltonian(temp_X,k);
    }
    for(k=1;k<N;k++){
	temp_X[k]=X[k]+temp_K[2][k];
    }
    for(k=1;k<N;k++){
	temp_K[3][k]=h*hamiltonian(temp_X,k);
    }
    for(k=1;k<N;k++){
	X[k]=X[k]+(1./6.)*(temp_K[0][k]+2.*(temp_K[1][k]+temp_K[2][k])+temp_K[3][k]);
    }
}


main(){
    double X[N+1];
    int k,i;

    init(X);
    bc(X);


    for(i=1;i<M;i++){
	rk4(X);

// A COMPLETER
		
	}
    }
}