function R=qmCN(N,M,T,x,ic,V);
% QMCN: R=qmCN(N,M,T,x,ic,V)
%    Solves the Schrodinger equation
%       i u_t = -u_xx + V(x) u,   u(x,0)=ic
%    using Crank--Nicholson finite diff method on N x-subintervals.
%    Note x is an N+1 length row vector of equally--spaced 
%    x-values and dt=T/M.  Allows initial condition specified in ic.
%    Requires size(x)==size(ic)==size(V)==[1 N+1].
%    Always uses zero Dirichlet boundary conditions (which is 
%         reasonable for confining potentials like V=x^2).
% Try "N=100; dx=20/N; x=-10:dx:10; IC=pi^(-.25)*exp(i*-5*x-(x-2).^2/2);
%      qmCN(N,100,5,x,IC,zeros(1,N+1));"
%    for a gaussian wave packet traveling freely leftward.
% Try "qmCN(N,100,5,x,IC,x.^2);"
%    for gaussian wave packet in harmonic oscillator.
% (ELB 5/5/02)

if not(isequal(size(x),[1 N+1])), error('x wrong dimension'); end;
if not(all(diff(diff(x))<1.0e-12)), error('x values not equally spaced.'); end;
if not(isequal(size(x),size(ic))), error('init cond wrong dimension'); end;
if not(isequal(size(x),size(V))), error('V wrong matrix dimension'); end;
dx=x(2)-x(1); dt=T/M; t=0:dt:T; R=dt/(dx)^2;

w=zeros(M+1,N+1); b=zeros(1,N-1); w(1,:)=ic;
c=.5*R*ones(N-1,1); % sub and super diagonal
d=i-R-.5*dt*V(2:N)'; % diagonal
%should work: A=spdiags([c d c],-1:1,N-1,N-1);

for k=2:M+1
   w(k,1)=0; w(k,N+1)=0; % boundary conditions
   for j=2:N
      b(j-1)=-.5*R*w(k-1,j-1)+(i+R+.5*dt*V(j))*w(k-1,j)-.5*R*w(k-1,j+1); end;
   w(k,2:N)=tridiag(c,d,c,b,N-1);
   %for some reason "w(k,2:N)=(A\ b')';" does not work.
end;

prob=w.*conj(w); PP=max(max(prob));
mesh(x,t,prob);
axis([x(1) x(N+1) t(1) t(M+1) 0 PP]);
xlabel('x axis'); ylabel('t axis'); zlabel('|u|^2 prob. axis');
title(['QMCN: Squared magnitude of solution of Schrodinger '...
      'equation with N = ' num2str(N)]);
hold on;
AA=axis;
plot3(x,T*ones(1,N+1),(V/(max(abs(V+.001))))*(AA(6)-AA(5)));
text(x(1), T, (AA(6)-AA(5))*.7, 'Graph of V (not to scale).');
plot3(x(N+1)*ones(1,M+1),t,5*sum(prob,2)/N);
text(x(N+1), T/2, (AA(6)-AA(5))*.3, 'Estimate of norm of |u|^2.');
hold off;