%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% sqrwell_matching.m
%
% An example application of the Matching 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 and potential V0=50
%   >> sqrwell_matching(1.1, 50)
%
% - Note that following units are used in the calculation:
%   L = m = hbar = 1
%
% - May 2013
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y = sqrwell_matching(energy_estimate, potential)

  % initilize some parameters
  global V0;
  global L;
  V0 = potential;
  L = 1.0;             % well width
  n = 1001;            % number of parts for divisiron of the well
  N = 3;               % for plotting axis [-N*L, N*L]
  dx= 2*N*L/(n - 1);   % dx steps
  xL= -N*L;            % left value of x
  xR= +N*L;            % right value of x
  xM= 0.0;             % matching point 
  jM= int32( n/2 );    % index of the matching point
  dE= 0.25;            % initial energy steps
  Eold = 0;            % old energy value
  E = energy_estimate; % new energy
  maxiter=999;         % maximum number of iterations
  diff0 = 0.0;         % derivative at the matching point
  % initilize the wave functions
  for i=1:n
    psiL(i) = 0.0;
    psiR(i) = 0.0;
  end
  % first derivatives at the left/right end points
  % the values can be (+,+), (+,-), (-,+) or (-,-) depending on parity
  psiL(2)  = +0.001*dx;
  psiR(n-1)= +0.001*dx;

  %--- main iteration loop -----------------------------------------
  for iter = 0:maxiter
      % fix the values of wave functions at end points
      psiL(1) = 0.0; psiL(n) = 0.0;
      psiR(1) = 0.0; psiR(n) = 0.0;
      % integrate the differential equations
      xLeft=xL; xRight=xR;
      for j=2:n-1 
         xLeft = xLeft + dx;
         psiL(j+1) = 2*psiL(j)*(1-dx^2*(E-V(xLeft)))-psiL(j-1);
         xRight = xRight-dx;
         psiR(n-j) = 2*psiR(n-j+1)*(1-dx^2*(E-V(xRight)))-psiR(n-j+2);
      end
      % normalize maximum value of wave functions to 1.0
      psiL = psiL/max(psiL);
      psiR = psiR/max(psiR);
      % Rescale the right hand side wave function such that values
      % of the wave functions at the matching point are the same
      psiR = psiR*abs(psiL(jM)/psiR(jM));
      % Difference between first derivatives at the matching point!
      diff1 =(psiL(jM+1)-psiL(jM))/(dx) - (psiR(jM+1)-psiR(jM))/(dx);
      % change energy
      if iter==1
          if diff1>0
             E = E + dE;
          else
             E = E - dE;
          end
          diff0 = diff1;
          Eold = E;
          continue
      else
          if diff1*diff0 < 0
             dE = dE/2;
          end
          if diff1>0
             E = E + dE;
          else
             E = E - dE;
          end
          diff0 = diff1;
      end
      % check if energy does not change so much
      if abs(Eold-E)<1.0e-5
          break
      end
      Eold = E;
        % Dynamic plot of wave functions for each iteration
        axis([-N*L N*L -2.0 2.0]);
        hold on;
        x = -N*L:dx:N*L;
        p1 = plot(x,psiL,'k:o');
        p2 = plot(x,psiR,'r');
        pause(0.1);
        delete(p1)
        delete(p2)
        hold off;
        fprintf ('iter = %d  E=%f\n',iter,E);
  end
  %--- end of main iteration loop ----------------------------------
 
  if iter==maxiter
     fprintf('*** CDA Algorithm does not converge: iter = maxiter ***\n');
  end

  % plot wave functions
  axis([-N*L N*L -2.0 2.0]);
  hold on;
  grid on;
  x = -N*L:dx:N*L;
  plot(x,psiL,'k:o');
  plot(x,psiR,'r');
  xlabel('x');
  ylabel('arbitrary units');
  title ('Mathing Method for Infinite Square Well');
  text(-0.9, 1.75,['V0 =  ',num2str(V0)],'BackgroundColor','y');
  text(+0.1, 1.75,'o psiLeft ','BackgroundColor','y');
  text(+0.1, 1.50,'- psiRight','BackgroundColor','y', 'Color','r');
  hold off;
  
  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
