% Shooting method example
% 5-7-2008
%
% See corresponding lecture for details on what is being solved here and
% what the algorithm's structure looks like.

clear all; close all;

TOLERANCE = 10^(-4);  % Tolerance for the shooting method below

n0 = 100;  % We are choosing a particular value for n0
col = ['r' 'g' 'b' 'k' 'y'];  % Different colors for each mode

A=1;  % This parameter actually stays constant in this implementation since
      % we are using beta to 'shoot' the boundary condition.  Usually this
      % would be the parameter we would vary
x0=[0 1];  % Initial condition
xspan=[-1 1];  % The span to solve for, in this case our domain [-1 1]

beta_start = n0;  % Our initial guess for beta each time
hold on;

for modes=1:5
    beta = beta_start;  % Set our initial guess to what this value
    dbeta = n0/100;  % Set the change in beta each time we guess wrong
    for j=1:1000
        [x,y] = ode45('shoot2',xspan,x0,[],n0,beta);  % Solve the ODE

        if abs(y(end,1)-0) < TOLERANCE  % Check to see if we have converged
            beta                        % and stop the shooting loop if we
            break;                      % have.
        end
        if (-1)^(modes+1)*y(end,1)>0    % Modify beta depending on how we
            beta=beta - dbeta;          % were off from the true solution.
        else                            % Note that we only shrink the       
            beta=beta + dbeta/2;        % amount we are changing beta if 
            dbeta=dbeta/2;              % the condition is not met
        end
    end

    beta_start = beta-0.1;  % Change the starting value of beta so we can 
                            % find the next eigenvalue

    norm = trapz(x,y(:,1).*y(:,1));  % In this example the true solution
                                     % needs to be normalized
    plot(x,y(:,1)/sqrt(norm),col(modes));  % Plot each mode we find
end
title('Eigen modes of the Schrodinger equation')
xlabel('x'); ylabel('\Psi_n')