// http://gcc.gnu.org/onlinedocs/gcc-4.0.4/gcc/Complex.html#Complex
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <complex.h>
#define span 80
typedef _Complex double c;
double const pi = 3.1415926535897932385;
static double sq(double x){return x*x;}
c im(c x){return I*creal(x)-cimag(x);}
double dt = .00005, dx = .05;
c psi[2*span+1];
static void pc(int j, c x) {
   printf("%9.3f %13.4f %13.4fi\n", (j-span)*dx, creal(x), cimag(x));}

int main(){double time = 0, dpt = .05, npt = -dt/2;
double sig = 1/sqrt(2); // standard deviation of position of particle.
void pt(char * m){printf(" time = %e, %s\n", time, m);}
void pv(char* m){int i; pt(m);
{int j; double s=0;
     for(j=1; j<2*span-1; ++j) s += sq(creal(psi[j])) + sq(cimag(psi[j]));
     printf("sum = %16.12e\n", dx*s);}
   for(i=0; i<2*span+1; ++i) {
   if(0)printf("%3d", i); pc(i, psi[i]);}}
{// Initialize psi.
  double f = 1./(sqrt(sqrt(2*pi)*sig)), scl=dx/(2*sig);
  int j=2*span+1; while(j--) psi[j] = f*exp(-sq((j-span)*scl));}
while(time<.2+dt) {c hold=0; pt("cx");
  if(time > npt) {npt += dpt; pv("during");}
  {int j; for(j=1; j<2*span; ++j) {
      c new = psi[j]+im((dt/sq(dx))*(2*psi[j]-(psi[j-1]+psi[j+1])));
      if((j > span-2 && j < span+1) ||j==span+40) pc(j, psi[j]);
      psi[j-1] = hold; hold = new;} 
    psi[2*span-1] = hold;
    time += dt;}
  psi[0] = psi[2*span-1]; psi[2*span] = psi[1];
  }
pt("end");
return 0;}

/*(dx, dt, t)
(.1, .001. 0)
    0.000        0.7511        0.0000i
    0.100        0.7474        0.0000i
    0.700        0.5879        0.0000i
    1.600        0.2088        0.0000i
    1.700        0.1771        0.0000i
(.1, .001. .04)
    0.000        0.7494        0.0299i
    0.100        0.7457        0.0294i
    0.700        0.5878        0.0120i
    1.600        0.2098       -0.0130i
    1.700        0.1779       -0.0133i
(.05, .0001. .04)
    0.000        0.7493        0.0299i
    0.100        0.7456        0.0295i
    0.700        0.5878        0.0120i
    1.600        0.2098       -0.0130i
    1.700        0.1779       -0.0134i
(.1, .001. .1)
    0.000        0.7403        0.0732i
    0.100        0.7369        0.0721i
    0.700        0.5869        0.0303i
    1.600        0.2148       -0.0318i
    1.700        0.1823       -0.0330i
(.05, .0001. .1)
    0.000        0.7452        0.0726i
    0.100        0.7419        0.0718i
    0.700        0.5947        0.0353i
    1.600        0.2104       -0.0304i
    1.700        0.1789       -0.0342i
(.1, .001. .2)
    0.000        0.6914        0.1570i
    0.100        0.7275        0.1152i
    0.700        0.5866        0.0530i
    1.600        0.2280       -0.0468i
    1.700        0.2047       -0.0737i
(.05, .0001. .2)
    0.000     -729.4081     2158.0605i
    0.100     -681.0391     2063.9470i
    0.700       31.6971     1855.9067i
    1.600     -555.1055     -825.4354i
    1.700     -600.0964    -1017.5029i
(.05, .00005. .2)
    0.000        0.6981        0.1150i
    0.100        0.6977        0.1142i
    0.700        0.5872        0.0583i
    1.600        0.2346       -0.0660i
    1.700        0.2075       -0.0712i
As if there were a requirement that dt < dx2
*/