%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sqrwell.m
%
% Matlab Scrödinger Equation Solution for a free particle in an
% "square well" of size 2L (range [-L,L]) with V = +V0.
% A finite difference approach (cda) is used and all bound state
% energy eigenvalues are returned.
%
%     V0
% ----+     +----
%     | .   |
%     |     |
%     +-----+
%    -L     +L

% - Example usage for potential V0=50
%   >> sqrwell(50)
%
% - Note that following units are used in the calculation:
%   L = m = hbar = 1
%
% - May 2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function energy = sqrwell(potential)

 global L;
 global V0;
 L = 1.0;        % well width
 V0= potential;  % height of the well
 n = 1000;       % number of parts
 xL= -5*L;       % left  end
 xR= +5*L;       % right end
 dx= (xR-xL)/n;
 x = xL;
 A = zeros(n-1); % matrix elements
 for i=1:n-1
    for j=1:n-1
      if i==j+1 | i==j-1
          A(i,j) = -0.5/dx^2;
      end
    end
    x = x + dx;
    A(i,i) = 1.0/dx^2+V(x);
 end
 [wfn E] = eig(A); % get eigenvalues and eigenvectors

 x = xL+dx:dx:xR-dx;
 fprintf('*** Energy eigenvalues for bound states ***\n',n);
 for i=1:n-1
    if E(i,i)>=V0
        break;
    end
    energy(i) = E(i,i);
    wfn(:,i) = wfn(:,i)/max(wfn(:,i)); % normalize the wfn to 1
    fprintf('n=%2d --> E=%15.4f\n',i,energy(i));
    plot(x,wfn(:,i));
    axis([-3*L 3*L -2.0 2.0]);
    line([-L -L],[-2 2],[0 0],'Color','r')
    line([+L +L],[-2 2],[0 0],'Color','r')
    pause(1);
 end

end


% potential energy function definition
function y = V(x)
  global L;
  global V0;
  if abs(x)<=L
      y = 0.0;
  else
      y = V0;
  end
end
