% QM PIB simulation w/ shooting method for E

% set parity
odd = false;
even = true;
parity = odd;
% parity = even;
fac = 1/2;

% initial guess for E
E = pi^2 * fac;
dE = 0.5;
divergeSign = 1; % arbitrary, will change later

% constants
L = 1;
A = 1;
Vo = 1e3;
dx = 0.01;
x = (0:dx:L*1.5)';

% trap divergence here
PsiLim = 2;

% IC & pre-allocate
if parity
    Psi = A*ones(length(x)+1,1);  % note extra place to left
else
    Psi = [-dx ; zeros(length(x),1)]; % note perturbed LHS so that we
                                      % end up w/ pure zeros in shoot
end  

while 1  % use CRTL+c to break out of this loop
    lastSign = divergeSign;

    % make a run through
    for ii = 3:length(x)+1
        if x(ii-1) <= L
            V = 0;
        else
            V = Vo;
        end
        Psi(ii) = QM_PIBstep(V,E,Psi(ii-2),Psi(ii-1),dx);
        if abs(Psi(ii)) > PsiLim
            break;
        end
    end

    divergeSign = sign(Psi(ii));

    QM_PIBplot(Psi,x,A,L,parity)
    title(['E = ' rats(E/pi^2) '\pi^2; {\Delta}E = ' num2str(dE)])
    pause

    if lastSign ~= divergeSign
        dE = -dE/2;
    end
    E = E + dE;

end

return