function superposition

T     = 500;
t     = linspace(0,50,T);
x_lim = [-5 5];

[x y1 E1] = schroedinger3(@potential,x_lim,1,0,[0,5]);
[x y2 E2] = schroedinger3(@potential,x_lim,0,1,[0 5]);
L = length(x);

c1 = 1/sqrt(2);
c2 = 1/sqrt(2);

phase = 0;

for k = 1:T
    for j = 1:L
        y(j) = c1*y1(j)*exp(-i*E1*t(k)) + c2*y2(j)*exp(-i*(E2*t(k)+phase)); 
    end
    plot(x,real(y)); axis([x_lim -1 1]); hold on;
    plot(x,y.*conj(y),'r'); hold off;
    drawnow;
end