%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sqrwell_shooting.m
%
% An example application of the Shooting Method to solve Schrödinger 
% Equation for an electron in an infinite well of size 2L (range [-L,L]).
% A finite difference approach (cda) is used for the integration
% of the Schrödinger equation. Energy eigenvalue is returned.
%
%     V0
% ----+     +----
%     | .   |
%     |     |
%     +-----+
%    -L     +L

% - Example usage for energy estimate E=1.1, potential V0=50 and parity P=1
%   >> sqrwell_shooting(1.1, 50, 1)
%
% - parity can be 0 (even parity) or 1 (odd parity)
%
% - Note that following units are used in the calculation:
%   L = m = hbar = 1
%
% - May 2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = sqrwell_shooting(energy_estimate, potential, parity)

  % initilize some parameters
  global V0;
  global L;
  V0 = potential;
  L = 1.0;        % well width
  n = 1000;       % number of parts for divisiron of the well
  N = 3;          % for plotting axis [-N*L, N*L]
  dx= N*L/n;      % dx steps
 
  dE = 0.5;
  Enew = energy_estimate;
  E = Enew + dE;
  psi_old = 0;
  psi = 0;
  X = 0:dx:N*L-dx;
  W = zeros(n,1);
  while 1
      x=0; psi=1-parity; z=parity;
      for i=1:n
          z = z + 2*(V(x)-E)*psi*dx;
          psi = psi + z*dx;
          x = x + dx;
          W(i) = psi;
      end
      p = plot(X,W,'.');
      if psi*psi_old<0
          dE = -dE/2;
      end
      E = E + dE;
      psi_old = psi;
      if abs(E-Enew)<eps
           break
      end
      Enew = E;            
      pause(0.01);
      delete(p);
      fprintf('Parity = %d and Energy = %f\n',parity, E);
  end
  y = E;
end

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

