Time dependent Schrodinger equation in FreeFem

General Discussion
1

Good morning. I am trying to solve the time dependent Schrodinger equation using finite elements method. The problem I have is that since this is a complex equation I do not know how to define or work with complex functions and spaces in FreeFem++. At the beginning I got the compile error “Error: Problem a complex problem with no complex FE function“ and I think I solved it by setting Vh u, v, uu; I dont know if it is well done. By the way I used Cranck Nicolson Method in order to solve the evolution of time. Anyway I cannot get any plot of this code and I dont know what’s wrong. Thank you in advance

// Finite Elements Method: Time dependent Schrödinger equation in 2D
// Infinite Square Well

real L = 1.0;  

real hbar = 1.0545718e-34; 
real m    = 9.10938356e-31; 
real mu   = (hbar*hbar) / (2*m);

complex imaginary = 0. +1.0i;
int npoints = 30;    // puntos para la malla
real dt = 0.1;       // paso temporal

// ----- BOUNDARY DEFINITIONS -----
border lower(t=0,L) { x = t; y = 0; label=1; };
border right(t=0,L) { x = L; y = t; label=1; };
border upper(t=L,0) { x = t; y = L; label=1; };
border left(t=L,0) { x = 0; y = t; label=1; };

// ----- MESH -----
mesh Th = buildmesh(lower(npoints) + right(npoints) + upper(npoints) + left(npoints));
plot(Th, wait=1);

fespace Vh(Th, P1);       // espacio real
Vh<complex> u, v, uu;     // funciones complejas

problem schrodinger(u,v)
    = int2d(Th)(u*v + (imaginary*mu*dt/2.0)*(dx(u)*dx(v)+dy(u)*dy(v)))
    + int2d(Th)(-uu*v +(imaginary*mu*dt/2.0)*(dx(uu)*dx(v)+dy(uu)*dy(v)))
    + on(1, u=0);


real t=0;
real x0 = L/2;      // centro en x
real y0 = L/2;      // centro en y
real sigma = 0.1;   // ancho del paquete
real kx = 50.0;     // número de onda en x
real ky = 0.0;      // número de onda en y

uu = exp(-((x-x0)^2 + (y-y0)^2)/(2*sigma^2)) * exp(imaginary*(kx*x + ky*y));
real norm = sqrt(int2d(Th)(abs(uu)^2));
uu = uu / norm;

for(int m=0; m<=3/dt; m++){
	//Update
	t=t+dt;
	uu=u;
	//Solve
	schrodinger;

	//plot
	plot(abs(u), wait=true);   // módulo
	plot(real(u), wait=true);  // parte real
	plot(imag(u), wait=true);  // parte imaginaria
}