%------------------------------------------------------------

% Schroedinger3

% Description: Similar to Schroedinger2; uses fixed-step RK4

%              to integrate the time-independent SE.

% Parameters (in order):

%   potential: Function handle to the potential function. Note

%               that the definition of V(x) sets the x-coords.

%   x:         Interval of integration, i.e. [-5 5]

%   psi0:      Value of wavefunction at x = 0;

%   d_psi0:    Wavefunction slope at x = 0;

%   E0:        The initial energy bracket to begin search.

%------------------------------------------------------------

function [x psi E] = schroedinger3(potential,x,psi0,d_psi0,E0)



% To avoid Matlab's warning

%------------------------------------------------------------

global V;

V = potential;



%------------------------------------------------------------

% Quality control: i.e. "knobs to turn"

%------------------------------------------------------------

global x_grid x_error x_interval e_grid e_error;

%iter_max = 200;



% Position (x) related:

%------------------------------------------------------------

x_grid     = 1000;    % Grid size

x_interval = x;       % Interval of integration



% Energy related:

%------------------------------------------------------------

e_grid  = 10;    % Grid size

e_error = 1E-20; % Maximum energy bracket width



E = find_energy(E0,psi0,d_psi0);



[x_r psi_r] = solve(0,x_interval(2),psi0,d_psi0,E);

[x_l psi_l] = solve(0,x_interval(1),psi0,d_psi0,E);



% Correctly combine the two directions

%------------------------------------------------------------

L     = length(x_l);

x_l   = fliplr(x_l);   x_l   = x_l([1:L-1]);

psi_l = fliplr(psi_l); psi_l = psi_l([1:L-1]);



x   = [x_l   x_r];

psi = [psi_l psi_r]; 



% Normalization of the wavefunction

%------------------------------------------------------------

I   = integral(x,psi.^2);

psi = psi/I;



% Results

%------------------------------------------------------------

plot(x,psi)       ; hold on;

plot(x,psi.^2,'r'); hold off;



%------------------------------------------------------------

% Function:    integral

% Description: Computes the integral of y over the given x

%              interval. (x and y should have same lengths)

%------------------------------------------------------------

function I = integral(x,y)

L  = length(x);

dx = x(2)-x(1);

I  = 0;



for k = 1:L

    I = I + y(k)*dx;

end



%------------------------------------------------------------

% Function: find_energy

% Desc: Given an energy "bracket" and the initial 

%       conditions, calculate the appropriate energy

%       that allows the wavefunction to converge at

%       the edges.

% Computational detail: Uses the RIGHT limit to test

%       for convergence. No promises about the left.

% Parameters:

%       E0:     Two-element vector that serves as a

%               bound for the search.

%       psi0:   Initial wavefunction

%       d_psi0: Initial wavefunction slope.

% Returns:

%       calcE:  Calculated energy

%------------------------------------------------------------

function calcE = find_energy(E0,psi0,d_psi0)



global e_grid e_error x_interval;

E         = linspace(E0(1),E0(2),e_grid);

converged = false;



while(~converged)

    endval = zeros(1,e_grid);

    for k = 1:e_grid

        [x psi]   = solve(0,x_interval(2),psi0,d_psi0,E(k));

        endval(k) = psi(length(psi));

    end

    

    % Now we zoom into a more specific energy interval

    k_crit = sign_change(endval);

    

    % Handle the error case

    if k_crit == -1

       error = sprintf('No convergent solution contained in the given energy interval.')

       break;

    end

    

    E = linspace(E(k_crit(1)),E(k_crit(2)),e_grid);

    if (E(2)-E(1)) < e_error

        converged = true;

    end

end



% Once the energy is small enough, return the energy of best convergence

k_crit = find_min(endval);

calcE  = E(k_crit);



%------------------------------------------------------------

% Function:    sign_change

% Description: Returns the index of a vector where its element

%              undergoes a sign change

% Error code:  (-1) is returned when no sign change occurs

%------------------------------------------------------------

function k_crit = sign_change(x)



s      = sign(x);

first  = s(1);

k_crit = -1;



for k = 2:length(x)

    if sign(s(k)) ~= first 

        k_crit = [k-1 k];

        break;

    end

end



%------------------------------------------------------------

% Function:    find_min

% Description: Returns the index of a vector where its element

%              assumes the smallest absolute value

%------------------------------------------------------------

function k_crit = find_min(x)



x = abs(x);

k_crit  = 1;



for k = 2:length(x)

    if x(k) < x(k_crit)

        k_crit  = k; 

    end

end



%------------------------------------------------------------

% Function: solve

% Desc: Solve Schroedinger equation

% Parameters:

%       x0:     Starting x

%       xF:     Terminal x

%       psi0:   Initial value of wavefunction

%       d_psi0: Initial slope of wavefunction

%       E:      Energy eigenvalue

% Returns:

%       psi:    Solution in row-vector form

%       x:      Corresponding values of x [x0:xF]

%------------------------------------------------------------

function [x psi] = solve(x0,xF,psi0,d_psi0,E)



global x_grid V;

x      = linspace(x0,xF,x_grid);

deltaX = x(2)-x(1);

psi    = zeros(x_grid,2);



% Input conditions

psi(1,1:2) = [psi0 d_psi0];



%------------------------------------------------------------

% RK4 integration: 

%  This iteration scheme is best visualized in

%  terms of vectors in (psi,d_psi) space.

%------------------------------------------------------------

for j = 2:x_grid

   eps  = E-V(x(j));   % Energy differential prefactor

   psi0 = psi(j-1,:);  % Current "position" in (psi,d_psi) space

   

   k1   = deltaX * [psi0(2) -eps*psi0(1)];

   psi1 = psi0 + 1/2*k1;

   

   k2   = deltaX * [psi1(2) -eps*psi1(1)];

   psi2 = psi0 + 1/2*k2;

   

   k3   = deltaX * [psi2(2) -eps*psi2(1)];

   psi3 = psi0 + k3;

   

   k4   = deltaX * [psi3(2) -eps*psi3(1)];

   

   psi(j,1:2) = psi(j-1,:) + 1/6*(k1 + 2*k2 + 2*k3 + k4);

end



% Return only the wavefunction, not its derivative

psi = psi(:,1).';