clear all;

format long;


%  Initialize all of the parameters


i_imag = sqrt(-1);

N = 100;

L = 10000*10^-9;	% in meters

h = L/(N-1);

x = h*(0:N-1) - L/2;

h_bar = 1.055*10^-34;

m_gaas = .067*9.11*10^-31;

m_algaas = .09*9.11*10^-31;

tau = 1e-14;

V0 = .4*1.6*10^-19;


%  Set up the Hameltonian matrix

ham = zeros(N);

coeff = -h_bar^2/(2*m_gaas*h^2);

coeff2 = -h_bar^2/(2*m_algaas*h^2);

% Define the x-directed Hamiltonian

for ii=2:1:N-1
    
    hamx(ii,ii-1) = coeff;
    
    hamx(ii,ii)= -2*coeff;
    
    hamx(ii,ii+1) = coeff;
    
end

%  Define the first and last points

hamx(1,1) = -2*coeff+0;

hamx(1,2) = coeff;

hamx(N,N) = -2*coeff+0;

hamx(N,N-1) = coeff;

hamx(N,1)=coeff;

hamx(1,N)=coeff;

%  Compute the Crank-Nicholson matrix

dCNx =(inv(eye(N)+.5*i_imag*tau/h_bar*hamx)*(eye(N)-.5*i_imag*tau/h_bar*hamx));

%  Initialize the wavefunction  in the x-direction

x0 = 0;

velocity = 0;    	% in m/s

k0 = m_gaas*velocity/h_bar;

sigma0=L/10;

psi=exp(i_imag*k0*x').*exp(-(x'-x0).^2/(2*sigma0^2));

Norm_psi=sum(psi.*conj(psi));

psi = (1/sqrt(Norm_psi))*psi;

%  Plot the initial wavefunction

figure(1);

clf;

plot(x,real(psi),'-',x,imag(psi),'--');

title('Initial wave function');

xlabel('x');

ylabel('\psi(x)');

legend('Real','Imag');

drawnow;

pause(1);


%  Initialize loop and plot variables

max_iter=1e4;

plot_iter=max_iter/100;

p_plot(:,1)=psi.*conj(psi);

iplot=1;

figure(2);

clf;

axisV=[-L/2 L/2 0 max(p_plot)];

pmax = 0.5 * max(p_plot);

iter = 1;

cpmax = 1000;


%  Loop over the desired number of steps

while (cpmax>pmax)
   
   psi=dCNx*psi;
   
   %  Record plotting values
   
   if (rem(iter,plot_iter)<1)
      
      iplot=iplot+1;
      
      p_plot(:,iplot)=psi.*conj(psi);
      
      cpmax=max(p_plot);
      
      plot(x,p_plot(:,iplot));
      
      xlabel('x');
      
      ylabel('P(x,t)');
      
      title(sprintf('Finished %g of %g iterations',iter,max_iter));
      
      axis(axisV);
      
      drawnow;
      
   end
   
   iter = iter + 1;
   
end