schro.m

%  schro - Program to solve the Schrodinger equation 
%  using Crank-Nicolson scheme (Free particle version)
clear;  help schro;   % Clear memory and print header
i_imag = sqrt(-1);    % Imaginary i
N = input('Enter 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(,) %%%%%
ham = zeros(N);  % Set all elements to zero
coeff = -h_bar^2/(2*mass*h^2);
for i=2:(N-1)
  ham(i,i-1) = coeff;
  ham(i,i) = -2*coeff;  % Set interior rows
  ham(i,i+1) = coeff;
end
% Periodic boundary conditions
ham(1,N) = coeff;   ham(1,1) = -2*coeff; ham(1,2) = coeff;
ham(N,N-1) = coeff; ham(N,N) = -2*coeff; ham(N,1) = coeff;
disp('Computing the Crank-Nicolson matrix')
dCN = ( inv(eye(N) + .5*i_imag*tau/h_bar*ham) * ...
             (eye(N) - .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));
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);    % The particle should circle system
plot_iter = max_iter/20;        % Produce 20 curves
p_plot(:,1) = psi.*conj(psi);   % Record initial condition
iplot = 2;
%%%%% MAIN LOOP %%%%%
disp('Entering main loop')
for iter=1:max_iter
  psi = dCN*psi;  % Crank-Nicolson scheme
  if( rem(iter,plot_iter) < 1 )  % Every plot_iter steps record 
    p_plot(:,iplot) = psi.*conj(psi); % amplitude of psi(i)
    plot(x,p_plot(:,iplot));
    xlabel('Position'); ylabel('Prob. density');
    title(sprintf('Finished %g of %g iterations',iter,max_iter));
    drawnow;
    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');