#include <stdio.h>
#include <complex.h>
#include <string.h>

#define N 100
#define NP 100

double V(double x){
return (x<5.31&&x>4.99)?50.:0.;
//return 0.;
}

int main(){
double dx=10./N,dt=.0001;
complex psi[N],fs[N]={0.},fs1[N]={0.};
double M=1.-dt*dt/dx/dx/dx/dx;
int ii,i;
FILE *fp;

fp=fopen("sch.dat","r");
fread(psi,sizeof(complex),N,fp);
fclose(fp);

for(ii=0;ii<NP;ii++){
  for(i=1;i<N-1;i++) fs1[i]=
          (I+dt/dx/dx)*(fs[i-1]+fs[i+1])*dt/2./dx/dx/M
          +2.*cexp(-I*V(i*dx)*dt)*(1-I*dt/dx/dx)*psi[i]/M;
  memcpy(fs,fs1,N*sizeof(complex));
  }
for(i=0;i<N;i++){
  psi[i]=fs[i]-cexp(-I*V(i*dx)*dt)*psi[i];
  printf("%lf %lf %lf\n",i*dx,creal(psi[i]),cimag(psi[i]));
  }

fp=fopen("sch.dat","w");
fwrite(psi,sizeof(complex),N,fp);
fclose(fp);

}
