schrot.m

%  Program to solve the Schrodinger equation using sparce matrix
%  Crank-Nicolson scheme (Particle-in-a-box version)
clear;  help schrot;  % Clear memory and print header
i_imag = sqrt(-1);
N = input('Enter the number of grid points - ');
L = 100;                % System extends from -L/2 to L/2
h = L/(N-1);            % Grid size
h_bar = 1;  mass = 1;   % Natural units
tau = input('input time step - ');
%%%%%  Set up the Hamiltonian operator matrix ham(,) %%%%%
coeff = -h_bar^2/(2*mass*h^2);
for i=2:(N-1)
  ham(i,1) = coeff;
  ham(i,2) = -2*coeff;
  ham(i,3) = coeff;
end
% Direchlet boundary conditions
ham(1,1)=0; ham(1,2)=0; ham(1,3)=0;
ham(N,1)=0; ham(N,2)=0; ham(N,3)=0;
%%%%% Set up the matrix Q %%%%%
tri_eye = zeros(N,3);
tri_eye(1:N,2) = ones(N,1);   % Identity matrix in packed format
Q = 0.5*( tri_eye + 0.5*i_imag*tau/h_bar * ham );
%%%%% Initialize wavefunction %%%%%
x0 = 0;         % Location of the center of the wavepacket
velocity = 0.5; % Average velocity of the packet
k0 = mass*velocity/h_bar;       % Average wavenumber
sigma0 = L/10;  % Standard deviation of the wavefunction
Norm_psi = 1/(sqrt(sigma0*sqrt(pi)));  % Normalization
x = h*(0:N-1) - L/2;
psi = Norm_psi * exp(i_imag*k0*x') .* ...
                      exp(-(x'-x0).^2/(2*sigma0^2));
% Set phi to zero as boundary condition, i.e. perfect reflector
psi(1) = 0;  psi(N)=0;
plot(x,real(psi),x,imag(psi));
title('Hit any key to continue');
xlabel('Position');  ylabel('Real(psi) and imag(psi)');
pause;
%%%%% Set loop and plot variables %%%%%
max_iter = L/(velocity*tau);    % Particle should circle the system
plot_iter = max_iter/50;  % Produce 50 curves
p_plot(:,1) = psi.*conj(psi);  % Record initial condition
iplot = 2;
%%%%% MAIN LOOP %%%%%
disp('Entering main loop')
for iter=1:max_iter
  chi = tri_ge(Q,psi);  % Solve for vector chi
  psi = chi - psi;      % Crank-Nicolson scheme
  if( rem(iter,plot_iter) < 1 )  % Every plot_iter steps record 
    p_plot(:,iplot) = psi.*conj(psi);  % psi(i) for ploting
    plot(x,p_plot(:,iplot));
    xlabel('Position'); ylabel('Prob. density');
    title(sprintf('Finished %g of %g iterations\n',iter,max_iter));
    drawnow;
    tplot(iplot) = iter*tau;
    iplot = iplot+1;
  end
end
plot(x,p_plot(:,1:3:iplot-2),x,p_plot(:,iplot-1));
xlabel('Position'); ylabel('Prob. density');
title('Probability density at various times');
pause;
mesh(tplot,x,p_plot);
view([90 30]);
ylabel('Position');  xlabel('Time');